## simulazione Monte Carlo di un albero di eventi 

p.E <- rbeta(10000,5,30)
p.O <- rbeta(10000,5,90)
p.N <- p.E*12/15 + p.O*3/15
hist(p.N,30)

mean(p.N)
sd(p.N)

## soglia al 95%

p.N.ord <- sort(p.N)
p.N.ord[9500]

## cambio di media da prior a posterior

mean(rbeta(10000,5,30))
mean(rbeta(10000,12,33))

## come i momenti dipendono da k

n=100
x1 <- rnorm(n,0,2)
x2 <- rnorm((n-2),0,2)
a<- 10
x2 <- c(x2,a,-a)

m1 <- 1:10
m2 <- 1:10
for (k in 1:10) {
m1[k] <- ((1/n)*sum((x1-mean(x1))^(2*k)))^(1/(2*k))
m2[k] <- ((1/n)*sum((x2-mean(x2))^(2*k)))^(1/(2*k))
}
plot(m2)
lines(m1,col="red")

## si osservi che il primo valore, la sd, è simile.