I apologize if this posting shows up again, for some reason I wasn't able to
post from a different account. So, here I am (reborn). Could I have some
suggestions as to how I can display my results and their respective CIs in an
aesthetically pleasing manner?
Below is the example code.
rm(list = ls())
set.seed(1)
func <- function(d,t,beta,lambda,alpha,p.gamma,delta,B){
d <- c(5,1,5,14,3,19,1,1,4,22)
t <- c(94.32,15.72,62.88,125.76,5.24,31.44,1.048,1.048,2.096,10.48)
post <- matrix(0, nrow = 11, ncol = B)
theta <- c(lambda,beta)
beta.hat <- 2.471546
for(j in 1:B){
for(i in 1:(B-1)){
c.lambda <- rgamma(10,d+alpha,t+beta.hat)
c.beta <- rgamma(1,10*alpha+p.gamma,delta+sum(lambda))
c.theta <- c(c.lambda,c.beta)
pi.func <-
prod((c.lambda/lambda)^(d+alpha-1)*exp(-t*(c.lambda-lambda)))*(c.beta/beta)^(10*alpha+p.gamma-1)*exp(-c.beta*(delta+sum(c.lambda))+beta*(delta+sum(lambda)))
g.x <-
prod((beta.hat+t)^(alpha+d)*lambda^(alpha+d-1)*exp(-lambda*(beta.hat+t))/gamma(alpha+d))*(sum(lambda)+delta)^(p.gamma+10*alpha)*beta^(p.gamma+10*alpha-1)*exp(-beta*(sum(lambda)+delta))/gamma(p.gamma+10*alpha)
g.y <-
prod((beta.hat+t)^(alpha+d)*c.lambda^(alpha+d-1)*exp(-lambda*(beta.hat+t))/gamma(alpha+d))*(sum(c.lambda)+delta)^(p.gamma+10*alpha)*c.beta^(p.gamma+10*alpha-1)*exp(-c.beta*(sum(c.lambda)+delta))/gamma(p.gamma+10*alpha)
a <- pi.func*(g.x/g.y)
if(a>1){
theta<-c.theta
}
else
theta <- theta+(c.theta-theta)*rbinom(1,1,alpha)
}
post[,j] <- theta
#print(post[,j])
}
mean <- apply(post,1,mean)
ci <- apply(post,1,quants)
return(list("The Means" = mean, "The Corresponding Confidence
Intervals" = ci))
}
quants <- function(x){
lo <- quantile(x,.025)
hi <- quantile(x,.975)
return(c(lo,hi))
}
(out <-
func(d,t,beta=1.5,lambda=rep(0.5,10),alpha=1.8,p.gamma=0.01,delta=1,B=100))
