# inst/doc/Constrained_MNL_Vignette.R In bayesm: Bayesian Inference for Marketing/Micro-Econometrics

```## ---- 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")

## ------------------------------------------------------------------------
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")
```

## Try the bayesm package in your browser

Any scripts or data that you put into this service are public.

bayesm documentation built on Dec. 21, 2018, 9:04 a.m.