# 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)