# Rottakoe. Esimerkki hierarkisesta mallista # Aineisto y<-c(rep(0,14),rep(1,8),rep(2,9),1,5,2,5,3,2,7,7,3,3,2,9,10,rep(4,7),10,4,4,4,5,11,12,5,5,6,5,6,6,6,6, 16,15,15,9,4) n<-c(rep(20,7),rep(19,4),18,18,17,rep(20,4),19,19,18,18,25,24,23,rep(20,6),10,49,19,46,27,17,49,47,20, 20,13,48,50,rep(20,7),48,19,19,19,22,46,49,20,20,23,19,22,20,20,20,52,47,46,24,14) # Logaritmoitu posteriorijakauma lposterior<-function(par,n,y) { J<-length(n) alpha<-exp(sum(par))/(1+exp(par[1])) beta<-exp(par[2])/(1+exp(par[1])) lprior<-(-5/2)*log(alpha+beta)+log(alpha)+log(beta) lusk <-J*(lgamma(alpha+beta)-lgamma(alpha)-lgamma(beta))+ sum(lgamma(alpha+y)+lgamma(beta+n-y)-lgamma(alpha+beta+n)) lprior+lusk } # Muodostetaan diskretoitu jakauma par1<-seq(-2.2,-1.4,len=200) par2<-seq(1,5,len=200) jakauma<-matrix(0,nr=200,nc=200) for(i in 1:200) for(j in 1:200) jakauma[i,j]<-lposterior(c(par1[i],par2[j]),n,y) jakauma<-exp(jakauma-max(jakauma)) # Korkeuskäyräkuvio contour(par1,par2,jakauma,levels=seq(0.05,0.95,by=0.1)*max(jakauma),drawlabels=FALSE,xlab="log(alpha/beta)",ylab="log(alpha+beta)") # Simuloinnit par1.tf<-apply(jakauma,1,sum) index <- sample(1:200,1000,replace=TRUE,prob=par1.tf) par1.sim <- par1[index] par2.sim<-numeric(1000) for(i in 1:1000) par2.sim[i] <- sample(par2,1,prob=jakauma[index[i],]) plot(par1.sim,par2.sim,xlim=c(-2.2,-1.4),ylim=c(1,5),xlab="log(alpha/beta)",ylab="log(alpha+beta)") alpha<-exp(par1.sim+par2.sim)/(1+exp(par1.sim)) beta<-exp(par2.sim)/(1+exp(par1.sim)) # Populaatioparametrien simuloinnit theta<-matrix(0,nr=71,nc=1000) for(i in 1:1000) for(j in 1:71) theta[j,i]<-rbeta(1,alpha[i]+y[j],beta[i]+n[j]-y[j]) matplot(y/n,theta,cex=0.1,col="black",xlim=c(0,0.4),ylim=c(0,0.4),xlab=expression(y[j]/n[j]),ylab=expression(theta[j])) lines(c(0,0.4),c(0,0.4)) matpoints(y/n,apply(theta,1,median),cex=4,pch=".")