Nothing
Gibbs.ANOVA <- function(data,it=5000,burnin=500,seed=0){
#R program for Gibbs sampling from conditionals
if(seed!=0){set.seed(seed)}
I=it+burnin #number of iterations
#read only observations with complete information
x=data$g
y=data$y
N=length(y)
G=length(unique(data$g))
#other file input
fit.lm <- lm(y~as.factor(x)-1) #lm, no intercept (!)
x <- model.matrix(fit.lm)[,,drop = FALSE]
#establish parameter vectors and constant quantities
s1=matrix(1,I); b1=matrix(0,I,G)
s2=matrix(1,I); b2=matrix(0,I,G)
xtxi=solve(t(x)%*%x)
pars=coefficients(lm(y~x-1))
#Gibbs sampling begins
for(i in 2:I){
#simulate beta from its multivariate normal conditional
b1[i,]=pars+t(rnorm(G,mean=0,sd=1))%*%chol((s1[i-1]^2)*xtxi) #choleski decomposition
#simulate sigma from its inverse gamma distribution
s1[i]=sqrt(1/rgamma(1,N/2,.5*t(y-x%*%(b1[i,]))%*%(y-x%*%(b1[i,]))))
}
for(i in 2:I){
#simulate beta from its multivariate normal conditional
b2[i,]=pars+t(rnorm(G,mean=0,sd=1))%*%chol((s2[i-1]^2)*xtxi) #choleski decomposition
#simulate sigma from its inverse gamma distribution
s2[i]=sqrt(1/rgamma(1,N/2,.5*t(y-x%*%(b2[i,]))%*%(y-x%*%(b2[i,]))))
}
## plot posterior density for parameters of interest with indications for credible interval
par1=cbind(b1,s1)[-c(1:burnin),]
par2=cbind(b2,s2)[-c(1:burnin),]
par=list(par1,par2)
posterior <- rbind(par1,par2)
colnames(posterior) <- paste("Mean",1:(G+1))
colnames(posterior)[G+1] <- "SD"
#save and restore user par
oldpar <- par(no.readonly = TRUE)
on.exit(par(oldpar))
#set par for output plots
par(mfrow = c(G+1, 2), mar = rep(2, 4))
for (j in 1:(G+1)) {
hist(posterior[,j], breaks="Scott", main=colnames(par)[j], xlab="Value")
abline(v=quantile(posterior[,j], probs=c(0.025, 0.975)), col="blue")
plot.ts(cbind(par1[,j],par2[,j]),col=c("blue","green"), plot.type="single",ylab="")
}
results <- cbind(round(apply(posterior,2,mean),4),round(apply(posterior,2,median),4),
round(apply(posterior,2,sd),4))
colnames(results) <- c("mean","median","sd")
output <- list("summary"=results,"posterior"=posterior)
#return(output)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.