# questo non e' un vero script da eseguire, ma una serie di
# comandi da eseguire in successione

library(coda)

molla.mcmc <- read.coda("bugs.out","bugs.ind",quiet=TRUE)
summary(molla.mcmc)
plot(molla.mcmc)

molla.m <- as.matrix(molla.mcmc)

mean(molla.m[,7])
sd(molla.m[,7])
mean(molla.m[,8])
sd(molla.m[,8])

cor(molla.m[,1],molla.m[,2])
cor(molla.m[,4],molla.m[,5])
cor(molla.m[,7],molla.m[,8])

# etc, etc: per altri risultati vedi il log di bugs

#Inferenza su k condizionata da g=9.8:
length(molla.m[,7][abs(molla.m[,8]-9.8)<0.05])  # 4819
mean(molla.m[,7][abs(molla.m[,8]-9.8)<0.05])    # 43.78099
sd(molla.m[,7][abs(molla.m[,8]-9.8)<0.05])      # 0.4387989
# usando una tolleranza di 0.01 il risultato di mean e sd non cambia

# nota: le distribuzioni di k e g non sono `molto normali'
# my.skew e my.kurt sono definite in lin_fit.R
my.skew(molla.m[,7])    #  0.3400837
my.kurt(molla.m[,7])    #  5.549563
my.skew(molla.m[,8])    #  5.308687
my.kurt(molla.m[,8])    #  0.3493722

hist(molla.m[,7],nc=100,freq=F)
curve(dnorm(x,mean(molla.m[,7]),sd(molla.m[,7])),35,50,add=TRUE)
hist(molla.m[,8],nc=100,freq=F)
curve(dnorm(x,mean(molla.m[,8]),sd(molla.m[,8])),8,11,add=TRUE)
hist(molla.m[,2],nc=100,freq=F)
curve(dnorm(x,mean(molla.m[,2]),sd(molla.m[,2])),0.21,0.24,add=TRUE)
hist(molla.m[,3],nc=100,xlim=c(0,0.0035),freq=F)
hist(molla.m[,6],nc=100,xlim=c(0,0.02),freq=F)