Nothing
## ---- echo = FALSE------------------------------------------------------------
library(bayesm)
knitr::opts_chunk$set(fig.align = "center",
fig.height = 3.5,
warning = FALSE,
error = FALSE,
message = FALSE)
## -----------------------------------------------------------------------------
# define function
drawprior <- function (mubar_betak, nvar, ncomp, a, nu, Amu, V, ndraw) {
betakstar <- double(ndraw)
betak <- double(ndraw)
otherbeta <- double(ndraw)
mubar <- c(rep(0, nvar-1), mubar_betak)
for(i in 1:ndraw) {
comps=list()
for(k in 1:ncomp) {
Sigma <- rwishart(nu,chol2inv(chol(V)))$IW
comps[[k]] <- list(mubar + t(chol(Sigma/Amu)) %*% rnorm(nvar),
backsolve(chol(Sigma), diag(1,nvar)) )
}
pvec <- rdirichlet(a)
beta <- rmixture(1,pvec,comps)$x
betakstar[i] <- beta[nvar]
betak[i] <- -exp(beta[nvar])
otherbeta[i] <- beta[1]
}
return(list(betakstar=betakstar, betak=betak, otherbeta=otherbeta))
}
set.seed(1234)
## -----------------------------------------------------------------------------
# specify rhierMnlRwMixture defaults
mubar_betak <- 0
nvar <- 10
ncomp <- 3
a <- rep(5, ncomp)
nu <- nvar + 3
Amu <- 0.01
V <- nu*diag(c(rep(1,nvar-1),1))
ndraw <- 10000
defaultprior <- drawprior(mubar_betak, nvar, ncomp, a, nu, Amu, V, ndraw)
## ---- fig.align='center', fig.height=3.5, results='hold'----------------------
# plot priors under defaults
par(mfrow=c(1,3))
trimhist <- -20
hist(defaultprior$betakstar, breaks=40, col="magenta",
main="Beta_k_star", xlab="", ylab="", yaxt="n")
hist(defaultprior$betak[defaultprior$betak>trimhist],
breaks=40, col="magenta", main="Beta_k",
xlab="", ylab="", yaxt="n", xlim=c(trimhist,0))
hist(defaultprior$otherbeta, breaks=40, col="magenta",
main="Other Beta", xlab="", ylab="", yaxt="n")
## -----------------------------------------------------------------------------
# adjust priors for constraints
mubar_betak <- 2
nvar <- 10
ncomp <- 3
a <- rep(5, ncomp)
nu <- nvar + 15
Amu <- 0.1
V <- nu*diag(c(rep(4,nvar-1),0.1))
ndraw <- 10000
tightprior <- drawprior(mubar_betak, nvar, ncomp, a, nu, Amu, V, ndraw)
## ---- fig.align='center', fig.height=3.5, results='hold'----------------------
# plot priors under adjusted values
par(mfrow=c(1,3))
trimhist <- -20
hist(tightprior$betakstar, breaks=40, col="magenta",
main="Beta_k_star", xlab="", ylab="", yaxt="n")
hist(tightprior$betak[tightprior$betak>trimhist],
breaks=40, col="magenta", main="Beta_k",
xlab="", ylab="", yaxt="n", xlim=c(trimhist,0))
hist(tightprior$otherbeta, breaks=40, col="magenta",
main="Other Beta", xlab="", ylab="", yaxt="n")
## -----------------------------------------------------------------------------
library(bayesm)
data(camera)
length(camera)
str(camera[[1]])
## -----------------------------------------------------------------------------
str(camera[[1]]$y)
str(as.data.frame(camera[[1]]$X))
## ---- echo = FALSE------------------------------------------------------------
nvar <- 10
mubar <- c(rep(0,nvar-1),2)
Amu <- 0.1
ncomp <- 5
a <- rep(5, ncomp)
nu <- 25
V <- nu*diag(c(rep(4,nvar-1),0.1))
## ---- eval = FALSE------------------------------------------------------------
# SignRes <- c(rep(0,nvar-1),-1)
#
# data <- list(lgtdata=camera, p=5)
# prior <- list(mubar=mubar, Amu=Amu, ncomp=ncomp, a=a, nu=nu, V=V, SignRes=SignRes)
# mcmc <- list(R=1e4, nprint=0)
#
# out <- rhierMnlRwMixture(Data=data, Prior=prior, Mcmc=mcmc)
## ---- echo = FALSE------------------------------------------------------------
temp <- capture.output(
{SignRes <- c(rep(0,nvar-1),-1);
data <- list(lgtdata = camera, p = 5);
prior <- list(mubar=mubar, Amu=Amu, ncomp=ncomp, a=a, nu=nu, V=V, SignRes=SignRes);
mcmc <- list(R = 1e4, nprint = 0);
out <- rhierMnlRwMixture(Data = data, Prior = prior, Mcmc = mcmc)},
file = NULL)
## -----------------------------------------------------------------------------
par(mfrow=c(1,3))
ind_hist <- function(mod, i) {
hist(mod$betadraw[i , 10, ], breaks = seq(-14,0,0.5),
col = "dodgerblue3", border = "grey", yaxt = "n",
xlim = c(-14,0), xlab = "", ylab = "", main = paste("Ind.",i))
}
ind_hist(out,1)
ind_hist(out,2)
ind_hist(out,3)
## -----------------------------------------------------------------------------
par(mfrow=c(1,1))
hist(apply(out$betadraw[ , 10, ], 1, mean),
xlim = c(-20,0), breaks = 20,
col = "firebrick2", border = "gray", xlab = "", ylab = "",
main = "Posterior Means for Ind. Price Params,
With Sign Constraint")
## -----------------------------------------------------------------------------
data0 <- list(lgtdata = camera, p = 5)
prior0 <- list(ncomp = 5)
mcmc0 <- list(R = 1e4, nprint = 0)
## ---- eval = FALSE------------------------------------------------------------
# out0 <- rhierMnlRwMixture(Data = data0, Prior = prior0, Mcmc = mcmc0)
## ---- echo = FALSE------------------------------------------------------------
temp <- capture.output(
{out0 <- rhierMnlRwMixture(Data = data0, Prior = prior0, Mcmc = mcmc0)},
file = NULL)
## -----------------------------------------------------------------------------
par(mfrow=c(1,3))
ind_hist <- function(mod, i) {
hist(mod$betadraw[i , 10, ], breaks = seq(-12,2,0.5),
col = "dodgerblue4", border = "grey", yaxt = "n",
xlim = c(-12,2), xlab = "", ylab = "", main = paste("Ind.",i))
}
ind_hist(out0,1)
ind_hist(out0,2)
ind_hist(out0,3)
## -----------------------------------------------------------------------------
par(mfrow=c(1,1))
hist(apply(out0$betadraw[ , 10, ], 1, mean),
xlim = c(-15,5), breaks = 20,
col = "firebrick3", border = "gray",
xlab = "", ylab = "",
main = "Posterior Means for Ind. Price Params,
No Sign Constraint")
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.