MCMC_attempt.R
In jackiemauro/hurdleIV:
# another attempt at mcmc
#### laplace demon ####
#http://www.sumsar.net/blog/2013/06/three-ways-to-run-bayesian-models-in-r/
rm(list = ls())
library(LaplacesDemon)
y = rnorm(1000,4,2)
# The model specification
model <- function(parm, data) {
mu <- parm[1]
sigma <- exp(parm[2])
log_lik <- sum(dnorm(data$y, mu, sigma, log = T))
post_lik <- log_lik + dnorm(mu, 0, 100, log = T) + dlnorm(sigma, 0, 4, log = T)
# This list is returned and has to follow a the format specified by
# LaplacesDemon.
list(LP = post_lik, Dev = -2 * log_lik, Monitor = c(post_lik, sigma), yhat = NA,
parm = parm)
}
# Running the model
data_list <- list(N = length(y), y = y, mon.names = c("post_lik", "sigma"),
parm.names = c("mu", "log.sigma"))
mcmc_samples <- LaplacesDemon(Model = model, Data = data_list, Iterations = 30000,
Algorithm = "HARM", Thinning = 1)
plot(mcmc_samples, BurnIn = 10000, data_list)
mcmc_samples$Summary1['mu','Mean'];mcmc_samples$Summary1['mu','SD']
mcmc_samples$Summary1['sigma','Mean'];mcmc_samples$Summary1['sigma','SD']
# basic IV
rm(list = ls())
require(mvtnorm); require(MASS)
N = 10000
z <- rnorm(N,2,1)
sig2y1<-2; sig2x<-3; t0<-.2
Pi1<-1; Pi2<-3; B1<-2; B2<-1
errs = mvrnorm(n = N, mu = c(0,0), Sigma = matrix(c(sig2y1,t0,t0,sig2x),ncol = 2, byrow = T))
x2 <- Pi1 + Pi2*z + errs[,2]
y1 <- B1 + B2*x2 + errs[,1]
model <- function(parm, data) {
b1 <- parm[1]; b2 <- parm[2]; p1 <- parm[3]; p2 <- parm[4]
sigy1 <- exp(parm[5]); sigx <- exp(parm[6]); t0 <- parm[7]
mu_x <- p1 + p2*data$z
mu_y <- b1 + b2*mu_x
Sig = matrix(c(sigy1,t0,t0,sigx),ncol = 2, byrow = T)
#Transformation matrix from (eta, u, v) to (y0*, log y1*, x2) (without the mean)
A = rbind(c(1,b2),c(0,1))
#Covariance of (y0*, log y1*, x2)
Sig2 = A%*%Sig%*%t(A)
mu_y_x = mu_y + (Sig2[1,2]/Sig2[2,2])*(data$x-mu_x)
sig2_y_x = Sig2[1,1] - (Sig2[1,2]^2)/Sig2[2,2]
if(sig2_y_x>0){
log_lik <- sum(dnorm(data$y,mu_y_x,sqrt(sig2_y_x),log = T) + dnorm(data$x,mu_x,sqrt(Sig2[2,2]),log =T))
}
else{
mat = data.frame(y=data$y,x=data$x,predy=mu_y,predx = mu_x)
singlelikelihoods = apply(mat,1,function(x) dmvnorm(x[1:2],mean = x[3:4],sigma = Sig,log = T))
log_lik <- sum(singlelikelihoods)
}
log_prior <- sum(dnorm(c(b1,b2,p1,p2,t0),0,100,log = T))+dlnorm(sigy1, 0, 4, log = T)+dlnorm(sigx, 0, 4, log = T)
post_lik <- log_lik + log_prior
# This list is returned and has to follow a the format specified by
# LaplacesDemon.
list(LP = post_lik, Dev = -2 * log_lik, Monitor = c(post_lik, sigx,sigy1), yhat = NA,parm = parm)
}
# Running the model
data_list <- list(N = length(y1), y = y1,x = x2, z = z , mon.names = c("post_lik", "sigx","sigy1"),
parm.names = c("B1","B2","Pi1","Pi2","log.sigy","log.sigx","tau0"))
mcmc_samples <- LaplacesDemon(Model = model, Data = data_list, Iterations = 50000,
Algorithm = "HARM", Thinning = 1)
plot(mcmc_samples, BurnIn = 10000, data_list)
parameters = c('B1','B2','Pi1','Pi2','tau0','sigx','sigy1')
true = c(B1,B2,Pi1,Pi2,t0,sig2x,sig2y1)
names(true)<-parameters
for(par in parameters){
print(true[par] )
print(mcmc_samples$Summary1[par,'Mean'])
print(mcmc_samples$Summary1[par,'SD'])
}
# lognormal IV
rm(list = ls())
N = 10000
z <- rnorm(N,2,1); x1 <- rnorm(N,-1,3)
sig2y1<-2; sig2y0<-1; rho <- .1; t0<-.2; t1 <- .3; sig2x<-3
Pi1<-1; Pi2<-3; B1<-2; B2<-1;G1<- -1; G2 <- 2
errs = mvrnorm(n = N, mu = c(0,0,0), Sigma = matrix(c(sig2y0,rho,t0,rho,sig2y1,t1,t0,t1,sig2x),ncol = 3, byrow = T))
x2 <- Pi1*x1 + Pi2*z + errs[,3]
y0Star <- G1*x1 + G2*z + errs[,2]
lny1Star <- B1*x1 + B2*x2 + errs[,1]
y1 <- as.numeric(y0Star>0)*exp(lny1Star)
model <- function(parm, data) {
sigx <- exp(parm[1]); sigy1 <- exp(parm[2]); t1 <- parm[3]; t0 <- parm[4]
b1 <- parm[5]; b2 <- parm[6]; g1 <- parm[7]; g2 <- parm[8]; p1<-parm[9];p2<-parm[10];rho <- parm[11]
log_lik <- loglik_lgnorm_rho(parm)
log_prior <- sum(dnorm(c(b1,b2,g1,g2,p1,p2,t1,t0,rho),0,100,log = T))+dlnorm(sigy1, 0, 4, log = T)+dlnorm(sigx, 0, 4, log = T)
post_lik <- log_lik + log_prior
print(post_lik)
# This list is returned and has to follow a the format specified by
# LaplacesDemon.
list(LP = post_lik, Dev = -2 * log_lik, Monitor = c(post_lik, sigx,sigy1), yhat = NA,parm = parm)
}
# Running the model
data_list <- list(N = length(y1), y = y1,x = x2, x1=x1,z = z , mon.names = c("post_lik", "sigx","sigy1"),
parm.names = c("log.sigx","log.sigy","tau1","tau0","B1","B2","G1","G2","Pi1","Pi2","rho"))
mcmc_samples <- LaplacesDemon(Model = model, Data = data_list, Iterations = 100000,
Algorithm = "HARM", Thinning = 1,
Initial.Values = c(1,1,0,0,B1,B2,G1,G2,Pi1,Pi2,rho))
plot(mcmc_samples, BurnIn = 10000, data_list)
True.vals <- c(sig2x,sig2y1,t1,t0,B1,B2,G1,G2,Pi1,Pi2,rho)
names <- c('sigx','sigy1','tau1','tau0','B1','B2','G1','G2','Pi1','Pi2','rho')
Compare <- cbind(mcmc_samples$Summary2[names,c("Mean","SD")],True.vals)
# run it within the hurdleIVMCMC code
rm(list = ls())
z <- rnorm(N,2,1); x1 <- rnorm(N,-1,3)
sig2y1<-2; sig2y0<-1; rho <- .1; t0<-.2; t1 <- .3; sig2x<-3
Pi1<-1; Pi2<-3; B1<-2; B2<-1;G1<- -1; G2 <- 2
bs = c(0,B1,B2); gs = c(0,G1,G2); pis = c(0, Pi1, Pi2)
sig2y1<-2; sig2y0<-1; rho <- .1; t0<-.2; t1 <- .3; sig2x<-3
inst_mean = 2; inst_sd = 1; exog_mean = -1; exog_sd = 3
nms = c('rho','sigy1','tau0','tau1','sigx','beta1','beta2','beta3','gamma1','gamma2','gamma3','pi1','pi2','pi3')
true.vals = c(rho,sig2y1,t0,t1,sig2x,bs,gs,pis)
dat = hurdle.IV.sim(pi = pis, gamma = gs, beta = bs,
exog_mean = exog_mean, exog_sd = exog_sd,
z_mean = inst_mean, z_sd = inst_sd,
endog_sd = sqrt(sig2x), y_sd = sqrt(sig2y1),
rho = rho, tau0 = t0, tau1 = t1, n = 10000)
out = hurdle.IV.MCMC(y~exog1 + endog,
inst = inst1,
endog = endog,
exog = exog1,
data = dat,
k = 50000)
Compare <- cbind(out$Summary2[nms,c("Mean","SD")],true.vals)
# ##### this is from: https://nicercode.github.io/guides/mcmc/ ####
# just to get mean
m = 0
f <- function(x) dnorm(x,m,1)
q <- function(x) rnorm(1,x,1)
step <- function(x,f,q){
xp = q(x) #propose a new x
alpha = min(1,f(xp)/f(x)) #calculate alpha
if(runif(1)<alpha){x<-xp}
x
}
run <- function(x, f, q, nsteps) {
res <- matrix(NA, nsteps, length(x))
for (i in seq_len(nsteps))
res[i,] <- x <- step(x, f, q)
drop(res)
}
res <- run(-10, f, q, 1000)
layout(matrix(c(1, 2), 1, 2), widths=c(4, 1))
par(mar=c(4.1, .5, .5, .5), oma=c(0, 4.1, 0, 0))
plot(res, type="s", xpd=NA, ylab="Parameter", xlab="Sample", las=1)
usr <- par("usr")
xx <- seq(usr[3], usr[4], length=301)
plot(f(xx), xx, type="l", yaxs="i", axes=FALSE, xlab="")
#### looking at these guys' example: https://theoreticalecology.wordpress.com/2010/09/17/metropolis-hastings-mcmc-in-r/####
rm(list = ls())
trueA <- 5
trueB <- 0
trueSd <- 10
sampleSize <- 31
# create independent x-values
x <- (-(sampleSize-1)/2):((sampleSize-1)/2)
# create dependent values according to ax + b + N(0,sd)
y <- trueA * x + trueB + rnorm(n=sampleSize,mean=0,sd=trueSd)
likelihood <- function(param){
a = param[1]
b = param[2]
sd = param[3]
pred = a*x + b
singlelikelihoods = dnorm(y, mean = pred, sd = sd, log = T)
sumll = sum(singlelikelihoods)
return(sumll)
}
prior <- function(param){
a = param[1]
b = param[2]
sd = param[3]
aprior = dunif(a, min=0, max=10, log = T)
bprior = dnorm(b, sd = 5, log = T)
sdprior = dunif(sd, min=0, max=30, log = T)
return(aprior+bprior+sdprior)
}
posterior <- function(param){
return (likelihood(param) + prior(param)) #being bayesian
#return(likelihood(param))#what if i don't be bayesian?
}
proposalfunction <- function(param){
return(rnorm(3,mean = param, sd= c(0.1,0.5,0.3)))
}
run_metropolis_MCMC <- function(startvalue, iterations){
chain = array(dim = c(iterations+1,3))
chain[1,] = startvalue
for (i in 1:iterations){
proposal = proposalfunction(chain[i,])
probab = exp(posterior(proposal) - posterior(chain[i,]))
if (runif(1) < probab){
chain[i+1,] = proposal
}else{
chain[i+1,] = chain[i,]
}
}
return(chain)
}
startvalue = c(4,0,10)
chain1 = run_metropolis_MCMC(startvalue, 10000) #bayesian
chain2 = run_metropolis_MCMC(startvalue, 10000) #not bayesian (change posterior)
burnIn = 5000
acceptance = 1-mean(duplicated(chain1[-(1:burnIn),]))
names = c('a','b','sig'); true = c(trueA,trueB,trueSd)
for(i in 1:dim(chain1)[2]){
hist(chain1[-(1:burnIn),i],nclass=30, main=paste("Posterior of ",names[i]), xlab="True value = red line" )
abline(v = mean(chain1[-(1:burnIn),i]));abline(v = true[i], col="red" )
}
#### basic IV #####
# even though dmvnorm and dnorm using conditionals give the same probabilities
# asking the algorithm to estimate the mean of x/y within the likelihood makes it
# perform less well
rm(list = ls())
require(mvtnorm)
N = 10000
z <- rnorm(N,2,1)
sig2y1<-2; sig2x<-3; t0<-.2
Pi1<-1; Pi2<-3; B1<-2; B2<-1
errs = mvrnorm(n = N, mu = c(0,0), Sigma = matrix(c(sig2y1,t0,t0,sig2x),ncol = 2, byrow = T))
x <- Pi1 + Pi2*z + errs[,2]
y <- B1 + B2*x + errs[,1]
likelihood <- function(param){
s1 = param[1];s2 = param[2]; t = param[3]; p1 = param[4]; p2 = param[5]; b1 = param[6]; b2 = param[7]
xpred = p1 + p2*z
ypred = b1 + b2*xpred
preds = matrix(c(ypred,xpred),ncol = 2, byrow = F)
obs = matrix(c(y,x),ncol = 2, byrow = F)
Sig = matrix(c(s1,t,t,s2),ncol = 2, byrow = T)
# These two will give the same likelihoods, except, I think, when sig2_y1_x2<0
# Sig_err = matrix(c(s1,t,t,s2),ncol = 2, byrow = T)
# A = matrix(c(1,b1,0,1),ncol = 2, byrow = T)
# Sig = A%*%Sig_err%*%t(A)
#
# mu_y1_x2 = ypred + Sig[1,2]/Sig[2,2]*(x-xpred)
# sig2_y1_x2 = Sig[1,1] - Sig[1,2]^2/Sig[2,2]
#
# ylik = sum(dnorm(y,mean = mu_y1_x2,sd = sqrt(sig2_y1_x2),log = T))
# xlik = sum(dnorm(x,mean = xpred,sd = sqrt(Sig[2,2]),log = T))
# sumll1 = xlik + ylik
#
mat = data.frame(y=y,x=x,predy=ypred,predx = xpred)
singlelikelihoods = apply(mat,1,function(x) dmvnorm(x[1:2],mean = x[3:4],sigma = Sig,log = T))
sumll = sum(singlelikelihoods)
#
# print(sumll)
# print(sumll1)
return(sumll)
}
likelihood <- function(param){
# this one has to guess at the true value of x in order to get the conditional
# mean of y--faster but appears not to be as accurate.
s1 = param[1];s2 = param[2]; t = param[3]; p1 = param[4]; p2 = param[5]; b1 = param[6]; b2 = param[7]
xpred = p1 + p2*z
ypred = b1 + b2*xpred
preds = matrix(c(xpred,ypred),ncol = 2, byrow = F)
obs = matrix(c(x,y),ncol = 2, byrow = F)
#Sig = matrix(c(s1,t,t,s2),ncol = 2, byrow = T)
Sig_err = matrix(c(s1,t,t,s2),ncol = 2, byrow = T)
A = matrix(c(1,b1,0,1),ncol = 2, byrow = T)
Sig = A%*%Sig_err%*%t(A)
mu_y1_x2 = ypred + Sig[1,2]/Sig[2,2]*(x-xpred)
sig2_y1_x2 = Sig[1,1] - Sig[1,2]^2/Sig[2,2]
if(sig2_y1_x2<=0 | Sig[2,2]<=0){return(-Inf)}
ylik = sum(dnorm(y,mean = mu_y1_x2,sd = sqrt(sig2_y1_x2),log = T))
xlik = sum(dnorm(x,mean = xpred,sd = sqrt(Sig[2,2]),log = T))
sumll = xlik + ylik
return(sumll)
}
prior <- function(param){
s1 = param[1];s2 = param[2]; t = param[3]; p1 = param[4]; p2 = param[5]; b1 = param[6]; b2 = param[7]
s1prior = dunif(s1, min=0, max=10, log = T)
s2prior = dunif(s2, min=0, max=10, log = T)
tprior = dnorm(t, sd = 5, log = T)
p1prior = dnorm(p1, sd = 5, log = T)
p2prior = dnorm(p1, sd = 5, log = T)
b1prior = dnorm(p1, sd = 5, log = T)
b2prior = dnorm(p1, sd = 5, log = T)
return(s1prior + s2prior + tprior + p1prior + p2prior + b1prior + b2prior)
}
posterior <- function(param){
return (likelihood(param) + prior(param)) #being bayesian
#return(likelihood(param))#what if i don't be bayesian?
}
proposalfunction <- function(param){
return(rnorm(7,mean = param, sd= c(0.1,0.5,0.3,.5,.5,.5,.5)))
}
run_metropolis_MCMC <- function(startvalue, iterations){
chain = array(dim = c(iterations+1,length(startvalue)))
chain[1,] = startvalue
for (i in 1:iterations){
proposal = proposalfunction(chain[i,])
probab = exp(posterior(proposal) - posterior(chain[i,]))
if (!is.na(probab) & (runif(1) < probab)){
chain[i+1,] = proposal
}else{
chain[i+1,] = chain[i,]
}
}
return(chain)
}
names = c('sig2y1','sig2x','tau','pi1','pi2','b1','b2')
true = c(sig2y1,sig2x,t0,Pi1,Pi2,B1,B2)
startvalue = c(rep(1,length(true)))
chain1 = run_metropolis_MCMC(startvalue, 1000)
chain2 = run_metropolis_MCMC(startvalue, 1000)
burnIn = 5000
acceptance = 1-mean(duplicated(chain1[-(1:burnIn),]))
for(i in 1:dim(chain1)[2]){
hist(chain1[,i],nclass=30, main=paste("Posterior of ",names[i]), xlab="True value = red line" )
abline(v = mean(chain1[,i]));abline(v = true[i], col="red" )
}
xmins1 = apply(chain1,2,min);xmins2 = apply(chain2,2,min);xmins = apply(cbind(xmins1,xmins2),1,min)
xmaxs1 = apply(chain1,2,max);xmaxs2 = apply(chain2,2,max);xmaxs = apply(cbind(xmaxs1,xmaxs2),1,max)
par(mfrow = c(1,2))
for(i in 1:dim(chain2)[2]){
hist(chain1[,i],nclass=30, main=paste("Posterior of ",names[i]," chain 1"), xlim = c(xmins[i],xmaxs[i]),xlab="True value = red line" )
abline(v = mean(chain1[,i]));abline(v = true[i], col="red" )
hist(chain2[,i],nclass=30, main=paste("Posterior of ",names[i], "chain 2"), xlim = c(xmins[i],xmaxs[i]), xlab="True value = red line" )
abline(v = mean(chain2[,i]));abline(v = true[i], col="red" )
}
par(mfrow = c(1,1))
#### with fake loglik lgnorm ####
rm(list = ls())
require(mvtnorm)
N = 10000
z <- rnorm(N,2,1)
sig2y1<-2; sig2y0<-1; rho <- .1; t0<-.2; t1 <- .3; sig2x<-3
Pi1<-1; Pi2<-3; B1<-2; B2<-1;G1<- -1; G2 <- 2
errs = mvrnorm(n = N, mu = c(0,0,0), Sigma = matrix(c(sig2y0,rho,t0,rho,sig2y1,t1,t0,t1,sig2x),ncol = 3, byrow = T))
x <- Pi1 + Pi2*z + errs[,3]
y0Star <- G1 + G2*z + errs[,2]
y1Star <- B1 + B2*x + errs[,1]
y1 <- as.numeric(y0Star>0)*y1Star
likelihood <- function(param){ #not done
sy0 = param[1]; sy1 = param[2]; sx2 = param[3];
r = param[4]; t0 = param[5]; t1 = param[6]
p1 = param[7]; p2 = param[8];g1 = param[9];g2 = param[10]; b1 = param[11]; b2 = param[12]
xpred = p1 + p2*z
y0pred = g1 + g2*xpred
y1pred = b1 + b2*xpred
preds = matrix(c(ypred,xpred),ncol = 2, byrow = F)
obs = matrix(c(y,x),ncol = 2, byrow = F)
Sig = matrix(c(sy0,r,t0,r,sy1,t1,t0,t1,sx2),ncol = 3, byrow = T)
mat = data.frame(y=y,x=x,predy=ypred,predx = xpred)
singlelikelihoods = apply(mat,1,function(x) dmvnorm(x[1:2],mean = x[3:4],sigma = Sig,log = T))
sumll = sum(singlelikelihoods)
return(sumll)
}
prior <- function(param){
s1 = param[1];s2 = param[2]; t = param[3]; p1 = param[4]; p2 = param[5]; b1 = param[6]; b2 = param[7]
s1prior = dunif(s1, min=0, max=10, log = T)
s2prior = dunif(s2, min=0, max=10, log = T)
tprior = dnorm(t, sd = 5, log = T)
p1prior = dnorm(p1, sd = 5, log = T)
p2prior = dnorm(p1, sd = 5, log = T)
b1prior = dnorm(p1, sd = 5, log = T)
b2prior = dnorm(p1, sd = 5, log = T)
return(s1prior + s2prior + tprior + p1prior + p2prior + b1prior + b2prior)
}
posterior <- function(param){
return (likelihood(param) + prior(param)) #being bayesian
#return(likelihood(param))#what if i don't be bayesian?
}
proposalfunction <- function(param){
return(rnorm(7,mean = param, sd= c(0.1,0.5,0.3,.5,.5,.5,.5)))
}
run_metropolis_MCMC <- function(startvalue, iterations){
chain = array(dim = c(iterations+1,length(startvalue)))
chain[1,] = startvalue
for (i in 1:iterations){
proposal = proposalfunction(chain[i,])
probab = exp(posterior(proposal) - posterior(chain[i,]))
if (!is.na(probab) & (runif(1) < probab)){
chain[i+1,] = proposal
}else{
chain[i+1,] = chain[i,]
}
}
return(chain)
}
names = c('sig2y1','sig2x','tau','pi1','pi2','b1','b2')
true = c(sig2y1,sig2x,t0,Pi1,Pi2,B1,B2)
startvalue = c(rep(1,length(true)))
chain1 = run_metropolis_MCMC(startvalue, 1000)
chain2 = run_metropolis_MCMC(startvalue, 1000)
burnIn = 5000
acceptance = 1-mean(duplicated(chain1[-(1:burnIn),]))
for(i in 1:dim(chain1)[2]){
hist(chain1[,i],nclass=30, main=paste("Posterior of ",names[i]), xlab="True value = red line" )
abline(v = mean(chain1[,i]));abline(v = true[i], col="red" )
}
#### with loglik_lgnorm ####
rm(list = ls())
detach(dat)
ps = c(0,-1,3); gs = c(0,.8,.07); bs = c(0,.06,.02); sigv = 5; sigu = 2;r=.2;t0=.3;t1=.1
true = c(sigv,sigu,t0,t1,bs,gs,ps,r)
dat = hurdle.IV.sim(pi=ps,gamma = gs,beta = bs,endog_sd = sigv,y_sd = sigu,rho =r,tau0 = t0,tau1=t1)
dat$y1 <-dat$y
dat$x2<-dat$endog;dat$x1<-dat$exog1;dat$z<-dat$inst1
attach(dat)
likelihood<-function(param) {-loglik_lgnorm_rho(param)}
prior <- function(params){
#need to get this to be wishart
invsds = dnorm(params[1:2],sd = 100,log = T)
sds = 1/invsds
lasts = Reduce("+",lapply(params[3:length(params)],function(x) dnorm(x,100,log=T)))
lasts + sum(sds)
}
posterior <- function(param){
return (likelihood(param) + prior(param)) #being bayesian
#return(likelihood(param))#what if i don't be bayesian?
}
proposalfunction <- function(param){
return(rnorm(length(param),mean = param, sd= rep(1,length(param))))
}
run_metropolis_MCMC <- function(startvalue, iterations){
chain = array(dim = c(iterations+1,length(startvalue)))
chain[1,] = startvalue
for (i in 1:iterations){
proposal = proposalfunction(chain[i,])
probab = exp(posterior(proposal) - posterior(chain[i,]))
if ((!is.na(probab)) & runif(1) < probab){
chain[i+1,] = proposal
}else{
chain[i+1,] = chain[i,]
}
}
return(chain)
}
chain1 = run_metropolis_MCMC(true, 1e8) #bayesian from true
#chain1 = run_metropolis_MCMC(c(rep(1,14)),10000)
burnIn = 5000
acceptance = 1-mean(duplicated(chain1[-(1:burnIn),]))
true = c(sigv,sigu,t0,t1,bs,gs,ps,r)
nm_pars = c('sigv','sigu','tau0','tau1','beta1','beta2','beta3','gamma1','gamma2','gamma3','pi1','pi2','pi3','rho')
for(i in 1:length(true)){
hist(chain1[-(1:burnIn),i],nclass=30, main=paste("Posterior of",nm_pars[i])
, xlab="True value = red line" )
abline(v = mean(chain1[-(1:burnIn),i]));abline(v = true[i], col="red" )
}
jackiemauro/hurdleIV documentation built on May 18, 2019, 7:56 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
# another attempt at mcmc
#### laplace demon ####
#http://www.sumsar.net/blog/2013/06/three-ways-to-run-bayesian-models-in-r/
rm(list = ls())
library(LaplacesDemon)
y = rnorm(1000,4,2)
# The model specification
model <- function(parm, data) {
mu <- parm[1]
sigma <- exp(parm[2])
log_lik <- sum(dnorm(data$y, mu, sigma, log = T))
post_lik <- log_lik + dnorm(mu, 0, 100, log = T) + dlnorm(sigma, 0, 4, log = T)
# This list is returned and has to follow a the format specified by
# LaplacesDemon.
list(LP = post_lik, Dev = -2 * log_lik, Monitor = c(post_lik, sigma), yhat = NA,
parm = parm)
}
# Running the model
data_list <- list(N = length(y), y = y, mon.names = c("post_lik", "sigma"),
parm.names = c("mu", "log.sigma"))
mcmc_samples <- LaplacesDemon(Model = model, Data = data_list, Iterations = 30000,
Algorithm = "HARM", Thinning = 1)
plot(mcmc_samples, BurnIn = 10000, data_list)
mcmc_samples$Summary1['mu','Mean'];mcmc_samples$Summary1['mu','SD']
mcmc_samples$Summary1['sigma','Mean'];mcmc_samples$Summary1['sigma','SD']
# basic IV
rm(list = ls())
require(mvtnorm); require(MASS)
N = 10000
z <- rnorm(N,2,1)
sig2y1<-2; sig2x<-3; t0<-.2
Pi1<-1; Pi2<-3; B1<-2; B2<-1
errs = mvrnorm(n = N, mu = c(0,0), Sigma = matrix(c(sig2y1,t0,t0,sig2x),ncol = 2, byrow = T))
x2 <- Pi1 + Pi2*z + errs[,2]
y1 <- B1 + B2*x2 + errs[,1]
model <- function(parm, data) {
b1 <- parm[1]; b2 <- parm[2]; p1 <- parm[3]; p2 <- parm[4]
sigy1 <- exp(parm[5]); sigx <- exp(parm[6]); t0 <- parm[7]
mu_x <- p1 + p2*data$z
mu_y <- b1 + b2*mu_x
Sig = matrix(c(sigy1,t0,t0,sigx),ncol = 2, byrow = T)
#Transformation matrix from (eta, u, v) to (y0*, log y1*, x2) (without the mean)
A = rbind(c(1,b2),c(0,1))
#Covariance of (y0*, log y1*, x2)
Sig2 = A%*%Sig%*%t(A)
mu_y_x = mu_y + (Sig2[1,2]/Sig2[2,2])*(data$x-mu_x)
sig2_y_x = Sig2[1,1] - (Sig2[1,2]^2)/Sig2[2,2]
if(sig2_y_x>0){
log_lik <- sum(dnorm(data$y,mu_y_x,sqrt(sig2_y_x),log = T) + dnorm(data$x,mu_x,sqrt(Sig2[2,2]),log =T))
}
else{
mat = data.frame(y=data$y,x=data$x,predy=mu_y,predx = mu_x)
singlelikelihoods = apply(mat,1,function(x) dmvnorm(x[1:2],mean = x[3:4],sigma = Sig,log = T))
log_lik <- sum(singlelikelihoods)
}
log_prior <- sum(dnorm(c(b1,b2,p1,p2,t0),0,100,log = T))+dlnorm(sigy1, 0, 4, log = T)+dlnorm(sigx, 0, 4, log = T)
post_lik <- log_lik + log_prior
# This list is returned and has to follow a the format specified by
# LaplacesDemon.
list(LP = post_lik, Dev = -2 * log_lik, Monitor = c(post_lik, sigx,sigy1), yhat = NA,parm = parm)
}
# Running the model
data_list <- list(N = length(y1), y = y1,x = x2, z = z , mon.names = c("post_lik", "sigx","sigy1"),
parm.names = c("B1","B2","Pi1","Pi2","log.sigy","log.sigx","tau0"))
mcmc_samples <- LaplacesDemon(Model = model, Data = data_list, Iterations = 50000,
Algorithm = "HARM", Thinning = 1)
plot(mcmc_samples, BurnIn = 10000, data_list)
parameters = c('B1','B2','Pi1','Pi2','tau0','sigx','sigy1')
true = c(B1,B2,Pi1,Pi2,t0,sig2x,sig2y1)
names(true)<-parameters
for(par in parameters){
print(true[par] )
print(mcmc_samples$Summary1[par,'Mean'])
print(mcmc_samples$Summary1[par,'SD'])
}
# lognormal IV
rm(list = ls())
N = 10000
z <- rnorm(N,2,1); x1 <- rnorm(N,-1,3)
sig2y1<-2; sig2y0<-1; rho <- .1; t0<-.2; t1 <- .3; sig2x<-3
Pi1<-1; Pi2<-3; B1<-2; B2<-1;G1<- -1; G2 <- 2
errs = mvrnorm(n = N, mu = c(0,0,0), Sigma = matrix(c(sig2y0,rho,t0,rho,sig2y1,t1,t0,t1,sig2x),ncol = 3, byrow = T))
x2 <- Pi1*x1 + Pi2*z + errs[,3]
y0Star <- G1*x1 + G2*z + errs[,2]
lny1Star <- B1*x1 + B2*x2 + errs[,1]
y1 <- as.numeric(y0Star>0)*exp(lny1Star)
model <- function(parm, data) {
sigx <- exp(parm[1]); sigy1 <- exp(parm[2]); t1 <- parm[3]; t0 <- parm[4]
b1 <- parm[5]; b2 <- parm[6]; g1 <- parm[7]; g2 <- parm[8]; p1<-parm[9];p2<-parm[10];rho <- parm[11]
log_lik <- loglik_lgnorm_rho(parm)
log_prior <- sum(dnorm(c(b1,b2,g1,g2,p1,p2,t1,t0,rho),0,100,log = T))+dlnorm(sigy1, 0, 4, log = T)+dlnorm(sigx, 0, 4, log = T)
post_lik <- log_lik + log_prior
print(post_lik)
# This list is returned and has to follow a the format specified by
# LaplacesDemon.
list(LP = post_lik, Dev = -2 * log_lik, Monitor = c(post_lik, sigx,sigy1), yhat = NA,parm = parm)
}
# Running the model
data_list <- list(N = length(y1), y = y1,x = x2, x1=x1,z = z , mon.names = c("post_lik", "sigx","sigy1"),
parm.names = c("log.sigx","log.sigy","tau1","tau0","B1","B2","G1","G2","Pi1","Pi2","rho"))
mcmc_samples <- LaplacesDemon(Model = model, Data = data_list, Iterations = 100000,
Algorithm = "HARM", Thinning = 1,
Initial.Values = c(1,1,0,0,B1,B2,G1,G2,Pi1,Pi2,rho))
plot(mcmc_samples, BurnIn = 10000, data_list)
True.vals <- c(sig2x,sig2y1,t1,t0,B1,B2,G1,G2,Pi1,Pi2,rho)
names <- c('sigx','sigy1','tau1','tau0','B1','B2','G1','G2','Pi1','Pi2','rho')
Compare <- cbind(mcmc_samples$Summary2[names,c("Mean","SD")],True.vals)
# run it within the hurdleIVMCMC code
rm(list = ls())
z <- rnorm(N,2,1); x1 <- rnorm(N,-1,3)
sig2y1<-2; sig2y0<-1; rho <- .1; t0<-.2; t1 <- .3; sig2x<-3
Pi1<-1; Pi2<-3; B1<-2; B2<-1;G1<- -1; G2 <- 2
bs = c(0,B1,B2); gs = c(0,G1,G2); pis = c(0, Pi1, Pi2)
sig2y1<-2; sig2y0<-1; rho <- .1; t0<-.2; t1 <- .3; sig2x<-3
inst_mean = 2; inst_sd = 1; exog_mean = -1; exog_sd = 3
nms = c('rho','sigy1','tau0','tau1','sigx','beta1','beta2','beta3','gamma1','gamma2','gamma3','pi1','pi2','pi3')
true.vals = c(rho,sig2y1,t0,t1,sig2x,bs,gs,pis)
dat = hurdle.IV.sim(pi = pis, gamma = gs, beta = bs,
exog_mean = exog_mean, exog_sd = exog_sd,
z_mean = inst_mean, z_sd = inst_sd,
endog_sd = sqrt(sig2x), y_sd = sqrt(sig2y1),
rho = rho, tau0 = t0, tau1 = t1, n = 10000)
out = hurdle.IV.MCMC(y~exog1 + endog,
inst = inst1,
endog = endog,
exog = exog1,
data = dat,
k = 50000)
Compare <- cbind(out$Summary2[nms,c("Mean","SD")],true.vals)
# ##### this is from: https://nicercode.github.io/guides/mcmc/ ####
# just to get mean
m = 0
f <- function(x) dnorm(x,m,1)
q <- function(x) rnorm(1,x,1)
step <- function(x,f,q){
xp = q(x) #propose a new x
alpha = min(1,f(xp)/f(x)) #calculate alpha
if(runif(1)<alpha){x<-xp}
x
}
run <- function(x, f, q, nsteps) {
res <- matrix(NA, nsteps, length(x))
for (i in seq_len(nsteps))
res[i,] <- x <- step(x, f, q)
drop(res)
}
res <- run(-10, f, q, 1000)
layout(matrix(c(1, 2), 1, 2), widths=c(4, 1))
par(mar=c(4.1, .5, .5, .5), oma=c(0, 4.1, 0, 0))
plot(res, type="s", xpd=NA, ylab="Parameter", xlab="Sample", las=1)
usr <- par("usr")
xx <- seq(usr[3], usr[4], length=301)
plot(f(xx), xx, type="l", yaxs="i", axes=FALSE, xlab="")
#### looking at these guys' example: https://theoreticalecology.wordpress.com/2010/09/17/metropolis-hastings-mcmc-in-r/####
rm(list = ls())
trueA <- 5
trueB <- 0
trueSd <- 10
sampleSize <- 31
# create independent x-values
x <- (-(sampleSize-1)/2):((sampleSize-1)/2)
# create dependent values according to ax + b + N(0,sd)
y <- trueA * x + trueB + rnorm(n=sampleSize,mean=0,sd=trueSd)
likelihood <- function(param){
a = param[1]
b = param[2]
sd = param[3]
pred = a*x + b
singlelikelihoods = dnorm(y, mean = pred, sd = sd, log = T)
sumll = sum(singlelikelihoods)
return(sumll)
}
prior <- function(param){
a = param[1]
b = param[2]
sd = param[3]
aprior = dunif(a, min=0, max=10, log = T)
bprior = dnorm(b, sd = 5, log = T)
sdprior = dunif(sd, min=0, max=30, log = T)
return(aprior+bprior+sdprior)
}
posterior <- function(param){
return (likelihood(param) + prior(param)) #being bayesian
#return(likelihood(param))#what if i don't be bayesian?
}
proposalfunction <- function(param){
return(rnorm(3,mean = param, sd= c(0.1,0.5,0.3)))
}
run_metropolis_MCMC <- function(startvalue, iterations){
chain = array(dim = c(iterations+1,3))
chain[1,] = startvalue
for (i in 1:iterations){
proposal = proposalfunction(chain[i,])
probab = exp(posterior(proposal) - posterior(chain[i,]))
if (runif(1) < probab){
chain[i+1,] = proposal
}else{
chain[i+1,] = chain[i,]
}
}
return(chain)
}
startvalue = c(4,0,10)
chain1 = run_metropolis_MCMC(startvalue, 10000) #bayesian
chain2 = run_metropolis_MCMC(startvalue, 10000) #not bayesian (change posterior)
burnIn = 5000
acceptance = 1-mean(duplicated(chain1[-(1:burnIn),]))
names = c('a','b','sig'); true = c(trueA,trueB,trueSd)
for(i in 1:dim(chain1)[2]){
hist(chain1[-(1:burnIn),i],nclass=30, main=paste("Posterior of ",names[i]), xlab="True value = red line" )
abline(v = mean(chain1[-(1:burnIn),i]));abline(v = true[i], col="red" )
}
#### basic IV #####
# even though dmvnorm and dnorm using conditionals give the same probabilities
# asking the algorithm to estimate the mean of x/y within the likelihood makes it
# perform less well
rm(list = ls())
require(mvtnorm)
N = 10000
z <- rnorm(N,2,1)
sig2y1<-2; sig2x<-3; t0<-.2
Pi1<-1; Pi2<-3; B1<-2; B2<-1
errs = mvrnorm(n = N, mu = c(0,0), Sigma = matrix(c(sig2y1,t0,t0,sig2x),ncol = 2, byrow = T))
x <- Pi1 + Pi2*z + errs[,2]
y <- B1 + B2*x + errs[,1]
likelihood <- function(param){
s1 = param[1];s2 = param[2]; t = param[3]; p1 = param[4]; p2 = param[5]; b1 = param[6]; b2 = param[7]
xpred = p1 + p2*z
ypred = b1 + b2*xpred
preds = matrix(c(ypred,xpred),ncol = 2, byrow = F)
obs = matrix(c(y,x),ncol = 2, byrow = F)
Sig = matrix(c(s1,t,t,s2),ncol = 2, byrow = T)
# These two will give the same likelihoods, except, I think, when sig2_y1_x2<0
# Sig_err = matrix(c(s1,t,t,s2),ncol = 2, byrow = T)
# A = matrix(c(1,b1,0,1),ncol = 2, byrow = T)
# Sig = A%*%Sig_err%*%t(A)
#
# mu_y1_x2 = ypred + Sig[1,2]/Sig[2,2]*(x-xpred)
# sig2_y1_x2 = Sig[1,1] - Sig[1,2]^2/Sig[2,2]
#
# ylik = sum(dnorm(y,mean = mu_y1_x2,sd = sqrt(sig2_y1_x2),log = T))
# xlik = sum(dnorm(x,mean = xpred,sd = sqrt(Sig[2,2]),log = T))
# sumll1 = xlik + ylik
#
mat = data.frame(y=y,x=x,predy=ypred,predx = xpred)
singlelikelihoods = apply(mat,1,function(x) dmvnorm(x[1:2],mean = x[3:4],sigma = Sig,log = T))
sumll = sum(singlelikelihoods)
#
# print(sumll)
# print(sumll1)
return(sumll)
}
likelihood <- function(param){
# this one has to guess at the true value of x in order to get the conditional
# mean of y--faster but appears not to be as accurate.
s1 = param[1];s2 = param[2]; t = param[3]; p1 = param[4]; p2 = param[5]; b1 = param[6]; b2 = param[7]
xpred = p1 + p2*z
ypred = b1 + b2*xpred
preds = matrix(c(xpred,ypred),ncol = 2, byrow = F)
obs = matrix(c(x,y),ncol = 2, byrow = F)
#Sig = matrix(c(s1,t,t,s2),ncol = 2, byrow = T)
Sig_err = matrix(c(s1,t,t,s2),ncol = 2, byrow = T)
A = matrix(c(1,b1,0,1),ncol = 2, byrow = T)
Sig = A%*%Sig_err%*%t(A)
mu_y1_x2 = ypred + Sig[1,2]/Sig[2,2]*(x-xpred)
sig2_y1_x2 = Sig[1,1] - Sig[1,2]^2/Sig[2,2]
if(sig2_y1_x2<=0 | Sig[2,2]<=0){return(-Inf)}
ylik = sum(dnorm(y,mean = mu_y1_x2,sd = sqrt(sig2_y1_x2),log = T))
xlik = sum(dnorm(x,mean = xpred,sd = sqrt(Sig[2,2]),log = T))
sumll = xlik + ylik
return(sumll)
}
prior <- function(param){
s1 = param[1];s2 = param[2]; t = param[3]; p1 = param[4]; p2 = param[5]; b1 = param[6]; b2 = param[7]
s1prior = dunif(s1, min=0, max=10, log = T)
s2prior = dunif(s2, min=0, max=10, log = T)
tprior = dnorm(t, sd = 5, log = T)
p1prior = dnorm(p1, sd = 5, log = T)
p2prior = dnorm(p1, sd = 5, log = T)
b1prior = dnorm(p1, sd = 5, log = T)
b2prior = dnorm(p1, sd = 5, log = T)
return(s1prior + s2prior + tprior + p1prior + p2prior + b1prior + b2prior)
}
posterior <- function(param){
return (likelihood(param) + prior(param)) #being bayesian
#return(likelihood(param))#what if i don't be bayesian?
}
proposalfunction <- function(param){
return(rnorm(7,mean = param, sd= c(0.1,0.5,0.3,.5,.5,.5,.5)))
}
run_metropolis_MCMC <- function(startvalue, iterations){
chain = array(dim = c(iterations+1,length(startvalue)))
chain[1,] = startvalue
for (i in 1:iterations){
proposal = proposalfunction(chain[i,])
probab = exp(posterior(proposal) - posterior(chain[i,]))
if (!is.na(probab) & (runif(1) < probab)){
chain[i+1,] = proposal
}else{
chain[i+1,] = chain[i,]
}
}
return(chain)
}
names = c('sig2y1','sig2x','tau','pi1','pi2','b1','b2')
true = c(sig2y1,sig2x,t0,Pi1,Pi2,B1,B2)
startvalue = c(rep(1,length(true)))
chain1 = run_metropolis_MCMC(startvalue, 1000)
chain2 = run_metropolis_MCMC(startvalue, 1000)
burnIn = 5000
acceptance = 1-mean(duplicated(chain1[-(1:burnIn),]))
for(i in 1:dim(chain1)[2]){
hist(chain1[,i],nclass=30, main=paste("Posterior of ",names[i]), xlab="True value = red line" )
abline(v = mean(chain1[,i]));abline(v = true[i], col="red" )
}
xmins1 = apply(chain1,2,min);xmins2 = apply(chain2,2,min);xmins = apply(cbind(xmins1,xmins2),1,min)
xmaxs1 = apply(chain1,2,max);xmaxs2 = apply(chain2,2,max);xmaxs = apply(cbind(xmaxs1,xmaxs2),1,max)
par(mfrow = c(1,2))
for(i in 1:dim(chain2)[2]){
hist(chain1[,i],nclass=30, main=paste("Posterior of ",names[i]," chain 1"), xlim = c(xmins[i],xmaxs[i]),xlab="True value = red line" )
abline(v = mean(chain1[,i]));abline(v = true[i], col="red" )
hist(chain2[,i],nclass=30, main=paste("Posterior of ",names[i], "chain 2"), xlim = c(xmins[i],xmaxs[i]), xlab="True value = red line" )
abline(v = mean(chain2[,i]));abline(v = true[i], col="red" )
}
par(mfrow = c(1,1))
#### with fake loglik lgnorm ####
rm(list = ls())
require(mvtnorm)
N = 10000
z <- rnorm(N,2,1)
sig2y1<-2; sig2y0<-1; rho <- .1; t0<-.2; t1 <- .3; sig2x<-3
Pi1<-1; Pi2<-3; B1<-2; B2<-1;G1<- -1; G2 <- 2
errs = mvrnorm(n = N, mu = c(0,0,0), Sigma = matrix(c(sig2y0,rho,t0,rho,sig2y1,t1,t0,t1,sig2x),ncol = 3, byrow = T))
x <- Pi1 + Pi2*z + errs[,3]
y0Star <- G1 + G2*z + errs[,2]
y1Star <- B1 + B2*x + errs[,1]
y1 <- as.numeric(y0Star>0)*y1Star
likelihood <- function(param){ #not done
sy0 = param[1]; sy1 = param[2]; sx2 = param[3];
r = param[4]; t0 = param[5]; t1 = param[6]
p1 = param[7]; p2 = param[8];g1 = param[9];g2 = param[10]; b1 = param[11]; b2 = param[12]
xpred = p1 + p2*z
y0pred = g1 + g2*xpred
y1pred = b1 + b2*xpred
preds = matrix(c(ypred,xpred),ncol = 2, byrow = F)
obs = matrix(c(y,x),ncol = 2, byrow = F)
Sig = matrix(c(sy0,r,t0,r,sy1,t1,t0,t1,sx2),ncol = 3, byrow = T)
mat = data.frame(y=y,x=x,predy=ypred,predx = xpred)
singlelikelihoods = apply(mat,1,function(x) dmvnorm(x[1:2],mean = x[3:4],sigma = Sig,log = T))
sumll = sum(singlelikelihoods)
return(sumll)
}
prior <- function(param){
s1 = param[1];s2 = param[2]; t = param[3]; p1 = param[4]; p2 = param[5]; b1 = param[6]; b2 = param[7]
s1prior = dunif(s1, min=0, max=10, log = T)
s2prior = dunif(s2, min=0, max=10, log = T)
tprior = dnorm(t, sd = 5, log = T)
p1prior = dnorm(p1, sd = 5, log = T)
p2prior = dnorm(p1, sd = 5, log = T)
b1prior = dnorm(p1, sd = 5, log = T)
b2prior = dnorm(p1, sd = 5, log = T)
return(s1prior + s2prior + tprior + p1prior + p2prior + b1prior + b2prior)
}
posterior <- function(param){
return (likelihood(param) + prior(param)) #being bayesian
#return(likelihood(param))#what if i don't be bayesian?
}
proposalfunction <- function(param){
return(rnorm(7,mean = param, sd= c(0.1,0.5,0.3,.5,.5,.5,.5)))
}
run_metropolis_MCMC <- function(startvalue, iterations){
chain = array(dim = c(iterations+1,length(startvalue)))
chain[1,] = startvalue
for (i in 1:iterations){
proposal = proposalfunction(chain[i,])
probab = exp(posterior(proposal) - posterior(chain[i,]))
if (!is.na(probab) & (runif(1) < probab)){
chain[i+1,] = proposal
}else{
chain[i+1,] = chain[i,]
}
}
return(chain)
}
names = c('sig2y1','sig2x','tau','pi1','pi2','b1','b2')
true = c(sig2y1,sig2x,t0,Pi1,Pi2,B1,B2)
startvalue = c(rep(1,length(true)))
chain1 = run_metropolis_MCMC(startvalue, 1000)
chain2 = run_metropolis_MCMC(startvalue, 1000)
burnIn = 5000
acceptance = 1-mean(duplicated(chain1[-(1:burnIn),]))
for(i in 1:dim(chain1)[2]){
hist(chain1[,i],nclass=30, main=paste("Posterior of ",names[i]), xlab="True value = red line" )
abline(v = mean(chain1[,i]));abline(v = true[i], col="red" )
}
#### with loglik_lgnorm ####
rm(list = ls())
detach(dat)
ps = c(0,-1,3); gs = c(0,.8,.07); bs = c(0,.06,.02); sigv = 5; sigu = 2;r=.2;t0=.3;t1=.1
true = c(sigv,sigu,t0,t1,bs,gs,ps,r)
dat = hurdle.IV.sim(pi=ps,gamma = gs,beta = bs,endog_sd = sigv,y_sd = sigu,rho =r,tau0 = t0,tau1=t1)
dat$y1 <-dat$y
dat$x2<-dat$endog;dat$x1<-dat$exog1;dat$z<-dat$inst1
attach(dat)
likelihood<-function(param) {-loglik_lgnorm_rho(param)}
prior <- function(params){
#need to get this to be wishart
invsds = dnorm(params[1:2],sd = 100,log = T)
sds = 1/invsds
lasts = Reduce("+",lapply(params[3:length(params)],function(x) dnorm(x,100,log=T)))
lasts + sum(sds)
}
posterior <- function(param){
return (likelihood(param) + prior(param)) #being bayesian
#return(likelihood(param))#what if i don't be bayesian?
}
proposalfunction <- function(param){
return(rnorm(length(param),mean = param, sd= rep(1,length(param))))
}
run_metropolis_MCMC <- function(startvalue, iterations){
chain = array(dim = c(iterations+1,length(startvalue)))
chain[1,] = startvalue
for (i in 1:iterations){
proposal = proposalfunction(chain[i,])
probab = exp(posterior(proposal) - posterior(chain[i,]))
if ((!is.na(probab)) & runif(1) < probab){
chain[i+1,] = proposal
}else{
chain[i+1,] = chain[i,]
}
}
return(chain)
}
chain1 = run_metropolis_MCMC(true, 1e8) #bayesian from true
#chain1 = run_metropolis_MCMC(c(rep(1,14)),10000)
burnIn = 5000
acceptance = 1-mean(duplicated(chain1[-(1:burnIn),]))
true = c(sigv,sigu,t0,t1,bs,gs,ps,r)
nm_pars = c('sigv','sigu','tau0','tau1','beta1','beta2','beta3','gamma1','gamma2','gamma3','pi1','pi2','pi3','rho')
for(i in 1:length(true)){
hist(chain1[-(1:burnIn),i],nclass=30, main=paste("Posterior of",nm_pars[i])
, xlab="True value = red line" )
abline(v = mean(chain1[-(1:burnIn),i]));abline(v = true[i], col="red" )
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.