[r-cran-bayesm] 27/44: Import Upstream version 2.1-2
Andreas Tille
tille at debian.org
Thu Sep 7 11:16:22 UTC 2017
This is an automated email from the git hooks/post-receive script.
tille pushed a commit to branch master
in repository r-cran-bayesm.
commit 2f48d42d57895c428be3c7cfd216d9d7df8ff0b2
Author: Andreas Tille <tille at debian.org>
Date: Thu Sep 7 13:09:35 2017 +0200
Import Upstream version 2.1-2
---
DESCRIPTION | 13 +-
NAMESPACE | 21 +++-
R/init.rmultiregfp.R | 28 -----
R/mixDen.R | 3 +-
R/plot.bayesm.hcoef.R | 43 +++++++
R/plot.bayesm.mat.R | 57 +++++++++
R/plot.bayesm.nmix.R | 80 +++++++++++++
R/rbiNormGibbs.R | 2 +
R/rbprobitGibbs.R | 10 +-
R/rhierBinLogit.R | 18 +--
R/rhierLinearMixture.R | 14 ++-
R/rhierLinearModel.R | 13 +-
R/rhierMnlRwMixture.R | 11 +-
R/rhierNegbinRw.R | 17 ++-
R/rivGibbs.R | 14 +++
R/rmnlIndepMetrop.R | 6 +-
R/rmnpGibbs.R | 11 +-
R/rmultiregfp.R | 53 ---------
R/rmvpGibbs.R | 8 ++
R/rnegbinRw.R | 17 ++-
R/rnmixGibbs.R | 6 +-
R/rordprobitGibbs.R | 290 +++++++++++++++++++++++++++++++++++++++++++++
R/rscaleUsage.R | 21 +++-
R/rsurGibbs.R | 8 +-
R/runireg.R | 7 ++
R/runiregGibbs.R | 5 +
R/simmnl.R | 51 --------
R/simmnlwX.R | 41 -------
R/simmnp.R | 27 -----
R/simmvp.R | 18 ---
R/summary.bayesm.mat.R | 48 ++++++++
R/summary.bayesm.nmix.R | 35 ++++++
R/summary.bayesm.var.R | 37 ++++++
data/orangeJuice.rda | Bin 0 -> 1156010 bytes
data/tuna.rda | Bin 0 -> 48101 bytes
inst/doc/bayesm-manual.pdf | Bin 441639 -> 464284 bytes
man/Scotch.Rd | 10 +-
man/bank.Rd | 32 +++--
man/cheese.Rd | 30 ++---
man/customerSat.Rd | 2 +-
man/detailing.Rd | 20 ++--
man/init.rmultiregfp.Rd | 46 -------
man/margarine.Rd | 20 ++--
man/orangeJuice.Rd | 165 ++++++++++++++++++++++++++
man/plot.bayesm.hcoef.Rd | 42 +++++++
man/plot.bayesm.mat.Rd | 51 ++++++++
man/plot.bayesm.nmix.Rd | 46 +++++++
man/rbprobitGibbs.Rd | 16 ++-
man/rhierBinLogit.Rd | 25 ++--
man/rhierLinearMixture.Rd | 120 +++----------------
man/rhierLinearModel.Rd | 31 ++---
man/rhierMnlRwMixture.Rd | 117 ++++++------------
man/rhierNegbinRw.Rd | 45 +++----
man/rivGibbs.Rd | 29 +++--
man/rmnlIndepMetrop.Rd | 38 +++++-
man/rmnpGibbs.Rd | 44 +++++--
man/rmultireg.Rd | 1 -
man/rmultiregfp.Rd | 48 --------
man/rmvpGibbs.Rd | 30 +++--
man/rnegbinRw.Rd | 32 +++--
man/rnmixGibbs.Rd | 74 +++---------
man/rordprobitGibbs.Rd | 115 ++++++++++++++++++
man/rscaleUsage.Rd | 8 +-
man/rsurGibbs.Rd | 26 ++--
man/runireg.Rd | 16 ++-
man/runiregGibbs.Rd | 22 ++--
man/simmnl.Rd | 42 -------
man/simmnlwX.Rd | 36 ------
man/simmnp.Rd | 44 -------
man/simmvp.Rd | 43 -------
man/summary.bayesm.mat.Rd | 42 +++++++
man/summary.bayesm.nmix.Rd | 38 ++++++
man/summary.bayesm.var.Rd | 43 +++++++
man/tuna.Rd | 110 +++++++++++++++++
74 files changed, 1752 insertions(+), 980 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index c831f0c..86448c6 100755
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,23 +1,24 @@
Package: bayesm
-Version: 2.0-9
-Date: 2007-01-15
+Version: 2.1-2
+Date: 2007-03-15
Title:Bayesian Inference for Marketing/Micro-econometrics
Author: Peter Rossi <peter.rossi at ChicagoGsb.edu>,
Rob McCulloch <robert.mcculloch at ChicagoGsb.edu>.
Maintainer: Peter Rossi <peter.rossi at chicagosb.edu>
-Depends: R (>= 1.8.1)
+Depends: R (>= 2.2.0)
Description: bayesm covers many important models used
in marketing and micro-econometrics applications.
The package includes:
Bayes Regression (univariate or multivariate dep var),
Bayes Seemingly Unrelated Regression (SUR),
+ Binary and Ordinal Probit,
Multinomial Logit (MNL) and Multinomial Probit (MNP),
Multivariate Probit,
Negative Binomial (Poisson) Regression,
Multivariate Mixtures of Normals (including clustering),
Hierarchical Linear Models with normal prior and covariates,
- Hierarchical Linear Models with mixture of normals prior and covariates,
- Hierarchical Multinomial Logits with mixture of normals prior
+ Hierarchical Linear Models with a mixture of normals prior and covariates,
+ Hierarchical Multinomial Logits with a mixture of normals prior
and covariates,
Hierarchical Negative Binomial Regression Models,
Bayesian analysis of choice-based conjoint data,
@@ -29,4 +30,4 @@ Description: bayesm covers many important models used
Marketing by Rossi, Allenby and McCulloch.
License: GPL (version 2 or later)
URL: http://faculty.chicagogsb.edu/peter.rossi/research/bsm.html
-Packaged: Wed Jan 17 16:27:59 2007; per
+Packaged: Thu Mar 15 17:10:45 2007; per
diff --git a/NAMESPACE b/NAMESPACE
index c713490..f640295 100755
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -2,12 +2,23 @@ useDynLib(bayesm)
export(breg,cgetC,createX,eMixMargDen,mixDen,fsh,llmnl,llmnp,llnhlogit,
lndIChisq,lndIWishart,lndMvn,lndMvst,mnlHess,momMix,nmat,numEff,rdirichlet,
-rmixture,rmultireg,rmultiregfp,rwishart,rmvst,rtrun,rbprobitGibbs,runireg,
-runiregGibbs,simmnl,simmnp,simmvp,simnhlogit,rmnpGibbs,rmixGibbs,rnmixGibbs,
-rmvpGibbs,rhierLinearModel,rhierMnlRwMixture,simmnlwX,rivGibbs,
-rmnlIndepMetrop,rscaleUsage,ghkvec,condMom,logMargDenNR,init.rmultiregfp,
+rmixture,rmultireg,rwishart,rmvst,rtrun,rbprobitGibbs,runireg,
+runiregGibbs,simnhlogit,rmnpGibbs,rmixGibbs,rnmixGibbs,
+rmvpGibbs,rhierLinearModel,rhierMnlRwMixture,rivGibbs,
+rmnlIndepMetrop,rscaleUsage,ghkvec,condMom,logMargDenNR,
rhierBinLogit,rnegbinRw,rhierNegbinRw,rbiNormGibbs,clusterMix,rsurGibbs,
-mixDenBi,mnpProb,rhierLinearMixture)
+mixDenBi,mnpProb,rhierLinearMixture,summary.bayesm.mat,plot.bayesm.mat,
+plot.bayesm.hcoef,plot.bayesm.nmix,rordprobitGibbs)
+
+
+## register S3 methods
+S3method(plot, bayesm.mat)
+S3method(plot, bayesm.nmix)
+S3method(plot, bayesm.hcoef)
+S3method(summary, bayesm.mat)
+S3method(summary, bayesm.var)
+S3method(summary, bayesm.nmix)
+
diff --git a/R/init.rmultiregfp.R b/R/init.rmultiregfp.R
deleted file mode 100755
index 2a61535..0000000
--- a/R/init.rmultiregfp.R
+++ /dev/null
@@ -1,28 +0,0 @@
-init.rmultiregfp=
-function(X,A,Bbar,nu,V)
-{
-#
-# revision history:
-# revised 1/11/05 by P. Rossi
-# purpose:
-# prepare parameters for rmultiregfp
-# rmultiregfp is designed to take advantage of the fact that posterior calculations
-# don't change X or prior parms even when the "data" or Y is changing
-#
-# use to create list for call of rmultiregfp as in:
-# Fparm=init.rmultiregfp(X,A,Bbar,nu,V)
-# rmultiregfp(Y,X,Fparm)
-#
-#
-# for multivariate regression model Y = XB + U Y is n x m, X, n x k, B is k x m
-# ith row of U is N(0,Sigma)
-#
-# prior parms: vec(B)| Sigma ~ N(vec(Bbar),Sigma (x) A^-1)
-# Sigma ~ IW(nu,V)
-#
-k=ncol(X)
-IR=backsolve(chol(crossprod(X)+A),diag(k))
-RA=chol(A)
-RABbar=RA%*%Bbar
-return(list(IR=IR,RA=RA,RABbar=RABbar,nu=nu,V=V))
-}
diff --git a/R/mixDen.R b/R/mixDen.R
index f911bdd..d1cfbcc 100755
--- a/R/mixDen.R
+++ b/R/mixDen.R
@@ -4,6 +4,7 @@ function(x,pvec,comps)
# Revision History:
# R. McCulloch 11/04
# P. Rossi 3/05 -- put in backsolve
+# P. Rossi 1/06 -- put in crossprod
#
# purpose: compute marginal densities for multivariate mixture of normals (given by p and comps) at x
#
@@ -36,7 +37,7 @@ for(i in 1:nc) {
mu[i,] = comps[[i]][[1]]
# root = solve(comps[[i]][[2]])
root= backsolve(comps[[i]][[2]],diag(rep(1,dim)))
- sigma[i,] = sqrt(diag(t(root)%*%root))
+ sigma[i,] = sqrt(diag(crossprod(root)))
}
return(list(mu=mu,sigma=sigma))
}
diff --git a/R/plot.bayesm.hcoef.R b/R/plot.bayesm.hcoef.R
new file mode 100755
index 0000000..b031a91
--- /dev/null
+++ b/R/plot.bayesm.hcoef.R
@@ -0,0 +1,43 @@
+plot.bayesm.hcoef=function(x,burnin=trunc(.1*R),...){
+#
+# S3 method to plot arrays of draws of coefs in hier models
+# 3 dimensional arrays: unit x var x draw
+# P. Rossi 2/07
+#
+ X=x
+ if(mode(X) == "list") stop("list entered \n Possible Fixup: extract from list \n")
+ if(mode(X) !="numeric") stop("Requires numeric argument \n")
+ d=dim(X)
+ if(length(d) !=3) stop("Requires 3-dim array \n")
+ op=par(no.readonly=TRUE)
+ on.exit(par(op))
+ nunits=d[1]
+ nvar=d[2]
+ R=d[3]
+ if(R < 100) {cat("fewer than 100 draws submitted \n"); return(invisible())}
+ #
+ # plot posterior distributions of nvar coef for 30 rand units
+ #
+ rsam=sort(sample(c(1:nunits),30)) # randomly sample 30 cross-sectional units
+ par(mfrow=c(1,1))
+ par(las=3) # horizontal labeling
+ for(var in 1:nvar){
+ ext=X[rsam,var,(burnin+1):R]; ext=data.frame(t(ext))
+ colnames(ext)=as.character(rsam)
+ out=boxplot(ext,plot=FALSE,...)
+ out$stats=apply(ext,2,quantile,probs=c(0,.05,.95,1))
+ bxp(out,xlab="Cross-sectional Unit",main=paste("Coefficients on Var ",var,sep=""),boxfill="magenta",...)
+ if(var==1) par(ask=dev.interactive())
+ }
+ #
+ # plot posterior means for each var
+ #
+ par(las=1)
+ pmeans=matrix(0,nrow=nunits,ncol=nvar)
+ for(i in 1:nunits) pmeans[i,]=apply(X[i,,(burnin+1):R],1,mean)
+ names=as.character(1:nvar)
+ attributes(pmeans)$class="bayesm.mat"
+ for(i in 1:nvar) names[i]=paste("Posterior Means of Coef ",i,sep="")
+ plot(pmeans,names,TRACEPLOT=FALSE,INT=FALSE,DEN=FALSE,...)
+invisible()
+}
diff --git a/R/plot.bayesm.mat.R b/R/plot.bayesm.mat.R
new file mode 100755
index 0000000..d714239
--- /dev/null
+++ b/R/plot.bayesm.mat.R
@@ -0,0 +1,57 @@
+plot.bayesm.mat=function(x,names,burnin=trunc(.1*nrow(X)),tvalues,TRACEPLOT=TRUE,DEN=TRUE,INT=TRUE,...){
+#
+# S3 method to print matrices of draws the object X is of class "bayesm.mat"
+#
+# P. Rossi 2/07
+#
+ X=x
+ if(mode(X) == "list") stop("list entered \n Possible Fixup: extract from list \n")
+ if(mode(X) !="numeric") stop("Requires numeric argument \n")
+ op=par(no.readonly=TRUE)
+ on.exit(par(op))
+ if(is.null(attributes(X)$dim)) X=as.matrix(X)
+ nx=ncol(X)
+ if(nrow(X) < 100) {cat("fewer than 100 draws submitted \n"); return(invisible())}
+ if(!missing(tvalues)){
+ if(mode(tvalues) !="numeric") {stop("tvalues must be a numeric vector \n")}
+ else
+ {if(length(tvalues)!=nx) stop("tvalues are wrong length \n")}
+ }
+ if(nx==1) par(mfrow=c(1,1))
+ if(nx==2) par(mfrow=c(2,1))
+ if(nx==3) par(mfrow=c(3,1))
+ if(nx==4) par(mfrow=c(2,2))
+ if(nx>=5) par(mfrow=c(3,2))
+
+ if(missing(names)) {names=as.character(1:nx)}
+ if (DEN) ylabtxt="density" else ylabtxt="freq"
+ for(index in 1:nx){
+ hist(X[(burnin+1):nrow(X),index],xlab="",ylab=ylabtxt,main=names[index],freq=!DEN,col="magenta",...)
+ if(!missing(tvalues)) abline(v=tvalues[index],lwd=2,col="blue")
+ if(INT){
+ quants=quantile(X[(burnin+1):nrow(X),index],prob=c(.025,.975))
+ mean=mean(X[(burnin+1):nrow(X),index])
+ semean=numEff(X[(burnin+1):nrow(X),index])$stderr
+ text(quants[1],0,"|",cex=3.0,col="green")
+ text(quants[2],0,"|",cex=3.0,col="green")
+ text(mean,0,"|",cex=3.0,col="red")
+ text(mean-2*semean,0,"|",cex=2,col="yellow")
+ text(mean+2*semean,0,"|",cex=2,col="yellow")
+ }
+ par(ask=dev.interactive())
+ }
+ if(TRACEPLOT){
+ if(nx==1) par(mfrow=c(1,2))
+ if(nx==2) par(mfrow=c(2,2))
+ if(nx>=3) par(mfrow=c(3,2))
+ for(index in 1:nx){
+ plot(as.vector(X[,index]),xlab="",ylab="",main=names[index],type="l",col="red")
+ if(!missing(tvalues)) abline(h=tvalues[index],lwd=2,col="blue")
+ if(var(X[,index])>1.0e-20) {acf(as.vector(X[,index]),xlab="",ylab="",main="")}
+ else
+ {plot.default(X[,index],xlab="",ylab="",type="n",main="No ACF Produced")}
+ }
+ }
+ invisible()
+}
+
diff --git a/R/plot.bayesm.nmix.R b/R/plot.bayesm.nmix.R
new file mode 100755
index 0000000..6c498bd
--- /dev/null
+++ b/R/plot.bayesm.nmix.R
@@ -0,0 +1,80 @@
+plot.bayesm.nmix=function(x,names,burnin=trunc(.1*nrow(probdraw)),Grid,bi.sel,nstd=2,...){
+#
+# S3 method to plot normal mixture marginal and bivariate densities
+# nmixlist is a list of 3 components, nmixlist[[1]]: array of mix comp prob draws,
+# mmixlist[[2]] is not used, nmixlist[[3]] is list of draws of components
+# P. Rossi 2/08
+#
+ nmixlist=x
+ if(mode(nmixlist) != "list") stop(" Argument must be a list \n")
+ probdraw=nmixlist[[1]]; compdraw=nmixlist[[3]]
+ if(!is.matrix(probdraw)) stop(" First element of list (probdraw) must be a matrix \n")
+ if(mode(compdraw) != "list") stop(" Third element of list (compdraw) must be a list \n")
+ op=par(no.readonly=TRUE)
+ on.exit(par(op))
+ R=nrow(probdraw)
+ if(R < 100) {cat(" fewer than 100 draws submitted \n"); return(invisible())}
+ datad=length(compdraw[[1]][[1]]$mu)
+ if(missing(bi.sel)) bi.sel=list(c(1,2)) # default to the first pair of variables
+ ind=as.integer(seq(from=(burnin+1),to=R,length.out=max(200,trunc(.05*R))))
+ if(missing(Grid)){
+ out=momMix(probdraw[ind,],compdraw[ind])
+ mu=out$mu
+ sd=out$sd
+ Grid=matrix(0,nrow=50,ncol=datad)
+ for(i in 1:datad ) Grid[,i]=c(seq(from=(mu[i]-nstd*sd[i]),to=(mu[i]+nstd*sd[i]),length=50))
+ }
+ #
+ # plot posterior mean of marginal densities
+ #
+ mden=eMixMargDen(Grid,probdraw[ind,],compdraw[ind])
+ nx=datad
+ if(nx==1) par(mfrow=c(1,1))
+ if(nx==2) par(mfrow=c(2,1))
+ if(nx==3) par(mfrow=c(3,1))
+ if(nx==4) par(mfrow=c(2,2))
+ if(nx>=5) par(mfrow=c(3,2))
+
+ if(missing(names)) {names=as.character(1:nx)}
+ for(index in 1:nx){
+ if(index == 2) par(ask=dev.interactive())
+ plot(range(Grid[,index]),c(0,1.1*max(mden[,index])),type="n",xlab="",ylab="density")
+ title(names[index])
+ lines(Grid[,index],mden[,index],col="black",lwd=1.5)
+ polygon(c(Grid[1,index],Grid[,index],Grid[nrow(Grid),index]),c(0,mden[,index],0),col="magenta")
+ }
+ #
+ # now plot bivariates in list bi.sel
+ #
+ par(ask=dev.interactive())
+ nsel=length(bi.sel)
+ ngrid=50
+ den=array(0,dim=c(nsel,ngrid,ngrid))
+ lstxixj=NULL
+ for(sel in 1:nsel){
+ i=bi.sel[[sel]][1]
+ j=bi.sel[[sel]][2]
+ rxi=range(Grid[,i])
+ rxj=range(Grid[,j])
+ xi=seq(from=rxi[1],to=rxi[2],length.out=ngrid)
+ xj=seq(from=rxj[1],to=rxj[2],length.out=ngrid)
+ lstxixj[[sel]]=list(xi,xj)
+ for(elt in ind){
+ den[sel,,]=den[sel,,]+mixDenBi(i,j,xi,xj,probdraw[elt,],compdraw[[elt]])
+ }
+ }
+ nx=nsel
+ par(mfrow=c(1,1))
+
+ for(index in 1:nx){
+ xi=unlist(lstxixj[[index]][1])
+ xj=unlist(lstxixj[[index]][2])
+ xlabtxt=paste("Var ",bi.sel[[index]][1],sep="")
+ ylabtxt=paste("Var ",bi.sel[[index]][2],sep="")
+ image(xi,xj,den[index,,],col=terrain.colors(100),xlab=xlabtxt,ylab=ylabtxt)
+ contour(xi,xj,den[index,,],add=TRUE,drawlabels=FALSE)
+ }
+
+ invisible()
+}
+
diff --git a/R/rbiNormGibbs.R b/R/rbiNormGibbs.R
index 7607c33..695246f 100755
--- a/R/rbiNormGibbs.R
+++ b/R/rbiNormGibbs.R
@@ -112,5 +112,7 @@ if(continue != "n" & continue != "No" & continue != "no"){
title(paste("IID draws: Rho =",rho))
points(idraws,pch=20,col="magenta",cex=.7)
}
+attributes(draws)$class=c("bayesm.mat","mcmc")
+attributes(draws)$mcpar=c(1,R,1)
return(draws)
}
diff --git a/R/rbprobitGibbs.R b/R/rbprobitGibbs.R
index c951916..6a4010d 100755
--- a/R/rbprobitGibbs.R
+++ b/R/rbprobitGibbs.R
@@ -4,6 +4,8 @@ function(Data,Prior,Mcmc)
#
# revision history:
# p. rossi 1/05
+# 3/07 added validity check of values of y and classes
+# 3/07 fixed error with betabar supplied
#
# purpose:
# draw from posterior for binary probit using Gibbs Sampler
@@ -69,6 +71,7 @@ nobs=length(y)
# check data for validity
#
if(length(y) != nrow(X) ) {pandterm("y and X not of same row dim")}
+if(sum(unique(y) %in% c(0:1)) < length(unique(y))) {pandterm("Invalid y, must be 0,1")}
#
# check for Prior
#
@@ -77,7 +80,7 @@ if(missing(Prior))
else
{
if(is.null(Prior$betabar)) {betabar=c(rep(0,nvar))}
- else {Deltabar=Prior$Deltabar}
+ else {betabar=Prior$betabar}
if(is.null(Prior$A)) {A=.01*diag(nvar)}
else {A=Prior$A}
}
@@ -103,6 +106,9 @@ else
#
cat(" ", fill=TRUE)
cat("Starting Gibbs Sampler for Binary Probit Model",fill=TRUE)
+cat(" with ",length(y)," observations",fill=TRUE)
+cat("Table of y Values",fill=TRUE)
+print(table(y))
cat(" ", fill=TRUE)
cat("Prior Parms: ",fill=TRUE)
cat("betabar",fill=TRUE)
@@ -148,5 +154,7 @@ for (rep in 1:R)
}
ctime = proc.time()[3]
cat(' Total Time Elapsed: ',round((ctime-itime)/60,2),'\n')
+attributes(betadraw)$class=c("bayesm.mat","mcmc")
+attributes(betadraw)$mcpar=c(1,R,keep)
return(list(betadraw=betadraw))
}
diff --git a/R/rhierBinLogit.R b/R/rhierBinLogit.R
index a778b7c..b501dfb 100755
--- a/R/rhierBinLogit.R
+++ b/R/rhierBinLogit.R
@@ -6,6 +6,8 @@ function(Data,Prior,Mcmc){
#
# revision history:
# changed 5/12/05 by Rossi to add error checking
+# 1/07 removed init.rmultiregfp
+# 3/07 added classes
#
# purpose: run binary heterogeneous logit model
#
@@ -17,7 +19,7 @@ function(Data,Prior,Mcmc){
# Z is a matrix of demographic variables nlgt*nz that have been
# mean centered so that the intercept is interpretable
# Prior contains a list of (nu,V,Deltabar,ADelta)
-# beta_i ~ N(Delta,Vbeta)
+# beta_i ~ N(Z%*%Delta,Vbeta)
# vec(Delta) ~ N(vec(Deltabar),Vbeta (x) ADelta^-1)
# Vbeta ~ IW(nu,V)
# Mcmc is a list of (sbeta,R,keep)
@@ -158,11 +160,6 @@ betan = array(0,dim=c(nvar))
reject = array(0,dim=c(R/keep))
llike=array(0,dim=c(R/keep))
-#
-# set up fixed parm for the draw of Vbeta, Delta=Delta
-#
-Fparm=init.rmultiregfp(Z,ADelta,Deltabar,nu,V)
-
itime=proc.time()[3]
cat("MCMC Iteration (est time to end - min)",fill=TRUE)
fsh()
@@ -195,7 +192,7 @@ logkold = -.5*(t(betad)-Z[i,]%*%oldDelta) %*% oldVbetai %*% (betad-t(Z[i,]%*%old
rej = rej+1 }
}
# Draw B-bar and V as a multivariate regression
- out=rmultiregfp(oldbetas,Z,Fparm)
+ out=rmultireg(oldbetas,Z,Deltabar,ADelta,nu,V)
oldDelta=out$B
oldVbeta=out$Sigma
oldVbetai=chol2inv(chol(oldVbeta))
@@ -218,6 +215,13 @@ logkold = -.5*(t(betad)-Z[i,]%*%oldDelta) %*% oldVbetai %*% (betad-t(Z[i,]%*%old
ctime=proc.time()[3]
cat(" Total Time Elapsed: ",round((ctime-itime)/60,2),fill=TRUE)
+
+attributes(betadraw)$class=c("bayesm.hcoef")
+attributes(Deltadraw)$class=c("bayesm.mat","mcmc")
+attributes(Deltadraw)$mcpar=c(1,R,keep)
+attributes(Vbetadraw)$class=c("bayesm.var","bayesm.mat","mcmc")
+attributes(Vbetadraw)$mcpar=c(1,R,keep)
+
return(list(betadraw=betadraw,Vbetadraw=Vbetadraw,Deltadraw=Deltadraw,llike=llike,reject=reject))
}
diff --git a/R/rhierLinearMixture.R b/R/rhierLinearMixture.R
index 94fa346..e4acf77 100755
--- a/R/rhierLinearMixture.R
+++ b/R/rhierLinearMixture.R
@@ -6,6 +6,7 @@ function(Data,Prior,Mcmc)
# changed 12/17/04 by rossi to fix bug in drawdelta when there is zero/one unit
# in a mixture component
# adapted to linear model by Vicky Chen 6/06
+# put in classes 3/07
#
# purpose: run hierarchical linear model with mixture of normals
#
@@ -300,7 +301,6 @@ for(rep in 1:R)
tau[reg]=regout$sigmasqdraw
}
#
- #
# print time to completion and draw # every 100th draw
#
if(((rep/100)*100) ==(floor(rep/100)*100))
@@ -322,8 +322,16 @@ for(rep in 1:R)
}
ctime=proc.time()[3]
cat(" Total Time Elapsed: ",round((ctime-itime)/60,2),fill=TRUE)
+attributes(taudraw)$class=c("bayesm.mat","mcmc")
+attributes(taudraw)$mcpar=c(1,R,keep)
+if(drawdelta){
+ attributes(Deltadraw)$class=c("bayesm.mat","mcmc")
+ attributes(Deltadraw)$mcpar=c(1,R,keep)}
+attributes(betadraw)$class=c("bayesm.hcoef")
+nmix=list(probdraw=probdraw,zdraw=NULL,compdraw=compdraw)
+attributes(nmix)$class="bayesm.nmix"
if(drawdelta)
- {return(list(taudraw=taudraw,Deltadraw=Deltadraw,betadraw=betadraw,probdraw=probdraw,compdraw=compdraw))}
+ {return(list(taudraw=taudraw,Deltadraw=Deltadraw,betadraw=betadraw,nmix=nmix))}
else
- {return(list(taudraw=taudraw,betadraw=betadraw,probdraw=probdraw,compdraw=compdraw))}
+ {return(list(taudraw=taudraw,betadraw=betadraw,nmix=nmix))}
}
diff --git a/R/rhierLinearModel.R b/R/rhierLinearModel.R
index 8f164bb..23d12b0 100755
--- a/R/rhierLinearModel.R
+++ b/R/rhierLinearModel.R
@@ -5,6 +5,7 @@ function(Data,Prior,Mcmc)
# Revision History
# 1/17/05 P. Rossi
# 10/05 fixed error in setting prior if Prior argument is missing
+# 3/07 added classes
#
# Purpose:
# run hiearchical regression model
@@ -204,8 +205,6 @@ betas=matrix(double(nreg*nvar),ncol=nvar)
# "k" = nz
# general model: Beta = Z Delta + U
#
-Fparm=init.rmultiregfp(Z,A,Deltabar,nu,V)
-#
# Create XpX elements of regdata and initialize tau
#
regdata=lapply(regdata,append)
@@ -235,7 +234,7 @@ for(rep in 1:R)
#
# draw Vbeta, Delta | {beta_i}
#
- rmregout=rmultiregfp(betas,Z,Fparm)
+ rmregout=rmultireg(betas,Z,Deltabar,A,nu,V)
Vbeta=as.vector(rmregout$Sigma)
Delta=as.vector(rmregout$B)
#
@@ -259,5 +258,13 @@ for(rep in 1:R)
ctime = proc.time()[3]
cat(' Total Time Elapsed: ',round((ctime-itime)/60,2),'\n')
+attributes(taudraw)$class=c("bayesm.mat","mcmc")
+attributes(taudraw)$mcpar=c(1,R,keep)
+attributes(Deltadraw)$class=c("bayesm.mat","mcmc")
+attributes(Deltadraw)$mcpar=c(1,R,keep)
+attributes(Vbetadraw)$class=c("bayesm.var","bayesm.mat","mcmc")
+attributes(Vbetadraw)$mcpar=c(1,R,keep)
+attributes(betadraw)$class=c("bayesm.hcoef")
+
return(list(Vbetadraw=Vbetadraw,Deltadraw=Deltadraw,betadraw=betadraw,taudraw=taudraw))
}
diff --git a/R/rhierMnlRwMixture.R b/R/rhierMnlRwMixture.R
index 397243e..062346e 100755
--- a/R/rhierMnlRwMixture.R
+++ b/R/rhierMnlRwMixture.R
@@ -7,6 +7,7 @@ function(Data,Prior,Mcmc)
# in a mixture component
# added loglike output, changed to reflect new argument order in llmnl, mnlHess 9/05
# changed weighting scheme to (1-w)logl_i + w*Lbar (normalized) 12/05
+# 3/07 added classes
#
# purpose: run hierarchical mnl logit model with mixture of normals
# using RW and cov(RW inc) = (hess_i + Vbeta^-1)^-1
@@ -373,8 +374,14 @@ for(rep in 1:R)
}
ctime=proc.time()[3]
cat(" Total Time Elapsed: ",round((ctime-itime)/60,2),fill=TRUE)
+if(drawdelta){
+ attributes(Deltadraw)$class=c("bayesm.mat","mcmc")
+ attributes(Deltadraw)$mcpar=c(1,R,keep)}
+attributes(betadraw)$class=c("bayesm.hcoef")
+nmix=list(probdraw=probdraw,zdraw=NULL,compdraw=compdraw)
+attributes(nmix)$class="bayesm.nmix"
if(drawdelta)
- {return(list(Deltadraw=Deltadraw,betadraw=betadraw,probdraw=probdraw,compdraw=compdraw,loglike=loglike))}
+ {return(list(Deltadraw=Deltadraw,betadraw=betadraw,nmix=nmix,loglike=loglike))}
else
- {return(list(betadraw=betadraw,probdraw=probdraw,compdraw=compdraw,loglike=loglike))}
+ {return(list(betadraw=betadraw,nmix=nmix,loglike=loglike))}
}
diff --git a/R/rhierNegbinRw.R b/R/rhierNegbinRw.R
index 9050b14..3761ca9 100755
--- a/R/rhierNegbinRw.R
+++ b/R/rhierNegbinRw.R
@@ -5,6 +5,7 @@ function(Data, Prior, Mcmc) {
# Sridhar Narayanan - 05/2005
# P. Rossi 6/05
# fixed error with nobs not specified and changed llnegbinFract 9/05
+# 3/07 added classes
#
# Model
# (y_i|lambda_i,alpha) ~ Negative Binomial(Mean = lambda_i, Overdispersion par = alpha)
@@ -49,8 +50,12 @@ llnegbin =
function(par,X,y, nvar) {
# Computes the log-likelihood
beta = par[1:nvar]
- alpha = exp(par[nvar+1])
- ll = sum(log(dnbinom(y,size=alpha,mu=exp(X%*%beta))))
+ alpha = exp(par[nvar+1])+1.0e-50
+ mean=exp(X%*%beta)
+ prob=alpha/(alpha+mean)
+ prob=ifelse(prob<1.0e-100,1.0e-100,prob)
+ out=.Internal(dnbinom(y,alpha,prob,TRUE))
+ return(sum(out))
}
llnegbinFract =
@@ -310,6 +315,14 @@ for (r in 1:R)
}
}
ctime = proc.time()[3]
+
+attributes(alphadraw)$class=c("bayesm.mat","mcmc")
+attributes(alphadraw)$mcpar=c(1,R,keep)
+attributes(Deltadraw)$class=c("bayesm.mat","mcmc")
+attributes(Deltadraw)$mcpar=c(1,R,keep)
+attributes(Vbetadraw)$class=c("bayesm.var","bayesm.mat","mcmc")
+attributes(Vbetadraw)$mcpar=c(1,R,keep)
+attributes(Betadraw)$class=c("bayesm.hcoef")
cat(' Total Time Elapsed: ',round((ctime-itime)/60,2),'\n')
return(list(llike=llike,Betadraw=Betadraw,alphadraw=alphadraw, Vbetadraw=Vbetadraw, Deltadraw=Deltadraw,
diff --git a/R/rivGibbs.R b/R/rivGibbs.R
index 92e6617..a2c4f68 100755
--- a/R/rivGibbs.R
+++ b/R/rivGibbs.R
@@ -7,6 +7,7 @@ function(Data,Prior,Mcmc)
# p. rossi 3/05
# p. rossi 1/06 -- fixed error in nins
# p. rossi 1/06 -- fixed def Prior settings for nu,V
+# 3/07 added classes
#
# purpose:
# draw from posterior for linear I.V. model
@@ -210,5 +211,18 @@ for(rep in 1:R) {
}
}
+ctime = proc.time()[3]
+cat(' Total Time Elapsed: ',round((ctime-itime)/60,2),'\n')
+
+attributes(deltadraw)$class=c("bayesm.mat","mcmc")
+attributes(deltadraw)$mcpar=c(1,R,keep)
+attributes(betadraw)$class=c("bayesm.mat","mcmc")
+attributes(betadraw)$mcpar=c(1,R,keep)
+attributes(gammadraw)$class=c("bayesm.mat","mcmc")
+attributes(gammadraw)$mcpar=c(1,R,keep)
+attributes(Sigmadraw)$class=c("bayesm.var","bayesm.mat","mcmc")
+attributes(Sigmadraw)$mcpar=c(1,R,keep)
+
+
return(list(deltadraw=deltadraw,betadraw=betadraw,gammadraw=gammadraw,Sigmadraw=Sigmadraw))
}
diff --git a/R/rmnlIndepMetrop.R b/R/rmnlIndepMetrop.R
index 0682bc4..a917124 100755
--- a/R/rmnlIndepMetrop.R
+++ b/R/rmnlIndepMetrop.R
@@ -86,8 +86,10 @@ else
#
cat(" ", fill=TRUE)
cat("Starting Independence Metropolis Sampler for Multinomial Logit Model",fill=TRUE)
-cat(" with ",p," alternatives",fill=TRUE)
+cat(" ",length(y)," obs with ",p," alternatives",fill=TRUE)
cat(" ", fill=TRUE)
+cat("Table of y Values",fill=TRUE)
+print(table(y))
cat("Prior Parms: ",fill=TRUE)
cat("betabar",fill=TRUE)
print(betabar)
@@ -157,5 +159,7 @@ for (rep in 1:R)
}
ctime = proc.time()[3]
cat(' Total Time Elapsed: ',round((ctime-itime)/60,2),'\n')
+attributes(betadraw)$class=c("bayesm.mat","mcmc")
+attributes(betadraw)$mcpar=c(1,R,keep)
return(list(betadraw=betadraw,loglike=loglike,acceptr=naccept/R))
}
diff --git a/R/rmnpGibbs.R b/R/rmnpGibbs.R
index e687257..ab90870 100755
--- a/R/rmnpGibbs.R
+++ b/R/rmnpGibbs.R
@@ -3,7 +3,8 @@ function(Data,Prior,Mcmc)
{
#
# Revision History:
-# modified by rossi 12/18/04 to include error checking`
+# modified by rossi 12/18/04 to include error checking
+# 3/07 added classes
#
# purpose: Gibbs MNP model with full covariance matrix
#
@@ -94,6 +95,8 @@ cat("Starting Gibbs Sampler for MNP",fill=TRUE)
cat(" ",n," obs; ",p," choice alternatives; ",k," indep vars (including intercepts)",fill=TRUE)
cat(" ",R," reps; keeping every ",keep,"th draw",fill=TRUE)
cat(" ",fill=TRUE)
+cat("Table of y values",fill=TRUE)
+print(table(y))
cat("Prior Parms:",fill=TRUE)
cat("betabar",fill=TRUE)
print(betabar)
@@ -192,6 +195,12 @@ for (rep in 1:R)
wold=wnew
betaold=betanew
}
+ctime = proc.time()[3]
+cat(' Total Time Elapsed: ',round((ctime-itime)/60,2),'\n')
+attributes(betadraw)$class=c("bayesm.mat","mcmc")
+attributes(betadraw)$mcpar=c(1,R,keep)
+attributes(sigmadraw)$class=c("bayesm.var","bayesm.mat","mcmc")
+attributes(sigmadraw)$mcpar=c(1,R,keep)
list(betadraw=betadraw,sigmadraw=sigmadraw)
}
diff --git a/R/rmultiregfp.R b/R/rmultiregfp.R
deleted file mode 100755
index d7fb186..0000000
--- a/R/rmultiregfp.R
+++ /dev/null
@@ -1,53 +0,0 @@
-rmultiregfp=
-function(Y,X,Fparm)
-{
-#
-# revision history:
-# changed 1/11/05 by P. Rossi to fix sum of squares error
-#
-# purpose:
-# draw from posterior for Multivariate Regression Model
-# arguments:
-# Y is n x m matrix
-# X is n x k
-# Fparm is a list of fixed parameters created by init.rmultiregfp
-# output:
-# list of B, Sigma draws of matrix of coefficients and Sigma matrix
-# model:
-# Y=XB+U cov(u_i) = Sigma
-# B is k x m matrix of coefficients
-# priors: beta|Sigma ~ N(betabar,Sigma (x) A^-1)
-# betabar=vec(Bbar)
-# beta = vec(B)
-# Sigma ~ IW(nu,V) or Sigma^-1 ~ W(nu, V^-1)
-n=nrow(Y)
-m=ncol(Y)
-k=ncol(X)
-#
-# first draw Sigma
-#
-W=rbind(X,Fparm$RA)
-Z=rbind(Y,Fparm$RABbar)
-# note: Y,X,A,Bbar must be matrices!
-IR=Fparm$IR
-# W'W = R'R & (W'W)^-1 = IRIR' -- this is the UL decomp!
-Btilde=crossprod(t(IR))%*%crossprod(W,Z)
-# IRIR'(W'Z) = (X'X+A)^-1(X'Y + ABbar)
-S=crossprod(Z-W%*%Btilde)
-# E'E
-rwout=rwishart(Fparm$nu+n,chol2inv(chol(Fparm$V+S)))
-#
-# now draw B given Sigma
-# note beta ~ N(vec(Btilde),Sigma (x) Covxxa)
-# Cov=(X'X + A)^-1 = IR t(IR)
-# Sigma=CICI'
-# therefore, cov(beta)= Omega = CICI' (x) IR IR' = (CI (x) IR) (CI (x) IR)'
-# so to draw beta we do beta= vec(Btilde) +(CI (x) IR)vec(Z_mk)
-# Z_mk is m x k matrix of N(0,1)
-# since vec(ABC) = (C' (x) A)vec(B), we have
-# B = Btilde + IR Z_mk CI'
-#
-B = Btilde + IR%*%matrix(rnorm(m*k),ncol=m)%*%t(rwout$CI)
-return(list(B=B,Sigma=rwout$IW))
-}
-
diff --git a/R/rmvpGibbs.R b/R/rmvpGibbs.R
index 3da0732..dd27f36 100755
--- a/R/rmvpGibbs.R
+++ b/R/rmvpGibbs.R
@@ -4,6 +4,7 @@ function(Data,Prior,Mcmc)
#
# Revision History:
# modified by rossi 12/18/04 to include error checking
+# 3/07 added classes
#
# purpose: Gibbs MVP model with full covariance matrix
#
@@ -189,6 +190,13 @@ for (rep in 1:R)
wold=wnew
betaold=betanew
}
+ctime = proc.time()[3]
+cat(' Total Time Elapsed: ',round((ctime-itime)/60,2),'\n')
+
+attributes(betadraw)$class=c("bayesm.mat","mcmc")
+attributes(betadraw)$mcpar=c(1,R,keep)
+attributes(sigmadraw)$class=c("bayesm.var","bayesm.mat","mcmc")
+attributes(sigmadraw)$mcpar=c(1,R,keep)
return(list(betadraw=betadraw,sigmadraw=sigmadraw))
}
diff --git a/R/rnegbinRw.R b/R/rnegbinRw.R
index 0b7d36d..79cf0a1 100755
--- a/R/rnegbinRw.R
+++ b/R/rnegbinRw.R
@@ -4,6 +4,7 @@ function(Data, Prior, Mcmc) {
# Revision History
# Sridhar Narayanan - 05/2005
# P. Rossi 6/05
+# 3/07 added classes
#
# Model
# (y|lambda,alpha) ~ Negative Binomial(Mean = lambda, Overdispersion par = alpha)
@@ -35,13 +36,16 @@ function(Data, Prior, Mcmc) {
#
# Definitions of functions used within rhierNegbinRw
#
-
llnegbin =
function(par,X,y, nvar) {
# Computes the log-likelihood
beta = par[1:nvar]
- alpha = exp(par[nvar+1])
- ll = sum(log(dnbinom(y,size=alpha,mu=exp(X%*%beta))))
+ alpha = exp(par[nvar+1])+1.0e-50
+ mean=exp(X%*%beta)
+ prob=alpha/(alpha+mean)
+ prob=ifelse(prob<1.0e-100,1.0e-100,prob)
+ out=.Internal(dnbinom(y,alpha,prob,TRUE))
+ return(sum(out))
}
lpostbetai =
@@ -203,8 +207,13 @@ if(r%%keep == 0) {
}
ctime = proc.time()[3]
-
cat(' Total Time Elapsed: ',round((ctime-itime)/60,2),'\n')
+
+
+attributes(betadraw)$class=c("bayesm.mat","mcmc")
+attributes(betadraw)$mcpar=c(1,R,keep)
+attributes(alphadraw)$class=c("bayesm.mat","mcmc")
+attributes(alphadraw)$mcpar=c(1,R,keep)
return(list(llike=llike,betadraw=betadraw,alphadraw=alphadraw,
acceptrbeta=nacceptbeta/R*100,acceptralpha=nacceptalpha/R*100))
}
diff --git a/R/rnmixGibbs.R b/R/rnmixGibbs.R
index cb2ed78..e1fbf0b 100755
--- a/R/rnmixGibbs.R
+++ b/R/rnmixGibbs.R
@@ -7,6 +7,7 @@ function(Data,Prior,Mcmc)
# add check to see if Mubar is a vector 9/05
# fixed bug in saving comps draw comps[[mkeep]]= 9/05
# fixed so that ncomp can be =1; added check that nobs >= 2*ncomp 12/06
+# 3/07 added classes
#
# purpose: do Gibbs sampling inference for a mixture of multivariate normals
#
@@ -164,5 +165,8 @@ for(rep in 1:R)
}
ctime = proc.time()[3]
cat(' Total Time Elapsed: ',round((ctime-itime)/60,2),'\n')
-return(list(probdraw=pdraw,zdraw=zdraw,compdraw=compdraw))
+
+nmix=list(probdraw=pdraw,zdraw=zdraw,compdraw=compdraw)
+attributes(nmix)$class="bayesm.nmix"
+return(nmix)
}
diff --git a/R/rordprobitGibbs.R b/R/rordprobitGibbs.R
new file mode 100755
index 0000000..f505d2a
--- /dev/null
+++ b/R/rordprobitGibbs.R
@@ -0,0 +1,290 @@
+rordprobitGibbs=
+ function(Data,Prior,Mcmc)
+{
+#
+# revision history:
+# 3/07 Hsiu-Wen Liu
+# 3/07 fixed naming of dstardraw rossi
+#
+# purpose:
+# draw from posterior for ordered probit using Gibbs Sampler
+# and metropolis RW
+#
+# Arguments:
+# Data - list of X,y,k
+# X is nobs x nvar, y is nobs vector of 1,2,.,k (ordinal variable)
+# Prior - list of A, betabar
+# A is nvar x nvar prior preci matrix
+# betabar is nvar x 1 prior mean
+# Ad is ndstar x ndstar prior preci matrix of dstar (ncut is number of cut-offs being estimated)
+# dstarbar is ndstar x 1 prior mean of dstar
+# Mcmc
+# R is number of draws
+# keep is thinning parameter
+# s is scale parameter of random work Metropolis
+#
+# Output:
+# list of betadraws and cutdraws
+#
+# Model:
+# z=Xbeta + e < 0 e ~N(0,1)
+# y=1,..,k, if z~c(c[k], c[k+1])
+#
+# cutoffs = c[1],..,c[k+1]
+# dstar = dstar[1],dstar[k-2]
+# set c[1]=-100, c[2]=0, ...,c[k+1]=100
+#
+# c[3]=exp(dstar[1]),c[4]=c[3]+exp(dstar[2]),...,
+# c[k]=c[k-1]+exp(datsr[k-2])
+#
+# Note: 1. length of dstar = length of cutoffs - 3
+# 2. Be careful in assessing prior parameter, Ad. .1 is too small for many applications.
+#
+# Prior: beta ~ N(betabar,A^-1)
+# dstar ~ N(dstarbar, Ad^-1)
+#
+#
+# ----------------------------------------------------------------------
+# define functions needed
+#
+breg1=
+function(root,X,y,Abetabar)
+{
+#
+# p.rossi 12/04
+#
+# Purpose: draw from posterior for linear regression, sigmasq=1.0
+#
+# Arguments:
+# root is chol((X'X+A)^-1)
+# Abetabar = A*betabar
+#
+# Output: draw from posterior
+#
+# Model: y = Xbeta + e e ~ N(0,I)
+# Prior: beta ~ N(betabar,A^-1)
+#
+cov=crossprod(root,root)
+betatilde=cov%*%(crossprod(X,y)+Abetabar)
+betatilde+t(root)%*%rnorm(length(betatilde))
+}
+
+#
+# dstartoc is a fuction to transfer dstar to its cut-off value
+
+ dstartoc=function(dstar) {c(-100, 0, cumsum(exp(dstar)), 100)}
+
+# compute conditional likelihood of data given cut-offs
+#
+ lldstar=function(dstar,y,mu){
+ gamma=dstartoc(dstar)
+ arg = pnorm(gamma[y+1]-mu)-pnorm(gamma[y]-mu)
+ epsilon=1.0e-50
+ arg=ifelse(arg < epsilon,epsilon,arg)
+ return(sum(log(arg)))
+ }
+
+
+dstarRwMetrop=
+function(y,mu,olddstar,s,inc.root,dstarbar,oldll,rootdi){
+#
+# function to execute rw metropolis for the dstar
+# y is n vector with element = 1,...,j
+# X is n x k matrix of x values
+# RW increments are N(0,s^2*t(inc.root)%*%inc.root)
+# prior on dstar is N(dstarbar,Sigma) Sigma^-1=rootdi*t(rootdi)
+# inc.root, rootdi are upper triangular
+# this means that we are using the UL decomp of Sigma^-1 for prior
+# olddstar is the current
+
+ stay=0
+ dstarc=olddstar + s*t(inc.root)%*%(matrix(rnorm(ncut),ncol=1))
+ cll=lldstar(dstarc,y,mu)
+ clpost=cll+lndMvn(dstarc,dstarbar,rootdi)
+ ldiff=clpost-oldll-lndMvn(olddstar,dstarbar,rootdi)
+ alpha=min(1,exp(ldiff))
+
+ if(alpha < 1) {unif=runif(1)} else {unif=0}
+ if (unif <= alpha)
+ {dstardraw=dstarc; oldll=cll}
+ else
+ {dstardraw=olddstar; stay=1}
+
+return(list(dstardraw=dstardraw,oldll=oldll, stay=stay))
+}
+
+pandterm=function(message) {stop(message,call.=FALSE)}
+#
+# ----------------------------------------------------------------------
+#
+# check arguments
+#
+if(missing(Data)) {pandterm("Requires Data argument -- list of y and X")}
+ if(is.null(Data$X)) {pandterm("Requires Data element X")}
+ X=Data$X
+ if(is.null(Data$y)) {pandterm("Requires Data element y")}
+ y=Data$y
+ if(is.null(Data$k)) {pandterm("Requires Data element k")}
+ k=Data$k
+
+nvar=ncol(X)
+nobs=length(y)
+ndstar = k-2 # number of dstar being estimated
+ncuts = k+1 # number of cut-offs (including zero and two ends)
+ncut = ncuts-3 # number of cut-offs being estimated c[1]=-100, c[2]=0, c[k+1]=100
+
+#
+# check data for validity
+#
+if(length(y) != nrow(X) ) {pandterm("y and X not of same row dim")}
+if( sum(unique(y) %in% (1:k) ) < length(unique(y)) )
+ {pandterm("some value of y is not vaild")}
+
+#
+
+#
+# check for Prior
+#
+if(missing(Prior))
+ { betabar=c(rep(0,nvar)); A=.01*diag(nvar); Ad=diag(ndstar); dstarbar=c(rep(0,ndstar))}
+else
+ {
+ if(is.null(Prior$betabar)) {betabar=c(rep(0,nvar))}
+ else {betabar=Prior$betabar}
+ if(is.null(Prior$A)) {A=.01*diag(nvar)}
+ else {A=Prior$A}
+ if(is.null(Prior$Ad)) {Ad=diag(ndstar)}
+ else {Ad=Prior$Ad}
+ if(is.null(Prior$dstarbar)) {dstarbar=c(rep(0,ndstar))}
+ else {dstarbar=Prior$dstarbar}
+ }
+#
+# check dimensions of Priors
+#
+
+if(ncol(A) != nrow(A) || ncol(A) != nvar || nrow(A) != nvar)
+ {pandterm(paste("bad dimensions for A",dim(A)))}
+if(length(betabar) != nvar)
+ {pandterm(paste("betabar wrong length, length= ",length(betabar)))}
+if(ncol(Ad) != nrow(Ad) || ncol(Ad) != ndstar || nrow(Ad) != ndstar)
+ {pandterm(paste("bad dimensions for Ad",dim(Ad)))}
+if(length(dstarbar) != ndstar)
+ {pandterm(paste("dstarbar wrong length, length= ",length(dstarbar)))}
+
+#
+# check MCMC argument
+#
+if(missing(Mcmc)) {pandterm("requires Mcmc argument")}
+else
+ {
+ if(is.null(Mcmc$R))
+ {pandterm("requires Mcmc element R")} else {R=Mcmc$R}
+ if(is.null(Mcmc$keep)) {keep=1} else {keep=Mcmc$keep}
+ if(is.null(Mcmc$s)) {s=2.93/sqrt(ndstar)} else {s=Mcmc$s}
+ }
+#
+# print out problem
+#
+cat(" ", fill=TRUE)
+cat("Starting Gibbs Sampler for Ordered Probit Model",fill=TRUE)
+cat(" with ",nobs,"observations",fill=TRUE)
+cat(" ", fill=TRUE)
+cat("Table of y values",fill=TRUE)
+print(table(y))
+cat(" ",fill=TRUE)
+cat("Prior Parms: ",fill=TRUE)
+cat("betabar",fill=TRUE)
+print(betabar)
+cat(" ", fill=TRUE)
+cat("A",fill=TRUE)
+print(A)
+cat(" ", fill=TRUE)
+cat("dstarbar",fill=TRUE)
+print(dstarbar)
+cat(" ", fill=TRUE)
+cat("Ad",fill=TRUE)
+print(Ad)
+cat(" ", fill=TRUE)
+cat("MCMC parms: ",fill=TRUE)
+cat("R= ",R," keep= ",keep,"s= ",s, fill=TRUE)
+cat(" ",fill=TRUE)
+
+betadraw=matrix(double(floor(R/keep)*nvar),ncol=nvar)
+cutdraw=matrix(double(floor(R/keep)*ncuts),ncol=ncuts)
+dstardraw=matrix(double(floor(R/keep)*ndstar),ncol=ndstar)
+staydraw=array(0,dim=c(R/keep))
+
+sigma=c(rep(1,nrow(X)))
+root=chol(chol2inv(chol((crossprod(X,X)+A))))
+Abetabar=crossprod(A,betabar)
+rootdi=chol(chol2inv(chol(Ad)))
+
+# use (-Hessian+Ad)^(-1) evaluated at betahat as the basis of the
+# covariance matrix for the random walk Metropolis increments
+
+ betahat = chol2inv(chol(crossprod(X,X)))%*% crossprod(X,y)
+ dstarini = c(cumsum(c( rep(0.1, ndstar)))) # set initial value for dstar
+ dstarout = optim(dstarini, lldstar, method = "BFGS", hessian=T,
+ control = list(fnscale = -1,maxit=500,
+ reltol = 1e-06, trace=0), mu=X%*%betahat, y=y)
+ inc.root=chol(chol2inv(chol((-dstarout$hessian+Ad)))) # chol((H+Ad)^-1)
+
+# set initial values for MCMC
+
+ olddstar = c(rep(0,ndstar))
+ beta = betahat
+ cutoffs = dstartoc (olddstar)
+ oldll = lldstar(olddstar,y,mu=X%*%betahat)
+
+#
+#
+# start main iteration loop
+#
+itime=proc.time()[3]
+cat("MCMC Iteration (est time to end - min) ",fill=TRUE)
+fsh()
+# print time to completion and draw # every 100th draw
+#
+for (rep in 1:R)
+{
+ # draw z given beta(i-1), sigma, y, cut-offs
+ z = rtrun (X%*%beta, sigma=sigma, a=cutoffs[y] , b=cutoffs[y+1])
+
+ # draw beta given z and rest
+ beta= breg1(root,X,z, Abetabar)
+
+ # draw gamma given z
+ metropout = dstarRwMetrop(y,X%*%beta,olddstar,s,inc.root,dstarbar,oldll,rootdi)
+ olddstar = metropout$dstardraw
+ oldll = metropout$oldll
+ cutoffs = dstartoc (olddstar)
+ stay = metropout$stay
+
+
+# print time to completion and draw # every 100th draw
+
+ if(rep%%100 == 0)
+ {ctime=proc.time()[3]
+ timetoend=((ctime-itime)/rep)*(R-rep)
+ cat(" ",rep," (",round(timetoend/60,1),")",fill=TRUE)
+ fsh()}
+
+ if(rep%%keep == 0)
+ {mkeep=rep/keep; cutdraw[mkeep,]=cutoffs; dstardraw[mkeep,]=olddstar;betadraw[mkeep,]=beta;staydraw[mkeep]=stay }
+
+}
+ accept=1-sum(staydraw)/(R/keep)
+
+ctime = proc.time()[3]
+cat(' Total Time Elapsed: ',round((ctime-itime)/60,2),'\n')
+
+cutdraw=cutdraw[,2:k]
+attributes(cutdraw)$class="bayesm.mat"
+attributes(betadraw)$class="bayesm.mat"
+attributes(dstardraw)$class="bayesm.mat"
+attributes(cutdraw)$mcpar=c(1,R,keep)
+attributes(betadraw)$mcpar=c(1,R,keep)
+attributes(dstardraw)$mcpar=c(1,R,keep)
+
+return(list(cutdraw=cutdraw,betadraw=betadraw, dstardraw=dstardraw, accept=accept))
+}
diff --git a/R/rscaleUsage.R b/R/rscaleUsage.R
index be7798e..5ebdcfe 100755
--- a/R/rscaleUsage.R
+++ b/R/rscaleUsage.R
@@ -5,6 +5,7 @@ function(Data,Prior,Mcmc)
# purpose: run scale-usage mcmc
# draws y,Sigma,mu,tau,sigma,Lambda,e
# R. McCulloch 12/28/04
+# added classes 3/07
#
# arguments:
# Data:
@@ -420,8 +421,24 @@ for(rep in 1:ndpost) {
}
ctime = proc.time()[3]
cat(' Total Time Elapsed: ',round((ctime-itime)/60,2),'\n')
-return(list(Sigmadraw=drSigma,mudraw=drmu,taudraw = drtau,
- sigmadraw=drsigma,Lambdadraw=drLamda,edraw=dre))
+
+R=ndpost
+mudraw=drmu; taudraw=drtau; sigmadraw=drsigma; Lambdadraw=drLamda;
+edraw=dre; Sigmadraw=drSigma
+attributes(mudraw)$class=c("bayesm.mat","mcmc")
+attributes(mudraw)$mcpar=c(1,R,keep)
+attributes(taudraw)$class=c("bayesm.mat","mcmc")
+attributes(taudraw)$mcpar=c(1,R,keep)
+attributes(sigmadraw)$class=c("bayesm.mat","mcmc")
+attributes(sigmadraw)$mcpar=c(1,R,keep)
+attributes(Lambdadraw)$class=c("bayesm.mat","mcmc")
+attributes(Lambdadraw)$mcpar=c(1,R,keep)
+attributes(edraw)$class=c("bayesm.mat","mcmc")
+attributes(edraw)$mcpar=c(1,R,keep)
+attributes(Sigmadraw)$class=c("bayesm.var","bayesm.mat","mcmc")
+attributes(Sigmadraw)$mcpar=c(1,R,keep)
+return(list(Sigmadraw=Sigmadraw,mudraw=mudraw,taudraw = taudraw,
+ sigmadraw=sigmadraw,Lambdadraw=Lambdadraw,edraw=edraw))
}
diff --git a/R/rsurGibbs.R b/R/rsurGibbs.R
index 666cc84..2308611 100755
--- a/R/rsurGibbs.R
+++ b/R/rsurGibbs.R
@@ -4,6 +4,7 @@ function(Data,Prior,Mcmc)
#
# revision history:
# P. Rossi 9/05
+# 3/07 added classes
# Purpose:
# implement Gibbs Sampler for SUR
#
@@ -93,8 +94,7 @@ if(length(betabar) != nvar)
if(missing(Mcmc)) {pandterm("requires Mcmc argument")}
else
{
- if(is.null(Mcmc$R))
- {pandterm("requires Mcmc element R")} else {R=Mcmc$R}
+ if(is.null(Mcmc$R)) {pandterm("requires Mcmc element R")} else {R=Mcmc$R}
if(is.null(Mcmc$keep)) {keep=1} else {keep=Mcmc$keep}
}
#
@@ -178,6 +178,10 @@ for (rep in 1:R)
ctime = proc.time()[3]
cat(' Total Time Elapsed: ',round((ctime-itime)/60,2),'\n')
+attributes(betadraw)$class=c("bayesm.mat","mcmc")
+attributes(betadraw)$mcpar=c(1,R,keep)
+attributes(Sigmadraw)$class=c("bayesm.var","bayesm.mat","mcmc")
+attributes(Sigmadraw)$mcpar=c(1,R,keep)
return(list(betadraw=betadraw,Sigmadraw=Sigmadraw))
}
diff --git a/R/runireg.R b/R/runireg.R
index 5f26e61..7084405 100755
--- a/R/runireg.R
+++ b/R/runireg.R
@@ -5,6 +5,7 @@ function(Data,Prior,Mcmc)
# revision history:
# P. Rossi 1/17/05
# revised 9/05 to put in Data,Prior,Mcmc calling convention
+# 3/07 added classes
# Purpose:
# perform iid draws from posterior of regression model using
# conjugate prior
@@ -140,5 +141,11 @@ beta = btilde + as.vector(sqrt(sigmasq))*IR%*%rnorm(nvar)
}
ctime = proc.time()[3]
cat(' Total Time Elapsed: ',round((ctime-itime)/60,2),'\n')
+
+attributes(betadraw)$class=c("bayesm.mat","mcmc")
+attributes(betadraw)$mcpar=c(1,R,keep)
+attributes(sigmasqdraw)$class=c("bayesm.mat","mcmc")
+attributes(sigmasqdraw)$mcpar=c(1,R,keep)
+
return(list(betadraw=betadraw,sigmasqdraw=sigmasqdraw))
}
diff --git a/R/runiregGibbs.R b/R/runiregGibbs.R
index 83d09ff..ef8020d 100755
--- a/R/runiregGibbs.R
+++ b/R/runiregGibbs.R
@@ -4,6 +4,7 @@ function(Data,Prior,Mcmc)
#
# revision history:
# P. Rossi 1/17/05
+# 3/07 added classes
# Purpose:
# perform Gibbs iterations for Univ Regression Model using
# prior with beta, sigma-sq indep
@@ -138,6 +139,10 @@ for (rep in 1:Mcmc$R)
ctime = proc.time()[3]
cat(' Total Time Elapsed: ',round((ctime-itime)/60,2),'\n')
+attributes(betadraw)$class=c("bayesm.mat","mcmc")
+attributes(betadraw)$mcpar=c(1,R,keep)
+attributes(sigmasqdraw)$class=c("bayesm.mat","mcmc")
+attributes(sigmasqdraw)$mcpar=c(1,R,keep)
return(list(betadraw=betadraw,sigmasqdraw=sigmasqdraw))
}
diff --git a/R/simmnl.R b/R/simmnl.R
deleted file mode 100755
index 32745c7..0000000
--- a/R/simmnl.R
+++ /dev/null
@@ -1,51 +0,0 @@
-simmnl=
-function(p,n,beta)
-{
-#
-# p. rossi 2004
-#
-# Purpose: simulate from MNL (including X values)
-#
-# Arguments:
-# p is number of alternatives
-# n is number of obs
-# beta is true parm value
-#
-# Output:
-# list of X (note: we include full set of intercepts and 2 unif(-1,1) X vars)
-# y (indicator of choice-- 1, ...,p
-# prob is a n x p matrix of choice probs
-#
-# note: first choice alternative has intercept set to zero
-#
-k=length(beta)
-x1=runif(n*p,min=-1,max=1)
-x2=runif(n*p,min=-1,max=1)
-I2=diag(rep(1,p-1))
-zero=rep(0,p-1)
-xadd=rbind(zero,I2)
-for(i in 2:n) {
- xadd=rbind(xadd,zero,I2)
-}
-
-X=cbind(xadd,x1,x2)
-
-# now construct probabilities
-Xbeta=X%*%beta
-p=nrow(Xbeta)/n
-Xbeta=matrix(Xbeta,byrow=T,ncol=p)
-Prob=exp(Xbeta)
-iota=c(rep(1,p))
-denom=Prob%*%iota
-Prob=Prob/as.vector(denom)
-# draw y
-y=vector("double",n)
-ind=1:p
-for (i in 1:n)
- {
- yvec=rmultinom(1,1,Prob[i,])
- y[i]=ind%*%yvec
- }
-
-return(list(y=y,X=X,beta=beta,prob=Prob))
-}
diff --git a/R/simmnlwX.R b/R/simmnlwX.R
deleted file mode 100755
index 4995c26..0000000
--- a/R/simmnlwX.R
+++ /dev/null
@@ -1,41 +0,0 @@
-simmnlwX=
-function(n,X,beta)
-{
-#
-# p. rossi 2004
-#
-# Purpose: simulate from MNL model conditional on X matrix
-#
-# Arguments:
-# j is number of alternatives
-# X is full matrix, including intercepts
-# beta is true parm value
-#
-# Output:
-# list of X (note: we include full set of intercepts and 2 unif(-1,1) X vars)
-# y (indicator of choice-- 1, ...,j
-# prob is a n x j matrix of choice probs
-#
-# note: first choice alternative has intercept set to zero
-#
-#
-k=length(beta)
-# now construct probabilities
-Xbeta=X%*%beta
-j=nrow(Xbeta)/n
-Xbeta=matrix(Xbeta,byrow=T,ncol=j)
-Prob=exp(Xbeta)
-iota=c(rep(1,j))
-denom=Prob%*%iota
-Prob=Prob/as.vector(denom)
-# draw y
-y=vector("double",n)
-ind=1:j
-for (i in 1:n)
- {
- yvec=rmultinom(1,1,Prob[i,])
- y[i]=ind%*%yvec
- }
-
-list(y=y,X=X,beta=beta,prob=Prob)
-}
diff --git a/R/simmnp.R b/R/simmnp.R
deleted file mode 100755
index d9188b9..0000000
--- a/R/simmnp.R
+++ /dev/null
@@ -1,27 +0,0 @@
-simmnp=
-function(X,p,n,beta,sigma) {
-#
-# simulate from MNP: w_i=X_i*beta + e e ~ MVN(0,Sigma)
-# y_i = j if max(w_i) = w_ij
-# if max(w_i) < 0 y_i = p
-# Sigma is p-1 x p-1
-# X is n(p-1) x k
-# beta is k x 1
-#
-# create functions needed
-#
-indmax=function(x) {
-ind=1:length(x)
-ind[x==max(x)]
-}
-
-Xbeta=X%*%beta
-w=as.vector(crossprod(chol(sigma),matrix(rnorm((p-1)*n),ncol=n)))+ Xbeta
-w=matrix(w,ncol=(p-1),byrow=T)
-maxw=apply(w,1,max)
-y=apply(w,1,indmax)
-y=ifelse(maxw < 0,p,y)
-
-return(list(y=y,X=X,beta=beta,sigma=sigma))
-}
-
diff --git a/R/simmvp.R b/R/simmvp.R
deleted file mode 100755
index ef2915a..0000000
--- a/R/simmvp.R
+++ /dev/null
@@ -1,18 +0,0 @@
-simmvp=
-function(X,p,n,beta,sigma)
-{
-#
-# simulate from MVP: w_i=X_i*beta + e_i e ~ MVN(0,Sigma)
-# y_ij = 1 if w_ij > 0
-#
-# Sigma is p x p
-# X is np x k
-# beta is k x 1
-
-#
-Xbeta=X%*%beta
-w=as.vector(crossprod(chol(sigma),matrix(rnorm(p*n),ncol=n)))+ Xbeta
-y=ifelse(w<0,0,1)
-
-return(list(y=y,X=X,beta=beta,sigma=sigma))
-}
diff --git a/R/summary.bayesm.mat.R b/R/summary.bayesm.mat.R
new file mode 100755
index 0000000..c09b1a4
--- /dev/null
+++ b/R/summary.bayesm.mat.R
@@ -0,0 +1,48 @@
+summary.bayesm.mat=function(object,names,burnin=trunc(.1*nrow(X)),tvalues,QUANTILES=TRUE,TRAILER=TRUE,...){
+#
+# S3 method to compute and print posterior summaries for a matrix of draws
+# P. Rossi 2/07
+#
+ X=object
+ if(mode(X) == "list") stop("list entered \n Possible Fixup: extract from list \n")
+ if(mode(X) !="numeric") stop("Requires numeric argument \n")
+ if(is.null(attributes(X)$dim)) X=as.matrix(X)
+ nx=ncol(X)
+ if(missing(names)) names=as.character(1:nx)
+ if(nrow(X) < 100) {cat("fewer than 100 draws submitted \n"); return(invisible())}
+ X=X[(burnin+1):nrow(X),,drop=FALSE]
+ mat=matrix(apply(X,2,mean),nrow=1)
+ mat=rbind(mat,sqrt(matrix(apply(X,2,var),nrow=1)))
+ num_se=double(nx); rel_eff=double(nx); eff_s_size=double(nx)
+ for(i in 1:nx)
+ {out=numEff(X[,i])
+ if(is.nan(out$stderr))
+ {num_se[i]=-9999; rel_eff[i]=-9999; eff_s_size[i]=-9999}
+ else
+ {num_se[i]=out$stderr; rel_eff[i]=out$f; eff_s_size[i]=nrow(X)/ceiling(out$f)}
+ }
+ mat=rbind(mat,num_se,rel_eff,eff_s_size)
+ colnames(mat)=names
+ rownames(mat)[1]="mean"
+ rownames(mat)[2]="std dev"
+ rownames(mat)[3]="num se"
+ rownames(mat)[4]="rel eff"
+ rownames(mat)[5]="s size"
+ if(!missing(tvalues))
+ {if(mode(tvalues)!="numeric") stop("true values arguments must be numeric \n")
+ if(length(tvalues) != nx) stop("true values argument is wrong length \n")
+ mat=rbind(tvalues,mat) }
+ cat("Summary of Posterior Marginal Distributions ")
+ cat("\nMoments \n")
+ print(t(mat),digits=2)
+ if(QUANTILES){
+ qmat=apply(X,2,quantile,probs=c(.025,.05,.5,.95,.975))
+ colnames(qmat)=names
+ if(!missing(tvalues))
+ { qmat=rbind(tvalues,qmat)}
+ cat("\nQuantiles \n")
+ print(t(qmat),digits=2)}
+ if(TRAILER) cat(paste(" based on ",nrow(X)," valid draws (burn-in=",burnin,") \n",sep=""))
+ invisible()
+}
+
diff --git a/R/summary.bayesm.nmix.R b/R/summary.bayesm.nmix.R
new file mode 100755
index 0000000..bbb2a6d
--- /dev/null
+++ b/R/summary.bayesm.nmix.R
@@ -0,0 +1,35 @@
+summary.bayesm.nmix=function(object,names,burnin=trunc(.1*nrow(probdraw)),...){
+ nmixlist=object
+ if(mode(nmixlist) != "list") stop(" Argument must be a list \n")
+ probdraw=nmixlist[[1]]; compdraw=nmixlist[[3]]
+ if(!is.matrix(probdraw)) stop(" First Element of List (probdraw) must be a matrix \n")
+ if(mode(compdraw) != "list") stop(" Third Element of List (compdraw) must be a list \n")
+ ncomp=length(compdraw[[1]])
+ if(ncol(probdraw) != ncomp) stop(" Dim of First Element of List not compatible with Dim of Second
+ \n")
+#
+# function to summarize draws of normal mixture components
+#
+R=nrow(probdraw)
+if(R < 100) {cat("fewer than 100 draws submitted \n"); return(invisible())}
+datad=length(compdraw[[1]][[1]]$mu)
+mumat=matrix(0,nrow=R,ncol=datad)
+sigmat=matrix(0,nrow=R,ncol=(datad*datad))
+if(missing(names)) names=as.character(1:datad)
+for(i in (burnin+1):R){
+ if(i%%500 ==0) cat("processing draw ",i,"\n",sep="");fsh()
+ out=momMix(probdraw[i,,drop=FALSE],compdraw[i])
+ mumat[i,]=out$mu
+ sigmat[i,]=out$sigma
+ }
+cat("\nNormal Mixture Moments\n Mean\n")
+attributes(mumat)$class="bayesm.mat"
+attributes(sigmat)$class="bayesm.var"
+summary(mumat,names,QUANTILES=FALSE,TRAILER=FALSE)
+cat(" \n")
+summary(sigmat)
+cat("note: 1st and 2nd Moments for a Normal Mixture \n")
+cat(" may not be interpretable, consider plots\n")
+invisible()
+}
+
diff --git a/R/summary.bayesm.var.R b/R/summary.bayesm.var.R
new file mode 100755
index 0000000..b17f4a6
--- /dev/null
+++ b/R/summary.bayesm.var.R
@@ -0,0 +1,37 @@
+summary.bayesm.var=function(object,names,burnin=trunc(.1*nrow(Vard)),tvalues,QUANTILES=FALSE,...){
+#
+# S3 method to summarize draws of var-cov matrix (stored as a vector)
+# Vard is R x d**2 array of draws
+# P. Rossi 2/07
+#
+ Vard=object
+ if(mode(Vard) == "list") stop("list entered \n Possible Fixup: extract from list \n")
+ if(!is.matrix(Vard)) stop("Requires matrix argument \n")
+ if(trunc(sqrt(ncol(Vard)))!=sqrt(ncol(Vard))) stop("Argument cannot be draws from a square matrix \n")
+ if(nrow(Vard) < 100) {cat("fewer than 100 draws submitted \n"); return(invisible())}
+ d=sqrt(ncol(Vard))
+ corrd=t(apply(Vard[(burnin+1):nrow(Vard),],1,nmat))
+ pmeancorr=apply(corrd,2,mean)
+ dim(pmeancorr)=c(d,d)
+ indexdiag=(0:(d-1))*d+1:d
+ var=Vard[(burnin+1):nrow(Vard),indexdiag]
+ sdd=sqrt(var)
+ pmeansd=apply(sdd,2,mean)
+ mat=cbind(pmeansd,pmeancorr)
+ if(missing(names)) names=as.character(1:d)
+
+ cat("Posterior Means of Std Deviations and Correlation Matrix \n")
+ rownames(mat)=names
+ colnames(mat)=c("Std Dev",names)
+ print(mat,digits=2)
+
+ cat("\nUpper Triangle of Var-Cov Matrix \n")
+ ind=as.vector(upper.tri(matrix(0,ncol=d,nrow=d),diag=TRUE))
+ labels=cbind(rep(c(1:d),d),rep(c(1:d),each=d))
+ labels=labels[ind,]
+ plabels=paste(labels[,1],labels[,2],sep=",")
+ uppertri=as.matrix(Vard[,ind])
+ attributes(uppertri)$class="bayesm.mat"
+ summary(uppertri,names=plabels,tvalues=tvalues,QUANTILES=QUANTILES)
+ invisible()
+}
diff --git a/data/orangeJuice.rda b/data/orangeJuice.rda
new file mode 100755
index 0000000..6c0c69e
Binary files /dev/null and b/data/orangeJuice.rda differ
diff --git a/data/tuna.rda b/data/tuna.rda
new file mode 100755
index 0000000..33eed3c
Binary files /dev/null and b/data/tuna.rda differ
diff --git a/inst/doc/bayesm-manual.pdf b/inst/doc/bayesm-manual.pdf
index 22a0791..f19cf6c 100755
Binary files a/inst/doc/bayesm-manual.pdf and b/inst/doc/bayesm-manual.pdf differ
diff --git a/man/Scotch.Rd b/man/Scotch.Rd
index 870680b..2ae15c2 100755
--- a/man/Scotch.Rd
+++ b/man/Scotch.Rd
@@ -50,7 +50,7 @@ print(mat)
##
## use Scotch data to run Multivariate Probit Model
##
-if(nchar(Sys.getenv("LONG_TEST")) != 0){
+if(0){
##
y=as.matrix(Scotch)
@@ -72,13 +72,13 @@ out=rmvpGibbs(Data=Data,Mcmc=Mcmc)
ind=(0:(p-1))*p + (1:p)
cat(" Betadraws ",fill=TRUE)
mat=apply(out$betadraw/sqrt(out$sigmadraw[,ind]),2,quantile,probs=c(.01,.05,.5,.95,.99))
-print(mat)
+attributes(mat)$class="bayesm.mat"
+summary(mat)
rdraw=matrix(double((R)*p*p),ncol=p*p)
rdraw=t(apply(out$sigmadraw,1,nmat))
+attributes(rdraw)$class="bayesm.var"
cat(" Draws of Correlation Matrix ",fill=TRUE)
-mat=apply(rdraw,2,quantile,probs=c(.01,.05,.5,.95,.99))
-## correlation matrix too large to print -- summarize
-quantile(round(mat,digits=2))
+summary(rdraw)
}
diff --git a/man/bank.Rd b/man/bank.Rd
index f83d280..0378be5 100755
--- a/man/bank.Rd
+++ b/man/bank.Rd
@@ -69,13 +69,13 @@ print(mat)
## example of processing for use with rhierBinLogit
##
-if(nchar(Sys.getenv("LONG_TEST")) != 0)
+if(0)
{
choiceAtt=bank$choiceAtt
Z=bank$demo
## center demo data so that mean of random-effects
-## distribution can be interpretted as the average respondents
+## distribution can be interpreted as the average respondent
Z[,1]=rep(1,nrow(Z))
Z[,2]=Z[,2]-mean(Z[,2])
@@ -101,13 +101,27 @@ Mcmc=list(R=10000,sbeta=0.2,keep=20)
set.seed(66)
out=rhierBinLogit(Data=Data,Mcmc=Mcmc)
-cat(" Deltadraws ",fill=TRUE)
-mat=apply(out$Deltadraw,2,quantile,probs=c(.01,.05,.5,.95,.99))
-print(mat)
-cat(" Vbetadraws ",fill=TRUE)
-mat=apply(out$Vbetadraw,2,quantile,probs=c(.01,.05,.5,.95,.99))
-print(mat)
-}
+begin=5000/20
+end=10000/20
+
+summary(out$Deltadraw,burnin=begin)
+summary(out$Vbetadraw,burnin=begin)
+
+if(0){
+## plotting examples
+## plot grand means of random effects distribution (first row of Delta)
+index=4*c(0:13)+1
+matplot(out$Deltadraw[,index],type="l",xlab="Iterations/20",ylab="",
+main="Average Respondent Part-Worths")
+
+## plot hierarchical coefs
+plot(out$betadraw)
+
+## plot log-likelihood
+plot(out$llike,type="l",xlab="Iterations/20",ylab="",main="Log Likelihood")
+
+}
+}
}
\keyword{datasets}
diff --git a/man/cheese.Rd b/man/cheese.Rd
index a9e81fb..7e8b2b7 100755
--- a/man/cheese.Rd
+++ b/man/cheese.Rd
@@ -34,7 +34,7 @@ print(mat)
##
## example of processing for use with rhierLinearModel
##
-if(nchar(Sys.getenv("LONG_TEST")) != 0)
+if(0)
{
retailer=levels(cheese$RETAILER)
@@ -64,32 +64,18 @@ R=2000
Data=list(regdata=regdata,Z=Z)
Mcmc=list(R=R,keep=1)
-betamean=array(double(nreg*nvar),dim=c(nreg,nvar))
-burnin=100
-
set.seed(66)
out=rhierLinearModel(Data=Data,Mcmc=Mcmc)
-for (k in 1:nvar) { betamean[,k]=apply(out$betadraw[,k,burnin:R],1,mean)}
-print(betamean)
-cat(" Deltadraws ",fill=TRUE)
-mat=apply(out$Deltadraw,2,quantile,probs=c(.01,.05,.5,.95,.99))
-print(mat)
-cat(" Vbetadraws ",fill=TRUE)
-mat=apply(out$Vbetadraw,2,quantile,probs=c(.01,.05,.5,.95,.99))
-print(mat)
+cat("Summary of Delta Draws",fill=TRUE)
+summary(out$Deltadraw)
+cat("Summary of Vbeta Draws",fill=TRUE)
+summary(out$Vbetadraw)
if(0){
-coefn=c("Intercept","Display","LnPrice")
-colors=c("blue","green","red","yellow")
-par(mfrow=c(nvar,1),mar=c(5.1,15,4.1,13))
-for (n in 1:nvar)
-{
- plot(range(betamean[,n]),range(betamean[,n]),
- type="n",main=coefn[n],xlab="ls coef",ylab="post mean")
- points(lscoef[,n],betamean[,n],pch=17,col="blue",cex=1.2)
- abline(c(0,1))
-}
+#
+# plot hier coefs
+plot(out$betadraw)
}
}
diff --git a/man/customerSat.Rd b/man/customerSat.Rd
index 7f7f899..d680195 100755
--- a/man/customerSat.Rd
+++ b/man/customerSat.Rd
@@ -1,7 +1,7 @@
\name{customerSat}
\alias{customerSat}
\docType{data}
-\title{ Customer Satifaction Data}
+\title{ Customer Satisfaction Data}
\description{
Responses to a satisfaction survey for a Yellow Pages advertising product.
All responses are on a 10 point scale from 1 to 10 (10 is "Excellent"
diff --git a/man/detailing.Rd b/man/detailing.Rd
index f76b99d..eb8840c 100755
--- a/man/detailing.Rd
+++ b/man/detailing.Rd
@@ -46,9 +46,8 @@ print(mat)
##
## example of processing for use with rhierNegbinRw
##
-if(nchar(Sys.getenv("LONG_TEST")) != 0)
+if(0)
{
-
data(detailing)
counts = detailing$counts
Z = detailing$demo
@@ -100,14 +99,19 @@ ndraws = length(out$alphadraw)
for (i in 1:nreg) { Mbeta[i,] = rowSums(out$Betadraw[i, , ])/ndraws }
cat(" Deltadraws ",fill=TRUE)
-mat=apply(out$Deltadraw,2,quantile,probs=c(.01,.05,.5,.95,.99))
-print(mat)
+summary(out$Deltadraw)
cat(" Vbetadraws ",fill=TRUE)
-mat=apply(out$Vbetadraw,2,quantile,probs=c(.01,.05,.5,.95,.99))
-print(mat)
+summary(out$Vbetadraw)
cat(" alphadraws ",fill=TRUE)
-mat=apply(matrix(out$alphadraw),2,quantile,probs=c(.01,.05,.5,.95,.99))
-print(mat)
+summary(out$alphadraw)
+
+if(0){
+## plotting examples
+plot(out$betadraw)
+plot(out$alphadraw)
+plot(out$Deltadraw)
+}
}
+
}
\keyword{datasets}
diff --git a/man/init.rmultiregfp.Rd b/man/init.rmultiregfp.Rd
deleted file mode 100755
index 4a33491..0000000
--- a/man/init.rmultiregfp.Rd
+++ /dev/null
@@ -1,46 +0,0 @@
-\name{init.rmultiregfp}
-\alias{init.rmultiregfp}
-
-\title{ Initialize Variables for Multivariate Regression Draw }
-\description{
- \code{init.rmultiregfp} initializes variables which can be pre-computed
- for draws from the posterior of a multivariate regression model.
- \code{init.rmultiregfp} should be called prior to use of \code{rmultiregfp}
-}
-\usage{
-init.rmultiregfp(X, A, Bbar, nu, V)
-}
-\arguments{
- \item{X}{ Design matrix }
- \item{A}{ Prior Precision matrix (m x m) }
- \item{Bbar}{ Prior mean matrix (m x k) }
- \item{nu}{ degrees of freedom parmeter for Sigma prior }
- \item{V}{ location parameter for Sigma prior }
-}
-\details{
- model: \eqn{Y = XB + U}. \eqn{u_i} \eqn{\sim}{~} \eqn{N(0,\Sigma)}. \eqn{u_i} is the ith row of U. \eqn{Y} is n x m. \eqn{X} is n x k. \eqn{B} is k x m.
-
- priors: \eqn{vec(B)} \eqn{\sim}{~} \eqn{N(vec(Bbar,\Sigma (x) A^{-1})} \cr
- \eqn{\Sigma} \eqn{\sim}{~} \eqn{IW(nu,V)}.
-}
-\value{
- A list containing \ldots
- \item{IR }{Inverse of Cholesky Root of \eqn{(X'X + A)}}
- \item{RA }{ Cholesky root of A }
- \item{RABbar}{RA \%*\% Bbar}
- \item{nu}{d.f. parm for IWishart prior}
- \item{V}{location matrix for IWishart prior}
-}
-\references{ For further discussion, see \emph{Bayesian Statistics and Marketing}
- by Rossi, Allenby and McCulloch, Chapter 2. \cr
- \url{http://faculty.chicagogsb.edu/peter.rossi/research/bsm.html}
-}
-
-\author{ Peter Rossi, Graduate School of Business, University of Chicago,
- \email{Peter.Rossi at ChicagoGsb.edu}.
-}
-
-\seealso{ \code{\link{rmultiregfp}}}
-
-\keyword{ utilities }
-\keyword{ regression }
diff --git a/man/margarine.Rd b/man/margarine.Rd
index c69cf1d..42066c4 100755
--- a/man/margarine.Rd
+++ b/man/margarine.Rd
@@ -68,7 +68,7 @@ print(mat)
##
## example of processing for use with rhierMnlRwMixture
##
-if(nchar(Sys.getenv("LONG_TEST")) != 0)
+if(0)
{
select= c(1:5,7) ## select brands
chPr=as.matrix(margarine$choicePrice)
@@ -119,18 +119,14 @@ R=20000
mcmc1=list(keep=keep,R=R)
out=rhierMnlRwMixture(Data=list(p=p,lgtdata=lgtdata,Z=Z),Prior=list(ncomp=1),Mcmc=mcmc1)
-begin=100/keep; end=R/keep
-cat(" Posterior Mean of Delta ",fill=TRUE)
-mat=apply(out$Deltadraw[begin:end,],2,mean)
-print(matrix(mat,ncol=6))
-pmom=momMix(out$probdraw[begin:end,],out$compdraw[begin:end])
-cat(" posterior expectation of mu",fill=TRUE)
-print(pmom$mu)
-cat(" posterior expectation of sd",fill=TRUE)
-print(pmom$sd)
-cat(" posterior expectation of correlations",fill=TRUE)
-print(pmom$corr)
+summary(out$Deltadraw)
+summary(out$nmix)
+if(0){
+## plotting examples
+plot(out$nmix)
+plot(out$Deltadraw)}
}
+
}
\keyword{datasets}
diff --git a/man/orangeJuice.Rd b/man/orangeJuice.Rd
new file mode 100755
index 0000000..c208105
--- /dev/null
+++ b/man/orangeJuice.Rd
@@ -0,0 +1,165 @@
+\name{orangeJuice}
+\alias{orangeJuice}
+\docType{data}
+\title{Store-level Panel Data on Orange Juice Sales}
+\description{
+ yx, weekly sales of refrigerated orange juice at 83 stores. \cr
+ storedemo, contains demographic information on those stores. \cr
+}
+\usage{data(orangeJuice)}
+\format{
+ This R object is a list of two data frames, list(yx,storedemo).\cr
+
+ List of 2 \cr
+ \$ yx :'data.frame': 106139 obs. of 19 variables:\cr
+ \ldots \$ store : int [1:106139] 2 2 2 2 2 2 2 2 2 2 \cr
+ \ldots \$ brand : int [1:106139] 1 1 1 1 1 1 1 1 1 1 \cr
+ \ldots \$ week : int [1:106139] 40 46 47 48 50 51 52 53 54 57 \cr
+ \ldots \$ logmove : num [1:106139] 9.02 8.72 8.25 8.99 9.09 \cr
+ \ldots \$ constant: int [1:106139] 1 1 1 1 1 1 1 1 1 1 \cr
+ \ldots \$ price1 : num [1:106139] 0.0605 0.0605 0.0605 0.0605 0.0605 \cr
+ \ldots \$ price2 : num [1:106139] 0.0605 0.0603 0.0603 0.0603 0.0603 \cr
+ \ldots \$ price3 : num [1:106139] 0.0420 0.0452 0.0452 0.0498 0.0436 \cr
+ \ldots \$ price4 : num [1:106139] 0.0295 0.0467 0.0467 0.0373 0.0311 \cr
+ \ldots \$ price5 : num [1:106139] 0.0495 0.0495 0.0373 0.0495 0.0495 \cr
+ \ldots \$ price6 : num [1:106139] 0.0530 0.0478 0.0530 0.0530 0.0530 \cr
+ \ldots \$ price7 : num [1:106139] 0.0389 0.0458 0.0458 0.0458 0.0466 \cr
+ \ldots \$ price8 : num [1:106139] 0.0414 0.0280 0.0414 0.0414 0.0414 \cr
+ \ldots \$ price9 : num [1:106139] 0.0289 0.0430 0.0481 0.0423 0.0423 \cr
+ \ldots \$ price10 : num [1:106139] 0.0248 0.0420 0.0327 0.0327 0.0327 \cr
+ \ldots \$ price11 : num [1:106139] 0.0390 0.0390 0.0390 0.0390 0.0382 \cr
+ \ldots \$ deal : int [1:106139] 1 0 0 0 0 0 1 1 1 1 \cr
+ \ldots \$ feat : num [1:106139] 0 0 0 0 0 0 0 0 0 0 \cr
+ \ldots \$ profit : num [1:106139] 38.0 30.1 30.0 29.9 29.9 \cr
+
+ 1 Tropicana Premium 64 oz; 2 Tropicana Premium 96 oz; 3 Florida's Natural 64 oz; \cr
+ 4 Tropicana 64 oz; 5 Minute Maid 64 oz; 6 Minute Maid 96 oz; \cr
+ 7 Citrus Hill 64 oz; 8 Tree Fresh 64 oz; 9 Florida Gold 64 oz; \cr
+ 10 Dominicks 64 oz; 11 Dominicks 128 oz. \cr
+
+ \$ storedemo:'data.frame': 83 obs. of 12 variables:\cr
+ \ldots \$ STORE : int [1:83] 2 5 8 9 12 14 18 21 28 32 \cr
+ \ldots \$ AGE60 : num [1:83] 0.233 0.117 0.252 0.269 0.178 \cr
+ \ldots \$ EDUC : num [1:83] 0.2489 0.3212 0.0952 0.2222 0.2534 \cr
+ \ldots \$ ETHNIC : num [1:83] 0.1143 0.0539 0.0352 0.0326 0.3807 \cr
+ \ldots \$ INCOME : num [1:83] 10.6 10.9 10.6 10.8 10.0 \cr
+ \ldots \$ HHLARGE : num [1:83] 0.1040 0.1031 0.1317 0.0968 0.0572 \cr
+ \ldots \$ WORKWOM : num [1:83] 0.304 0.411 0.283 0.359 0.391 \cr
+ \ldots \$ HVAL150 : num [1:83] 0.4639 0.5359 0.0542 0.5057 0.3866 \cr
+ \ldots \$ SSTRDIST: num [1:83] 2.11 3.80 2.64 1.10 9.20 \cr
+ \ldots \$ SSTRVOL : num [1:83] 1.143 0.682 1.500 0.667 1.111 \cr
+ \ldots \$ CPDIST5 : num [1:83] 1.93 1.60 2.91 1.82 0.84 \cr
+ \ldots \$ CPWVOL5 : num [1:83] 0.377 0.736 0.641 0.441 0.106 \cr
+}
+\details{
+ \describe{
+ \item{\code{store}}{store number}
+ \item{\code{brand}}{brand indicator}
+ \item{\code{week}}{week number}
+ \item{\code{logmove}}{log of the number of units sold}
+ \item{\code{constant}}{a vector of 1}
+ \item{\code{price1}}{price of brand 1}
+ \item{\code{deal}}{in-store coupon activity}
+ \item{\code{feature}}{feature advertisement}
+ \item{\code{STORE}}{store number}
+ \item{\code{AGE60}}{percentage of the population that is aged 60 or older}
+ \item{\code{EDUC}}{percentage of the population that has a college degree}
+ \item{\code{ETHNIC}}{percent of the population that is black or Hispanic}
+ \item{\code{INCOME}}{median income}
+ \item{\code{HHLARGE}}{percentage of households with 5 or more persons}
+ \item{\code{WORKWOM}}{percentage of women with full-time jobs}
+ \item{\code{HVAL150}}{percentage of households worth more than \$150,000}
+ \item{\code{SSTRDIST}}{distance to the nearest warehouse store}
+ \item{\code{SSTRVOL}}{ratio of sales of this store to the nearest warehouse store}
+ \item{\code{CPDIST5}}{average distance in miles to the nearest 5 supermarkets}
+ \item{\code{CPWVOL5}}{ratio of sales of this store to the average of the nearest five stores}
+ }
+}
+\source{
+ Alan L. Montgomery (1997), "Creating Micro-Marketing Pricing Strategies Using Supermarket Scanner Data,"
+ \emph{Marketing Science} 16(4) 315-337.
+}
+\references{
+ Chapter 5, \emph{Bayesian Statistics and Marketing} by Rossi et al.\cr
+ \url{http://faculty.chicagogsb.edu/peter.rossi/research/bsm.html}
+}
+\examples{
+
+## Example
+## load data
+data(orangeJuice)
+
+## print some quantiles of yx data
+cat("Quantiles of the Variables in yx data",fill=TRUE)
+mat=apply(as.matrix(orangeJuice$yx),2,quantile)
+print(mat)
+
+## print some quantiles of storedemo data
+cat("Quantiles of the Variables in storedemo data",fill=TRUE)
+mat=apply(as.matrix(orangeJuice$storedemo),2,quantile)
+print(mat)
+
+
+## Example 2 processing for use with rhierLinearModel
+##
+##
+if(0)
+{
+
+## select brand 1 for analysis
+brand1=orangeJuice$yx[(orangeJuice$yx$brand==1),]
+
+store = sort(unique(brand1$store))
+nreg = length(store)
+nvar=14
+
+regdata=NULL
+for (reg in 1:nreg) {
+ y=brand1$logmove[brand1$store==store[reg]]
+ iota=c(rep(1,length(y)))
+ X=cbind(iota,log(brand1$price1[brand1$store==store[reg]]),
+ log(brand1$price2[brand1$store==store[reg]]),
+ log(brand1$price3[brand1$store==store[reg]]),
+ log(brand1$price4[brand1$store==store[reg]]),
+ log(brand1$price5[brand1$store==store[reg]]),
+ log(brand1$price6[brand1$store==store[reg]]),
+ log(brand1$price7[brand1$store==store[reg]]),
+ log(brand1$price8[brand1$store==store[reg]]),
+ log(brand1$price9[brand1$store==store[reg]]),
+ log(brand1$price10[brand1$store==store[reg]]),
+ log(brand1$price11[brand1$store==store[reg]]),
+ brand1$deal[brand1$store==store[reg]],
+ brand1$feat[brand1$store==store[reg]])
+ regdata[[reg]]=list(y=y,X=X)
+ }
+
+## storedemo is standardized to zero mean.
+
+Z=as.matrix(orangeJuice$storedemo[,2:12])
+dmean=apply(Z,2,mean)
+for (s in 1:nreg){
+ Z[s,]=Z[s,]-dmean
+}
+iotaz=c(rep(1,nrow(Z)))
+Z=cbind(iotaz,Z)
+nz=ncol(Z)
+
+
+Data=list(regdata=regdata,Z=Z)
+Mcmc=list(R=R,keep=1)
+
+out=rhierLinearModel(Data=Data,Mcmc=Mcmc)
+
+summary(out$Deltadraw)
+summary(out$Vbetadraw)
+
+if(0){
+## plotting examples
+plot(out$betadraw)
+}
+}
+
+}
+\keyword{datasets}
+
+
diff --git a/man/plot.bayesm.hcoef.Rd b/man/plot.bayesm.hcoef.Rd
new file mode 100755
index 0000000..44346ec
--- /dev/null
+++ b/man/plot.bayesm.hcoef.Rd
@@ -0,0 +1,42 @@
+\name{plot.bayesm.hcoef}
+\alias{plot.bayesm.hcoef}
+\concept{MCMC}
+\concept{S3 method}
+\concept{plot}
+\concept{hierarchical model}
+\title{Plot Method for Hierarchical Model Coefs }
+\description{
+ \code{plot.bayesm.hcoef} is an S3 method to plot 3 dim arrays of hierarchical coefficients.
+ Arrays are of class bayesm.hcoef with dimensions: cross-sectional unit x coef x MCMC draw.
+}
+\usage{
+\method{plot}{bayesm.hcoef}(x,burnin,...)
+}
+\arguments{
+ \item{x}{ An object of S3 class, bayesm.hcoef }
+ \item{burnin}{ no draws to burnin, def: .1*R }
+ \item{...}{ standard graphics parameters }
+}
+\details{
+ Typically, \code{plot.bayesm.hcoef} will be invoked by a call to the generic plot function as in
+ \code{plot(object)} where object is of class bayesm.hcoef. All of the \code{bayesm} hierarchical routines
+ return draws of hierarchical coefficients in this class (see example below). One can also simply invoke
+ \code{plot.bayesm.hcoef} on any valid 3-dim array as in \code{plot.bayesm.hcoef(betadraws)}
+ \cr
+ \cr
+ \code{plot.bayesm.hcoef} is also exported for use as a standard function, as in
+ \code{plot.bayesm.hcoef(array)}.
+}
+\author{ Peter Rossi, Graduate School of Business, University of Chicago,
+ \email{Peter.Rossi at ChicagoGsb.edu}.
+}
+\seealso{ \code{\link{rhierMnlRwMixture}},\code{\link{rhierLinearModel}},
+ \code{\link{rhierLinearMixture}},\code{\link{rhierNegbinRw}} }
+\examples{
+##
+## not run
+# out=rhierLinearModel(Data,Prior,Mcmc)
+# plot(out$betadraws)
+#
+}
+\keyword{ hplot }
diff --git a/man/plot.bayesm.mat.Rd b/man/plot.bayesm.mat.Rd
new file mode 100755
index 0000000..259f490
--- /dev/null
+++ b/man/plot.bayesm.mat.Rd
@@ -0,0 +1,51 @@
+\name{plot.bayesm.mat}
+\alias{plot.bayesm.mat}
+\concept{MCMC}
+\concept{S3 method}
+\concept{plot}
+\title{Plot Method for Arrays of MCMC Draws}
+\description{
+ \code{plot.bayesm.mat} is an S3 method to plot arrays of MCMC draws. The columns in the
+ array correspond to parameters and the rows to MCMC draws.
+}
+\usage{
+\method{plot}{bayesm.mat}(x,names,burnin,tvalues,TRACEPLOT,DEN,INT, ...)
+}
+\arguments{
+ \item{x}{ An object of either S3 class, bayesm.mat, or S3 class, mcmc }
+ \item{names}{optional character vector of names for coefficients}
+ \item{burnin}{number of draws to discard for burn-in, def: .1*nrow(X)}
+ \item{tvalues}{vector of true values}
+ \item{TRACEPLOT}{ logical, TRUE provide sequence plots of draws and acfs, def: TRUE }
+ \item{DEN}{ logical, TRUE use density scale on histograms, def: TRUE }
+ \item{INT}{ logical, TRUE put various intervals and points on graph, def: TRUE }
+ \item{...}{ standard graphics parameters }
+}
+\details{
+ Typically, \code{plot.bayesm.mat} will be invoked by a call to the generic plot function as in
+ \code{plot(object)} where object is of class bayesm.mat. All of the \code{bayesm} MCMC routines
+ return draws in this class (see example below). One can also simply invoke
+ \code{plot.bayesm.mat} on any valid 2-dim array as in \code{plot.bayesm.mat(betadraws)}. \cr
+ \cr
+ \code{plot.bayesm.mat} paints (by default) on the histogram: \cr
+ \cr
+ green "[]" delimiting 95\% Bayesian Credibility Interval \cr
+ yellow "()" showing +/- 2 numerical standard errors \cr
+ red "|" showing posterior mean
+ \cr
+ \cr
+ \code{plot.bayesm.mat} is also exported for use as a standard function, as in
+ \code{plot.bayesm.mat(matrix)}
+}
+\author{ Peter Rossi, Graduate School of Business, University of Chicago,
+ \email{Peter.Rossi at ChicagoGsb.edu}.
+}
+\examples{
+##
+## not run
+# out=runiregGibbs(Data,Prior,Mcmc)
+# plot(out$betadraw)
+#
+}
+\keyword{ hplot }
+
diff --git a/man/plot.bayesm.nmix.Rd b/man/plot.bayesm.nmix.Rd
new file mode 100755
index 0000000..489377d
--- /dev/null
+++ b/man/plot.bayesm.nmix.Rd
@@ -0,0 +1,46 @@
+\name{plot.bayesm.nmix}
+\alias{plot.bayesm.nmix}
+\concept{MCMC}
+\concept{S3 method}
+\concept{plot}
+\title{Plot Method for MCMC Draws of Normal Mixtures}
+\description{
+ \code{plot.bayesm.nmix} is an S3 method to plot aspects of the fitted density from a list
+ of MCMC draws of normal mixture components. Plots of marginal univariate and bivariate densities
+ are produced.
+}
+\usage{
+\method{plot}{bayesm.nmix}(x,names,burnin,Grid,bi.sel,nstd, ...)
+}
+\arguments{
+ \item{x}{ An object of either S3 class, bayesm.nmix, or S3 class, mcmc }
+ \item{names}{optional character vector of names for each of the dimensions}
+ \item{burnin}{number of draws to discard for burn-in, def: .1*nrow(X)}
+ \item{Grid}{matrix of grid points for densities, def: mean +/- nstd std deviations}
+ \item{bi.sel}{list of vectors, each giving pairs for bivariate distributions, def: list(c(1,2))}
+ \item{nstd}{number of standard deviations for default Grid, def: 2}
+ \item{...}{ standard graphics parameters }
+}
+\details{
+ Typically, \code{plot.bayesm.nmix} will be invoked by a call to the generic plot function as in
+ \code{plot(object)} where object is of class bayesm.nmix. These objects are lists of three components.
+ The first component is an array of draws of mixture component probabilties. The second component
+ is not used. The third is a lists of lists of lists with draws of each of the normal components.
+ \cr
+ \cr
+ \code{plot.bayesm.nmix} can also be used as a standard function, as in \code{plot.bayesm.nmix(list)}.
+}
+\author{ Peter Rossi, Graduate School of Business, University of Chicago,
+ \email{Peter.Rossi at ChicagoGsb.edu}.
+}
+\seealso{ \code{\link{rnmixGibbs}}, \code{\link{rhierMnlRwMixture}}, \code{\link{rhierLinearMixture}}}
+\examples{
+##
+## not run
+# out=rnmix(Data,Prior,Mcmc)
+# plot(out,bi.sel=list(c(1,2),c(3,4),c(1,3)))
+# plot bivariate distributions for dimension 1,2; 3,4; and 1,3
+#
+
+}
+\keyword{ hplot }
diff --git a/man/rbprobitGibbs.Rd b/man/rbprobitGibbs.Rd
index fe3144a..4e7b02d 100755
--- a/man/rbprobitGibbs.Rd
+++ b/man/rbprobitGibbs.Rd
@@ -70,13 +70,17 @@ beta=c(0,1,-1)
nvar=ncol(X)
simout=simbprobit(X,beta)
-Data=list(X=simout$X,y=simout$y)
-Mcmc=list(R=R,keep=1)
+Data1=list(X=simout$X,y=simout$y)
+Mcmc1=list(R=R,keep=1)
-out=rbprobitGibbs(Data=Data,Mcmc=Mcmc)
+out=rbprobitGibbs(Data=Data1,Mcmc=Mcmc1)
+
+summary(out$betadraw,tvalues=beta)
+
+if(0){
+## plotting example
+plot(out$betadraw,tvalues=beta)
+}
-cat(" Betadraws ",fill=TRUE)
-mat=apply(out$betadraw,2,quantile,probs=c(.01,.05,.5,.95,.99))
-mat=rbind(beta,mat); rownames(mat)[1]="beta"; print(mat)
}
\keyword{ models }
diff --git a/man/rhierBinLogit.Rd b/man/rhierBinLogit.Rd
index 0dc6c44..e3afbea 100755
--- a/man/rhierBinLogit.Rd
+++ b/man/rhierBinLogit.Rd
@@ -101,27 +101,18 @@ for (i in 1:nlgt)
lgtdata[[i]]=list(y=y,X=X,beta=beta)
}
-Data=list(Dat=lgtdata,Demo=Z)
out=rhierBinLogit(Data=list(lgtdata=lgtdata,Z=Z),Mcmc=list(R=R))
-cat(" Deltadraws ",fill=TRUE)
-mat=apply(out$Deltadraw,2,quantile,probs=c(.01,.05,.5,.95,.99))
-mat=rbind(as.vector(Delta),mat); rownames(mat)[1]="delta"; print(mat)
-cat(" Vbetadraws ",fill=TRUE)
-mat=apply(out$Vbetadraw,2,quantile,probs=c(.01,.05,.5,.95,.99))
-mat=rbind(as.vector(Vbeta),mat); rownames(mat)[1]="Vbeta"; print(mat)
+cat("Summary of Delta draws",fill=TRUE)
+summary(out$Deltadraw,tvalues=as.vector(Delta))
+cat("Summary of Vbeta draws",fill=TRUE)
+summary(out$Vbetadraw,tvalues=as.vector(Vbeta[upper.tri(Vbeta,diag=TRUE)]))
if(0){
-td=as.vector(Delta)
-par(mfrow=c(2,2))
-matplot(out$Deltadraw[,(1:nvar)],type="l")
-abline(h=td[1:nvar],col=(1:nvar))
-matplot(out$Deltadraw[,((nvar+1):(2*nvar))],type="l")
-abline(h=td[(nvar+1):(2*nvar)],col=(1:nvar))
-matplot(out$Vbetadraw[,c(1,7,13,19,25)],type="l")
-abline(h=1.5)
-matplot(out$Vbetadraw[,-c(1,7,13,19,25)],type="l")
-abline(h=.5)
+## plotting examples
+plot(out$Deltadraw,tvalues=as.vector(Delta))
+plot(out$betadraw)
+plot(out$Vbetadraw,tvalues=as.vector(Vbeta[upper.tri(Vbeta,diag=TRUE)]))
}
}
diff --git a/man/rhierLinearMixture.Rd b/man/rhierLinearMixture.Rd
index ffe571a..09aa9b6 100755
--- a/man/rhierLinearMixture.Rd
+++ b/man/rhierLinearMixture.Rd
@@ -64,17 +64,16 @@ rhierLinearMixture(Data, Prior, Mcmc)
\item{taudraw}{R/keep x nreg array of error variance draws}
\item{betadraw}{nreg x nvar x R/keep array of individual regression coef draws}
\item{Deltadraw}{R/keep x nz x nvar array of Deltadraws}
- \item{probdraw}{ R/keep x ncomp matrix of draws of probs of mixture components (pvec)}
- \item{compdraw}{ list of list of lists (length R/keep)}
+ \item{nmix}{list of three elements, (probdraw, NULL, compdraw)}
}
\note{
- More on \code{compdraw} component of return value list: \cr
- \itemize{
+ More on \code{probdraw} component of nmix return value list: \cr
+ this is an R/keep by ncomp array of draws of mixture component probs (pvec)\cr
+ More on \code{compdraw} component of nmix return value list: \cr
\item{compdraw[[i]]}{ the ith draw of components for mixtures}
\item{compdraw[[i]][[j]]}{ ith draw of the jth normal mixture comp}
\item{compdraw[[i]][[j]][[1]]}{ ith draw of jth normal mixture comp mean vector}
\item{compdraw[[i]][[j]][[2]]}{ ith draw of jth normal mixture cov parm (rooti) }
- }
Note: Z should \strong{not} include an intercept and should be centered for ease of interpretation.\cr
@@ -117,8 +116,6 @@ regdata=NULL # simulated data with Z
betas=matrix(double(nreg*nvar),ncol=nvar)
tind=double(nreg)
-# simulate datasets with/without Z, using same components and tau
-
for (reg in 1:nreg) {
tempout=rmixture(1,tpvec,tcomps)
betas[reg,]=Delta\%*\%Z[reg,]+as.vector(tempout$x)
@@ -131,110 +128,23 @@ regdata[[reg]]=list(y=y,X=X,beta=betas[reg,],tau=tau)
## run rhierLinearMixture
-Data=list(regdata=regdata,Z=Z)
-Prior=list(ncomp=3)
-Mcmc=list(R=R,keep=1)
-
-out1=rhierLinearMixture(Data=Data,Prior=Prior,Mcmc=Mcmc)
+Data1=list(regdata=regdata,Z=Z)
+Prior1=list(ncomp=3)
+Mcmc1=list(R=R,keep=1)
-if(R>1000) {begin=501} else {begin=1}
+out1=rhierLinearMixture(Data=Data1,Prior=Prior1,Mcmc=Mcmc1)
-apply(out1$Deltadraw[begin:R,],2,mean)
-cat(" Deltadraws ",fill=TRUE)
-mat=apply(out1$Deltadraw,2,quantile,probs=c(.01,.05,.5,.95,.99))
-mat=rbind(as.vector(Delta),mat)
-rownames(mat)[1]="delta"
-print(mat)
+cat("Summary of Delta draws",fill=TRUE)
+summary(out1$Deltadraw,tvalues=as.vector(Delta))
+cat("Summary of Normal Mixture Distribution",fill=TRUE)
+summary(out1$nmix)
if(0){
## plotting examples
-## plot histograms of draws of betas for nth regression
-
-betahist=function(n)
-{
-par(mfrow=c(1,3))
-hist(out1$betadraw[n,1,begin:R],breaks=30,main="beta1 with Z",xlab=" ",ylab=" ")
-abline(v=hhbetamean1[n,1],col="red",lwd=2)
-abline(v=regdata[[n]]$beta[1],col="blue",lwd=2)
-hist(out1$betadraw[n,2,begin:R],breaks=30,main="beta2 with Z",xlab=" ",ylab=" ")
-abline(v=hhbetamean1[n,2],col="red",lwd=2)
-abline(v=regdata[[n]]$beta[2],col="blue",lwd=2)
-hist(out1$betadraw[n,3,begin:R],breaks=30,main="beta3 with Z",xlab=" ",ylab=" ")
-abline(v=hhbetamean1[n,3],col="red",lwd=2)
-abline(v=regdata[[n]]$beta[3],col="blue",lwd=2)
-}
-
-
-betahist(10) ## plot betas for 10th regression, using regdata
-betahist(20)
-betahist(30)
-
-## plot univariate marginal density of betas
-
-grid=NULL
-for (i in 1:nvar){
- rgi=range(betas[,i])
- gr=seq(from=rgi[1],to=rgi[2],length.out=50)
- grid=cbind(grid,gr)
-}
-
-tmden=mixDen(grid,tpvec,tcomps)
-pmden=eMixMargDen(grid,as.matrix(out1$probdraw[begin:R,]),out1$compdraw[begin:R])
-
-par(mfrow=c(1,3))
-
-for (i in 1:nvar){
-plot(range(grid[,i]),c(0,1.1*max(tmden[,i],pmden[,i])),type="n",xlab="",ylab="density")
-lines(grid[,i],tmden[,i],col="blue",lwd=2)
-lines(grid[,i],pmden[,i],col="red",lwd=2)
+plot(out1$betadraw)
+plot(out1$nmix)
+plot(out1$Deltadraw)
}
-
-# plot bivariate marginal density of betas
-
-end=R
-rx1=range(betas[,1])
-rx2=range(betas[,2])
-rx3=range(betas[,3])
-x1=seq(from=rx1[1],to=rx1[2],length.out=50)
-x2=seq(from=rx2[1],to=rx2[2],length.out=50)
-x3=seq(from=rx3[1],to=rx3[2],length.out=50)
-den12=matrix(0,ncol=length(x1),nrow=length(x2))
-den23=matrix(0,ncol=length(x2),nrow=length(x3))
-den13=matrix(0,ncol=length(x1),nrow=length(x3))
-
-for(ind in as.integer(seq(from=begin,to=end,length.out=100))){
-den12=den12+mixDenBi(1,2,x1,x2,as.matrix(out1$probdraw[ind,]),out1$compdraw[[ind]])
-den23=den23+mixDenBi(2,3,x2,x3,as.matrix(out1$probdraw[ind,]),out1$compdraw[[ind]])
-den13=den13+mixDenBi(1,3,x1,x3,as.matrix(out1$probdraw[ind,]),out1$compdraw[[ind]])
-}
-
-tden12=matrix(0,ncol=length(x1),nrow=length(x2))
-tden23=matrix(0,ncol=length(x2),nrow=length(x3))
-tden13=matrix(0,ncol=length(x1),nrow=length(x3))
-tden12=mixDenBi(1,2,x1,x2,tpvec,tcomps)
-tden23=mixDenBi(2,3,x2,x3,tpvec,tcomps)
-tden13=mixDenBi(1,3,x1,x3,tpvec,tcomps)
-
-par(mfrow=c(2,3))
-image(x1,x2,tden12,col=terrain.colors(100),xlab="",ylab="")
-contour(x1,x2,tden12,add=TRUE,drawlabels=FALSE)
-title("True vars 1&2")
-image(x2,x3,tden23,col=terrain.colors(100),xlab="",ylab="")
-contour(x2,x3,tden23,add=TRUE,drawlabels=FALSE)
-title("True vars 2&3")
-image(x1,x3,tden13,col=terrain.colors(100),xlab="",ylab="")
-contour(x1,x3,tden13,add=TRUE,drawlabels=FALSE)
-title("True vars 1&3")
-image(x1,x2,den12,col=terrain.colors(100),xlab="",ylab="")
-contour(x1,x2,den12,add=TRUE,drawlabels=FALSE)
-title("Posterior vars 1&2")
-image(x2,x3,den23,col=terrain.colors(100),xlab="",ylab="")
-contour(x2,x3,den23,add=TRUE,drawlabels=FALSE)
-title("Posterior vars 2&3")
-image(x1,x3,den13,col=terrain.colors(100),xlab="",ylab="")
-contour(x1,x3,den13,add=TRUE,drawlabels=FALSE)
-title("Posterior vars 1&3")
-}
}
\keyword{ regression }
diff --git a/man/rhierLinearModel.Rd b/man/rhierLinearModel.Rd
index 6e82793..f38b1f3 100755
--- a/man/rhierLinearModel.Rd
+++ b/man/rhierLinearModel.Rd
@@ -63,6 +63,7 @@ rhierLinearModel(Data, Prior, Mcmc)
\author{ Peter Rossi, Graduate School of Business, University of Chicago,
\email{Peter.Rossi at ChicagoGsb.edu}.
}
+\seealso{ \code{\link{rhierLinearMixture}} }
\examples{
##
if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=2000} else {R=10}
@@ -80,26 +81,20 @@ regdata=NULL
for (reg in 1:nreg) { X=cbind(iota,matrix(runif(nobs*(nvar-1)),ncol=(nvar-1)))
y=X\%*\%Beta[reg,]+sqrt(tau)*rnorm(nobs); regdata[[reg]]=list(y=y,X=X) }
-nu.e=3
-ssq=NULL
-for(reg in 1:nreg) {ssq[reg]=var(regdata[[reg]]$y)}
-nu=nvar+3
-V=nu*diag(c(rep(1,nvar)))
-A=diag(c(rep(.01,nz)),ncol=nz)
-Deltabar=matrix(c(rep(0,nz*nvar)),nrow=nz)
+Data1=list(regdata=regdata,Z=Z)
+Mcmc1=list(R=R,keep=1)
+out=rhierLinearModel(Data=Data1,Mcmc=Mcmc1)
+cat("Summary of Delta draws",fill=TRUE)
+summary(out$Deltadraw,tvalues=as.vector(Delta))
+cat("Summary of Vbeta draws",fill=TRUE)
+summary(out$Vbetadraw,tvalues=as.vector(Vbeta[upper.tri(Vbeta,diag=TRUE)]))
-Data=list(regdata=regdata,Z=Z)
-Prior=list(Deltabar=Deltabar,A=A,nu.e=nu.e,ssq=ssq,nu=nu,V=V)
-Mcmc=list(R=R,keep=1)
-out=rhierLinearModel(Data=Data,Mcmc=Mcmc)
-
-cat(" Deltadraws ",fill=TRUE)
-mat=apply(out$Deltadraw,2,quantile,probs=c(.01,.05,.5,.95,.99))
-mat=rbind(as.vector(Delta),mat); rownames(mat)[1]="delta"; print(mat)
-cat(" Vbetadraws ",fill=TRUE)
-mat=apply(out$Vbetadraw,2,quantile,probs=c(.01,.05,.5,.95,.99))
-mat=rbind(as.vector(Vbeta),mat); rownames(mat)[1]="Vbeta"; print(mat)
+if(0){
+## plotting examples
+plot(out$betadraw)
+plot(out$Deltadraw)
+}
}
\keyword{ regression }
diff --git a/man/rhierMnlRwMixture.Rd b/man/rhierMnlRwMixture.Rd
index 2bfba08..ae4a6e8 100755
--- a/man/rhierMnlRwMixture.Rd
+++ b/man/rhierMnlRwMixture.Rd
@@ -19,7 +19,7 @@ rhierMnlRwMixture(Data, Prior, Mcmc)
}
\arguments{
\item{Data}{ list(p,lgtdata,Z) ( Z is optional) }
- \item{Prior}{ list(deltabar,Ad,mubar,Amu,nu,V,ncomp) (all but ncomp are optional)}
+ \item{Prior}{ list(a,deltabar,Ad,mubar,Amu,nu,V,ncomp) (all but ncomp are optional)}
\item{Mcmc}{ list(s,w,R,keep) (R required)}
}
\details{
@@ -44,6 +44,7 @@ rhierMnlRwMixture(Data, Prior, Mcmc)
\item{\code{lgtdata}}{list of lists with each cross-section unit MNL data}
\item{\code{lgtdata[[i]]$y}}{ \eqn{n_i} vector of multinomial outcomes (1,\ldots,m)}
\item{\code{lgtdata[[i]]$X}}{ \eqn{n_i} by nvar design matrix for ith unit}
+ \item{\code{a}}{vector of length ncomp of Dirichlet prior parms (def: rep(5,ncomp))}
\item{\code{deltabar}}{nz*nvar vector of prior means (def: 0)}
\item{\code{Ad}}{ prior prec matrix for vec(D) (def: .01I)}
\item{\code{mubar}}{ nvar x 1 prior mean vector for normal comp mean (def: 0)}
@@ -61,11 +62,12 @@ rhierMnlRwMixture(Data, Prior, Mcmc)
a list containing:
\item{Deltadraw}{R/keep x nz*nvar matrix of draws of Delta, first row is initial value}
\item{betadraw}{ nlgt x nvar x R/keep array of draws of betas}
- \item{probdraw}{ R/keep x ncomp matrix of draws of probs of mixture components (pvec)}
- \item{compdraw}{ list of list of lists (length R/keep)}
+ \item{nmix}{ list of 3 components, probdraw, NULL, compdraw }
\item{loglike}{ log-likelihood for each kept draw (length R/keep)}
}
\note{
+ More on \code{probdraw} component of nmix list:\cr
+ R/keep x ncomp matrix of draws of probs of mixture components (pvec) \cr
More on \code{compdraw} component of return value list: \cr
\itemize{
\item{compdraw[[i]]}{ the ith draw of components for mixtures}
@@ -74,7 +76,7 @@ rhierMnlRwMixture(Data, Prior, Mcmc)
\item{compdraw[[i]][[j]][[2]]}{ ith draw of jth normal mixture cov parm (rooti) }
}
- Note: Z does \strong{not} include an intercept and is centered for ease of interpretation.\cr
+ Note: Z should \strong{not} include an intercept and is centered for ease of interpretation.\cr
Be careful in assessing prior parameter, Amu. .01 is too small for many applications. See
Rossi et al, chapter 5 for full discussion.\cr
@@ -101,8 +103,7 @@ rhierMnlRwMixture(Data, Prior, Mcmc)
\seealso{ \code{\link{rmnlIndepMetrop}} }
\examples{
##
-if(nchar(Sys.getenv("LONG_TEST")) != 0)
-{
+if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=10000} else {R=10}
set.seed(66)
p=3 # num of choice alterns
@@ -119,15 +120,30 @@ comps[[2]]=list(mu=c(0,-1,-2)*2,rooti=diag(rep(1,3)))
comps[[3]]=list(mu=c(0,-1,-2)*4,rooti=diag(rep(1,3)))
pvec=c(.4,.2,.4)
+simmnlwX= function(n,X,beta) {
+ ## simulate from MNL model conditional on X matrix
+ k=length(beta)
+ Xbeta=X\%*\%beta
+ j=nrow(Xbeta)/n
+ Xbeta=matrix(Xbeta,byrow=TRUE,ncol=j)
+ Prob=exp(Xbeta)
+ iota=c(rep(1,j))
+ denom=Prob\%*\%iota
+ Prob=Prob/as.vector(denom)
+ y=vector("double",n)
+ ind=1:j
+ for (i in 1:n)
+ {yvec=rmultinom(1,1,Prob[i,]); y[i]=ind\%*\%yvec}
+ return(list(y=y,X=X,beta=beta,prob=Prob))
+}
+
## simulate data
simlgtdata=NULL
ni=rep(50,300)
for (i in 1:nlgt)
{ betai=Delta\%*\%Z[i,]+as.vector(rmixture(1,pvec,comps)$x)
- X=NULL
- for(j in 1:ni[i])
- { Xone=cbind(rbind(c(rep(0,p-1)),diag(p-1)),runif(p,min=-1.5,max=0))
- X=rbind(X,Xone) }
+ Xa=matrix(runif(ni[i]*p,min=-1.5,max=0),ncol=p)
+ X=createX(p,na=1,nd=NULL,Xa=Xa,Xd=NULL,base=1)
outa=simmnlwX(ni[i],X,betai)
simlgtdata[[i]]=list(y=outa$y,X=X,beta=betai)
}
@@ -142,80 +158,25 @@ for(i in 1:ncoef) hist(bmat[,i],breaks=30,col="magenta")
}
## set parms for priors and Z
-nu=ncoef+3
-V=nu*diag(rep(1,ncoef))
-Ad=.01*(diag(rep(1,nz*ncoef)))
-mubar=matrix(rep(0,ncoef),nrow=1)
-deltabar=rep(0,ncoef*nz)
-Amu=matrix(.01,ncol=1)
-a=rep(5,ncomp)
-
-R=10000
+Prior1=list(ncomp=5)
+
keep=5
-w=.1
-s=2.93/sqrt(ncoef)
+Mcmc1=list(R=R,keep=keep)
Data1=list(p=p,lgtdata=simlgtdata,Z=Z)
-Prior1=list(ncomp=ncomp,nu=nu,V=V,Amu=Amu,mubar=mubar,a=a,Ad=Ad,deltabar=deltabar)
-Mcmc1=list(s=s,w=w,R=R,keep=keep)
+
out=rhierMnlRwMixture(Data=Data1,Prior=Prior1,Mcmc=Mcmc1)
-if(R < 1000) {begin=1} else {begin=1000/keep}
-end=R/keep
-cat(" Deltadraws ",fill=TRUE)
-mat=apply(out$Deltadraw[begin:end,],2,quantile,probs=c(.01,.05,.5,.95,.99))
-mat=rbind(as.vector(Delta),mat); rownames(mat)[1]="delta"; print(mat)
-
-tmom=momMix(matrix(pvec,nrow=1),list(comps))
-pmom=momMix(out$probdraw[begin:end,],out$compdraw[begin:end])
-mat=rbind(tmom$mu,pmom$mu)
-rownames(mat)=c("true","post expect")
-cat(" mu and posterior expectation of mu",fill=TRUE)
-print(mat)
-mat=rbind(tmom$sd,pmom$sd)
-rownames(mat)=c("true","post expect")
-cat(" std dev and posterior expectation of sd",fill=TRUE)
-print(mat)
-mat=rbind(as.vector(tmom$corr),as.vector(pmom$corr))
-rownames(mat)=c("true","post expect")
-cat(" corr and posterior expectation of corr",fill=TRUE)
-print(t(mat))
+cat("Summary of Delta draws",fill=TRUE)
+summary(out$Deltadraw,tvalues=as.vector(Delta))
+cat("Summary of Normal Mixture Distribution",fill=TRUE)
+summary(out$nmix)
if(0) {
-## set if(1) to produce plots
-par(mfrow=c(4,1))
-plot(out$betadraw[1,1,])
-abline(h=simlgtdata[[1]]$beta[1])
-plot(out$betadraw[2,1,])
-abline(h=simlgtdata[[2]]$beta[1])
-plot(out$betadraw[100,3,])
-abline(h=simlgtdata[[100]]$beta[3])
-plot(out$betadraw[101,3,])
-abline(h=simlgtdata[[101]]$beta[3])
-par(mfrow=c(4,1))
-plot(out$Deltadraw[,1])
-abline(h=Delta[1,1])
-plot(out$Deltadraw[,2])
-abline(h=Delta[2,1])
-plot(out$Deltadraw[,3])
-abline(h=Delta[3,1])
-plot(out$Deltadraw[,6])
-abline(h=Delta[3,2])
-begin=1000/keep
-end=R/keep
-ngrid=50
-grid=matrix(0,ngrid,ncoef)
-rgm=matrix(c(-3,-7,-10,3,1,0),ncol=2)
-for(i in 1:ncoef) {rg=rgm[i,]; grid[,i]=rg[1] + ((1:ngrid)/ngrid)*(rg[2]-rg[1])}
-mden=eMixMargDen(grid,out$probdraw[begin:end,],out$compdraw[begin:end])
-par(mfrow=c(2,ncoef))
-for(i in 1:ncoef)
-{plot(grid[,i],mden[,i],type="l")}
-for(i in 1:ncoef)
-tden=mixDen(grid,pvec,comps)
-for(i in 1:ncoef)
-{plot(grid[,i],tden[,i],type="l")}
-}
+## plotting examples
+plot(out$betadraw)
+plot(out$nmix)
}
}
-\keyword{ models }
+
+\keyword{models}
diff --git a/man/rhierNegbinRw.Rd b/man/rhierNegbinRw.Rd
index fd41f69..3f09a34 100755
--- a/man/rhierNegbinRw.Rd
+++ b/man/rhierNegbinRw.Rd
@@ -81,8 +81,7 @@ rhierNegbinRw(Data, Prior, Mcmc)
}
\examples{
##
-if(nchar(Sys.getenv("LONG_TEST")) != 0)
-{
+if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=2000} else {R=10}
##
set.seed(66)
simnegbin =
@@ -119,37 +118,23 @@ for (i in 1:nreg) {
Beta = NULL
for (i in 1:nreg) {Beta=rbind(Beta,matrix(simnegbindata[[i]]$beta,nrow=1))}
-Data = list(regdata=simnegbindata, Z=Z)
-Deltabar = matrix(rep(0,nvar*nz),nrow=nz)
-Vdelta = 0.01 * diag(nvar)
-nu = nvar+3
-V = 0.01*diag(nvar)
-a = 0.5
-b = 0.1
-Prior = list(Deltabar=Deltabar, Vdelta=Vdelta, nu=nu, V=V, a=a, b=b)
+Data1 = list(regdata=simnegbindata, Z=Z)
+Mcmc1 = list(R=R)
-R=10000
-keep =1
-s_beta=2.93/sqrt(nvar)
-s_alpha=2.93
-c=2
-Mcmc = list(R=R, keep = keep, s_beta=s_beta, s_alpha=s_alpha, c=c)
-out = rhierNegbinRw(Data, Prior, Mcmc)
+out = rhierNegbinRw(Data=Data1, Mcmc=Mcmc1)
-# Unit level mean beta parameters
-Mbeta = matrix(rep(0,nreg*nvar),nrow=nreg)
-ndraws = length(out$alphadraw)
-for (i in 1:nreg) { Mbeta[i,] = rowSums(out$Betadraw[i, , ])/ndraws }
+cat("Summary of Delta draws",fill=TRUE)
+summary(out$Deltadraw,tvalues=as.vector(Delta))
+cat("Summary of Vbeta draws",fill=TRUE)
+summary(out$Vbetadraw,tvalues=as.vector(Vbeta[upper.tri(Vbeta,diag=TRUE)]))
+cat("Summary of alpha draws",fill=TRUE)
+summary(out$alpha,tvalues=alpha)
-cat(" Deltadraws ",fill=TRUE)
-mat=apply(out$Deltadraw,2,quantile,probs=c(.01,.05,.5,.95,.99))
-mat=rbind(as.vector(Delta),mat); rownames(mat)[1]="Delta"; print(mat)
-cat(" Vbetadraws ",fill=TRUE)
-mat=apply(out$Vbetadraw,2,quantile,probs=c(.01,.05,.5,.95,.99))
-mat=rbind(as.vector(Vbeta),mat); rownames(mat)[1]="Vbeta"; print(mat)
-cat(" alphadraws ",fill=TRUE)
-mat=apply(matrix(out$alphadraw),2,quantile,probs=c(.01,.05,.5,.95,.99))
-mat=rbind(as.vector(alpha),mat); rownames(mat)[1]="alpha"; print(mat)
+if(0){
+## plotting examples
+plot(out$betadraw)
+plot(out$alpha,tvalues=alpha)
+plot(out$Deltadraw,tvalues=as.vector(Delta))
}
}
diff --git a/man/rivGibbs.Rd b/man/rivGibbs.Rd
index 83c249c..35693b2 100755
--- a/man/rivGibbs.Rd
+++ b/man/rivGibbs.Rd
@@ -79,21 +79,20 @@ Sigma = matrix(c(1,rho,rho,1),ncol=2)
delta = c(1,4); beta = .5; gamma = c(1)
simiv = simIV(delta,beta,Sigma,n,z,w,gamma)
-Mcmc=list(); Prior=list(); Data = list()
-Data$z = z; Data$w=w; Data$x=simiv$x; Data$y=simiv$y
-Mcmc$R = R
-Mcmc$keep=1
-out=rivGibbs(Data=Data,Prior=Prior,Mcmc=Mcmc)
+Mcmc1=list(); Data1 = list()
+Data1$z = z; Data1$w=w; Data1$x=simiv$x; Data1$y=simiv$y
+Mcmc1$R = R
+Mcmc1$keep=1
+out=rivGibbs(Data=Data1,Mcmc=Mcmc1)
-cat(" deltadraws ",fill=TRUE)
-mat=apply(out$deltadraw,2,quantile,probs=c(.01,.05,.5,.95,.99))
-mat=rbind(delta,mat); rownames(mat)[1]="delta"; print(mat)
-cat(" betadraws ",fill=TRUE)
-qout=quantile(out$betadraw,probs=c(.01,.05,.5,.95,.99))
-mat=matrix(qout,ncol=1)
-mat=rbind(beta,mat); rownames(mat)=c("beta",names(qout)); print(mat)
-cat(" Sigma draws",fill=TRUE)
-mat=apply(out$Sigmadraw,2,quantile,probs=c(.01,.05,.5,.95,.99))
-mat=rbind(as.vector(Sigma),mat); rownames(mat)[1]="Sigma"; print(mat)
+cat("Summary of Beta draws",fill=TRUE)
+summary(out$betadraw,tvalues=beta)
+cat("Summary of Sigma draws",fill=TRUE)
+summary(out$Sigmadraw,tvalues=as.vector(Sigma[upper.tri(Sigma,diag=TRUE)]))
+
+if(0){
+## plotting examples
+plot(out$betadraw)
+}
}
\keyword{ models }
diff --git a/man/rmnlIndepMetrop.Rd b/man/rmnlIndepMetrop.Rd
index d289a02..f233d03 100755
--- a/man/rmnlIndepMetrop.Rd
+++ b/man/rmnlIndepMetrop.Rd
@@ -54,14 +54,40 @@ if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=2000} else {R=10}
set.seed(66)
n=200; p=3; beta=c(1,-1,1.5,.5)
+
+simmnl= function(p,n,beta) {
+ # note: create X array with 2 alt.spec vars
+ k=length(beta)
+ X1=matrix(runif(n*p,min=-1,max=1),ncol=p)
+ X2=matrix(runif(n*p,min=-1,max=1),ncol=p)
+ X=createX(p,na=2,nd=NULL,Xd=NULL,Xa=cbind(X1,X2),base=1)
+ Xbeta=X\%*\%beta # now do probs
+ p=nrow(Xbeta)/n
+ Xbeta=matrix(Xbeta,byrow=TRUE,ncol=p)
+ Prob=exp(Xbeta)
+ iota=c(rep(1,p))
+ denom=Prob\%*\%iota
+ Prob=Prob/as.vector(denom)
+ # draw y
+ y=vector("double",n)
+ ind=1:p
+ for (i in 1:n)
+ { yvec=rmultinom(1,1,Prob[i,]); y[i]=ind\%*\%yvec }
+ return(list(y=y,X=X,beta=beta,prob=Prob))
+}
+
simout=simmnl(p,n,beta)
-A=diag(c(rep(.01,length(beta)))); betabar=rep(0,length(beta))
-Data=list(y=simout$y,X=simout$X,p=p); Mcmc=list(R=R,keep=1) ; Prior=list(A=A,betabar=betabar)
-out=rmnlIndepMetrop(Data=Data,Prior=Prior,Mcmc=Mcmc)
-cat(" Betadraws ",fill=TRUE)
-mat=apply(out$betadraw,2,quantile,probs=c(.01,.05,.5,.95,.99))
-mat=rbind(beta,mat); rownames(mat)[1]="beta"; print(mat)
+Data1=list(y=simout$y,X=simout$X,p=p); Mcmc1=list(R=R,keep=1)
+out=rmnlIndepMetrop(Data=Data1,Mcmc=Mcmc1)
+
+cat("Summary of beta draws",fill=TRUE)
+summary(out$betadraw,tvalues=beta)
+
+if(0){
+## plotting examples
+plot(out$betadraw)
+}
}
\keyword{ models }
diff --git a/man/rmnpGibbs.Rd b/man/rmnpGibbs.Rd
index d46f601..ab3aee2 100755
--- a/man/rmnpGibbs.Rd
+++ b/man/rmnpGibbs.Rd
@@ -79,23 +79,41 @@ n=500
beta=c(-1,1,1,2)
Sigma=matrix(c(1,.5,.5,1),ncol=2)
k=length(beta)
-x1=runif(n*(p-1),min=-1,max=1); x2=runif(n*(p-1),min=-1,max=1)
-I2=diag(rep(1,p-1)); xadd=rbind(I2)
-for(i in 2:n) { xadd=rbind(xadd,I2)}
-X=cbind(xadd,x1,x2)
+X1=matrix(runif(n*p,min=0,max=2),ncol=p); X2=matrix(runif(n*p,min=0,max=2),ncol=p)
+X=createX(p,na=2,nd=NULL,Xa=cbind(X1,X2),Xd=NULL,DIFF=TRUE,base=p)
+
+simmnp= function(X,p,n,beta,sigma) {
+ indmax=function(x) {which(max(x)==x)}
+ Xbeta=X\%*\%beta
+ w=as.vector(crossprod(chol(sigma),matrix(rnorm((p-1)*n),ncol=n)))+ Xbeta
+ w=matrix(w,ncol=(p-1),byrow=TRUE)
+ maxw=apply(w,1,max)
+ y=apply(w,1,indmax)
+ y=ifelse(maxw < 0,p,y)
+ return(list(y=y,X=X,beta=beta,sigma=sigma))
+}
+
simout=simmnp(X,p,500,beta,Sigma)
-Data=list(p=p,y=simout$y,X=simout$X)
-Mcmc=list(R=R,keep=1)
+Data1=list(p=p,y=simout$y,X=simout$X)
+Mcmc1=list(R=R,keep=1)
-out=rmnpGibbs(Mcmc=Mcmc,Data=Data)
+out=rmnpGibbs(Data=Data1,Mcmc=Mcmc1)
-cat(" Betadraws ",fill=TRUE)
-mat=apply(out$betadraw/sqrt(out$sigmadraw[,1]),2,quantile,probs=c(.01,.05,.5,.95,.99))
-mat=rbind(beta,mat); rownames(mat)[1]="beta"; print(mat)
-cat(" Sigmadraws ",fill=TRUE)
-mat=apply(out$sigmadraw/out$sigmadraw[,1],2,quantile,probs=c(.01,.05,.5,.95,.99))
-mat=rbind(as.vector(Sigma),mat); rownames(mat)[1]="sigma"; print(mat)
+cat(" Summary of Betadraws ",fill=TRUE)
+betatilde=out$betadraw/sqrt(out$sigmadraw[,1])
+attributes(betatilde)$class="bayesm.mat"
+summary(betatilde,tvalues=beta)
+cat(" Summary of Sigmadraws ",fill=TRUE)
+sigmadraw=out$sigmadraw/out$sigmadraw[,1]
+attributes(sigmadraw)$class="bayesm.var"
+summary(sigmadraw,tvalues=as.vector(Sigma[upper.tri(Sigma,diag=TRUE)]))
+
+
+if(0){
+## plotting examples
+plot(betatilde,tvalues=beta)
+}
}
\keyword{ models }
diff --git a/man/rmultireg.Rd b/man/rmultireg.Rd
index 51ae600..b8dd997 100755
--- a/man/rmultireg.Rd
+++ b/man/rmultireg.Rd
@@ -48,7 +48,6 @@ rmultireg(Y, X, Bbar, A, nu, V)
input arguments for proper dimensions and type.
}
-\seealso{ \code{\link{rmultiregfp}},\code{\link{init.rmultiregfp}}}
\examples{
##
if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=2000} else {R=10}
diff --git a/man/rmultiregfp.Rd b/man/rmultiregfp.Rd
deleted file mode 100755
index 6889e96..0000000
--- a/man/rmultiregfp.Rd
+++ /dev/null
@@ -1,48 +0,0 @@
-\name{rmultiregfp}
-\alias{rmultiregfp}
-
-\title{ Draw from the Posterior of a Multivariate Regression }
-\description{
- \code{ rmultiregfp} draws from the posterior of a Multivariate Regression model with a natural conjugate
- prior.
-}
-\usage{
-rmultiregfp(Y, X, Fparm)
-}
-
-\arguments{
- \item{Y}{ n x m matrix of observations on m dep vars }
- \item{X}{ n x k matrix of observations on indep vars (supply intercept) }
- \item{Fparm}{ a list of prior parameters prepared by \code{init.rmultiregfp}}
-}
-\details{
- Model: \eqn{Y=XB+U}. \eqn{cov(u_i) = Sigma}. \eqn{B} is k x m matrix of coefficients.
- \eqn{Sigma} is an m x m covariance matrix.
-
- Priors: \eqn{beta} given \eqn{Sigma} \eqn{\sim}{~} \eqn{N(betabar,Sigma (x) A^{-1})}.
- \eqn{betabar=vec(Bbar)}; \eqn{beta = vec(B)}. \cr
- \eqn{Sigma} \eqn{\sim}{~} IW(nu,V).
-
- prepare Fparm by call \code{\link{init.rmultiregfp}}
-}
-\value{
- A list of the components of a draw from the posterior
- \item{B }{ draw of regression coefficient matrix }
- \item{Sigma }{ draw of \eqn{Sigma} }
-}
-\references{ For further discussion, see \emph{Bayesian Statistics and Marketing}
- by Rossi, Allenby and McCulloch. \cr
- \url{http://faculty.chicagogsb.edu/peter.rossi/research/bsm.html}
-}
-
-\author{ Peter Rossi, Graduate School of Business, University of Chicago,
- \email{Peter.Rossi at ChicagoGsb.edu}.
-}
-
-\section{Warning}{
- This routine is a utility routine that does \strong{not} check the
- input arguments for proper dimensions and type.
-}
-
-\seealso{ \code{\link{rmultireg}},\code{\link{init.rmultiregfp}}}
-\keyword{ regression }
diff --git a/man/rmvpGibbs.Rd b/man/rmvpGibbs.Rd
index 6d01487..4f59406 100755
--- a/man/rmvpGibbs.Rd
+++ b/man/rmvpGibbs.Rd
@@ -81,23 +81,39 @@ Sigma=matrix(c(1,.5,.5,.5,1,.5,.5,.5,1),ncol=3)
k=length(beta)
I2=diag(rep(1,p)); xadd=rbind(I2)
for(i in 2:n) { xadd=rbind(xadd,I2)}; X=xadd
+
+simmvp= function(X,p,n,beta,sigma) {
+ w=as.vector(crossprod(chol(sigma),matrix(rnorm(p*n),ncol=n)))+ X\%*\%beta
+ y=ifelse(w<0,0,1)
+ return(list(y=y,X=X,beta=beta,sigma=sigma))
+}
+
simout=simmvp(X,p,500,beta,Sigma)
-Data=list(p=p,y=simout$y,X=simout$X)
-Mcmc=list(R=R,keep=1)
-out=rmvpGibbs(Data=Data,Mcmc=Mcmc)
+Data1=list(p=p,y=simout$y,X=simout$X)
+Mcmc1=list(R=R,keep=1)
+out=rmvpGibbs(Data=Data1,Mcmc=Mcmc1)
ind=seq(from=0,by=p,length=k)
inda=1:3
ind=ind+inda
cat(" Betadraws ",fill=TRUE)
-mat=apply(out$betadraw/sqrt(out$sigmadraw[,ind]),2,quantile,probs=c(.01,.05,.5,.95,.99))
-mat=rbind(beta,mat); rownames(mat)[1]="beta"; print(mat)
+betatilde=out$betadraw/sqrt(out$sigmadraw[,ind])
+attributes(betatilde)$class="bayesm.mat"
+summary(betatilde,tvalues=beta/sqrt(diag(Sigma)))
+
rdraw=matrix(double((R)*p*p),ncol=p*p)
rdraw=t(apply(out$sigmadraw,1,nmat))
+attributes(rdraw)$class="bayesm.var"
+tvalue=nmat(as.vector(Sigma))
+dim(tvalue)=c(p,p)
+tvalue=as.vector(tvalue[upper.tri(tvalue,diag=TRUE)])
cat(" Draws of Correlation Matrix ",fill=TRUE)
-mat=apply(rdraw,2,quantile,probs=c(.01,.05,.5,.95,.99))
-mat=rbind(as.vector(Sigma),mat); rownames(mat)[1]="Sigma"; print(mat)
+summary(rdraw,tvalues=tvalue)
+
+if(0){
+plot(betatilde,tvalues=beta/sqrt(diag(Sigma)))
+}
}
\keyword{ models }
diff --git a/man/rnegbinRw.Rd b/man/rnegbinRw.Rd
index 4fd6434..f9df16a 100755
--- a/man/rnegbinRw.Rd
+++ b/man/rnegbinRw.Rd
@@ -93,23 +93,21 @@ simnegbindata = NULL
beta = c(0.6,0.2)
X = cbind(rep(1,nobs),rnorm(nobs,mean=2,sd=0.5))
simnegbindata = list(y=simnegbin(X,beta,alpha), X=X, beta=beta)
-Data = simnegbindata
-betabar = rep(0,nvar)
-A = 0.01 * diag(nvar)
-a = 0.5; b = 0.1
-Prior = list(betabar=betabar, A=A, a=a, b=b)
-
-keep =1
-s_beta=2.93/sqrt(nvar); s_alpha=2.93
-Mcmc = list(R=R, keep = keep, s_beta=s_beta, s_alpha=s_alpha)
-out = rnegbinRw(Data, Prior, Mcmc)
-
-cat(" betadraws ",fill=TRUE)
-mat=apply(out$betadraw,2,quantile,probs=c(.01,.05,.5,.95,.99))
-mat=rbind(beta,mat); rownames(mat)[1]="beta"; print(mat)
-cat(" alphadraws ",fill=TRUE)
-mat=apply(matrix(out$alphadraw),2,quantile,probs=c(.01,.05,.5,.95,.99))
-mat=rbind(as.vector(alpha),mat); rownames(mat)[1]="alpha"; print(mat)
+
+Data1 = simnegbindata
+Mcmc1 = list(R=R)
+
+out = rnegbinRw(Data=Data1,Mcmc=Mcmc1)
+
+cat("Summary of alpha/beta draw",fill=TRUE)
+summary(out$alphadraw,tvalues=alpha)
+summary(out$betadraw,tvalues=beta)
+
+if(0){
+## plotting examples
+plot(out$betadraw)
+}
+
}
\keyword{ models }
diff --git a/man/rnmixGibbs.Rd b/man/rnmixGibbs.Rd
index 598b03a..5f906da 100755
--- a/man/rnmixGibbs.Rd
+++ b/man/rnmixGibbs.Rd
@@ -49,13 +49,15 @@ rnmixGibbs(Data, Prior, Mcmc)
}
}
\value{
- a list containing:
+ \item{nmix} {a list containing: probdraw,zdraw,compdraw}
+}
+
+\note{
+ more details on contents of nmix: \cr
\item{probdraw}{R/keep x ncomp array of mixture prob draws}
\item{zdraw}{R/keep x nobs array of indicators of mixture comp identity for each obs}
\item{compdraw}{R/keep lists of lists of comp parm draws}
-}
-\note{
In this model, the component normal parameters are not-identified due to label-switching.
However, the fitted mixture of normals density is identified as it is invariant to label-switching.
See Allenby et al, chapter 5 for details. Use \code{eMixMargDen} or \code{momMix} to compute
@@ -90,71 +92,25 @@ pvec=(1:k)/sum(1:k)
nobs=500
dm = rmixture(nobs,pvec,comps)
-Data=list(y=dm$x)
+Data1=list(y=dm$x)
ncomp=9
-Prior=list(ncomp=ncomp)
-Mcmc=list(R=R,keep=1)
-out=rnmixGibbs(Data=Data,Prior=Prior,Mcmc=Mcmc)
+Prior1=list(ncomp=ncomp)
+Mcmc1=list(R=R,keep=1)
+out=rnmixGibbs(Data=Data1,Prior=Prior1,Mcmc=Mcmc1)
+cat("Summary of Normal Mixture Distribution",fill=TRUE)
+summary(out)
tmom=momMix(matrix(pvec,nrow=1),list(comps))
-if(R < 1000) {begin=1} else {begin=500}
-pmom=momMix(out$probdraw[begin:R,],out$compdraw[begin:R])
-mat=rbind(tmom$mu,pmom$mu)
-rownames(mat)=c("true","post expect")
-cat(" mu and posterior expectation of mu",fill=TRUE)
-print(mat)
-mat=rbind(tmom$sd,pmom$sd)
-rownames(mat)=c("true","post expect")
-cat(" std dev and posterior expectation of sd",fill=TRUE)
+mat=rbind(tmom$mu,tmom$sd)
+cat(" True Mean/Std Dev",fill=TRUE)
print(mat)
-mat=rbind(as.vector(tmom$corr),as.vector(pmom$corr))
-rownames(mat)=c("true","post expect")
-cat(" corr and posterior expectation of corr",fill=TRUE)
-print(t(mat))
if(0){
##
## plotting examples
##
-## check true and estimated marginal densities
-end=R
-grid=NULL
-for (i in 1:dim){
- rgi=range(dm$x[,i])
- gr=seq(from=rgi[1],to=rgi[2],length.out=50)
- grid=cbind(grid,gr)
+plot(out)
}
-
-tmden=mixDen(grid,pvec,comps)
-pmden=eMixMargDen(grid,out$probdraw[begin:end,],out$compdraw[begin:end])
-
-## plot the marginal on third variable
-plot(range(grid[,3]),c(0,1.1*max(tmden[,3],pmden[,3])),type="n",xlab="",ylab="density")
-lines(grid[,3],tmden[,3],col="blue",lwd=2)
-lines(grid[,3],pmden[,3],col="red",lwd=2)
-
-## compute implied bivariate marginal distributions
-i=1
-j=2
-rxi=range(dm$x[,1])
-rxj=range(dm$x[,2])
-xi=seq(from=rxi[1],to=rxi[2],length.out=50)
-xj=seq(from=rxj[1],to=rxj[2],length.out=50)
-den=matrix(0,ncol=length(xi),nrow=length(xj))
-for(ind in as.integer(seq(from=begin,to=end,length.out=100))){
- den=den+mixDenBi(i,j,xi,xj,out$probdraw[ind,],out$compdraw[[ind]])
-}
-tden=matrix(0,ncol=length(xi),nrow=length(xj))
-tden=mixDenBi(i,j,xi,xj,pvec,comps)
-
-par(mfrow=c(2,1))
-image(xi,xj,tden,col=terrain.colors(100),xlab="",ylab="")
-contour(xi,xj,tden,add=TRUE,drawlabels=FALSE)
-title("True Bivariate Marginal")
-image(xi,xj,den,col=terrain.colors(100),xlab="",ylab="")
-contour(xi,xj,den,add=TRUE,drawlabels=FALSE)
-title("Posterior Mean of Bivariate Marginal")
-}
-
+
}
\keyword{ multivariate }
diff --git a/man/rordprobitGibbs.Rd b/man/rordprobitGibbs.Rd
new file mode 100755
index 0000000..6dce63b
--- /dev/null
+++ b/man/rordprobitGibbs.Rd
@@ -0,0 +1,115 @@
+\name{rordprobitGibbs}
+\alias{rordprobitGibbs}
+\concept{bayes}
+\concept{MCMC}
+\concept{probit}
+\concept{Gibbs Sampling}
+
+\title{ Gibbs Sampler for Ordered Probit }
+\description{
+ \code{rordprobitGibbs} implements a Gibbs Sampler for the ordered probit model.
+
+}
+\usage{
+rordprobitGibbs(Data, Prior, Mcmc)
+}
+
+\arguments{
+ \item{Data}{ list(X, y, k)}
+ \item{Prior}{ list(betabar, A, dstarbar, Ad)}
+ \item{Mcmc}{ list(R, keep, s, change, draw) }
+}
+
+\details{
+ Model: \eqn{z = X\beta + e}. \eqn{e} \eqn{\sim}{~} \eqn{N(0,I)}.
+ y=1,..,k. cutoff=c( c [1] ,..c [k+1] ). \cr
+ y=k, if c [k] <= z < c [k+1] .
+
+ Prior: \eqn{\beta} \eqn{\sim}{~} \eqn{N(betabar,A^{-1})}.
+ \eqn{dstar} \eqn{\sim}{~} \eqn{N(dstarbar,Ad^{-1})}.
+
+ List arguments contain
+ \describe{
+ \item{\code{X}}{n x nvar Design Matrix}
+ \item{\code{y}}{n x 1 vector of observations, (1,...,k)}
+ \item{\code{k}}{the largest possible value of y}
+ \item{\code{betabar}}{nvar x 1 prior mean (def: 0)}
+ \item{\code{A}}{nvar x nvar prior precision matrix (def: .01I)}
+ \item{\code{dstarbar}}{ndstar x 1 prior mean, ndstar=k-2 (def: 0)}
+ \item{\code{Ad}}{ndstar x ndstar prior precision matrix (def:I)}
+ \item{\code{s}}{ scaling parm for RW Metropolis (def: 2.93/sqrt(nvar))}
+ \item{\code{R}}{ number of MCMC draws }
+ \item{\code{keep}}{ thinning parameter - keep every keepth draw (def: 1)}
+ }
+}
+
+\value{
+ \item{betadraw }{R/keep x k matrix of betadraws}
+ \item{cutdraw }{R/keep x (k-1) matrix of cutdraws}
+ \item{dstardraw }{R/keep x (k-2) matrix of dstardraws}
+ \item{accept }{a value of acceptance rate in RW Metropolis}
+}
+\note{
+ set c[1]=-100. c[k+1]=100. c[2] is set to 0 for identification. \cr
+
+ The relationship between cut-offs and dstar is \cr
+ c[3] = exp(dstar[1]), c[4]=c[3]+exp(dstar[2]),..., c[k] = c[k-1] + exp(datsr[k-2])
+
+ Be careful in assessing prior parameter, Ad. .1 is too small for many applications.
+}
+
+\references{ \emph{Bayesian Statistics and Marketing}
+ by Rossi, Allenby and McCulloch\cr
+ \url{http://faculty.chicagogsb.edu/peter.rossi/research/bsm.html}
+}
+
+\author{ Peter Rossi, Graduate School of Business, University of Chicago,
+ \email{Peter.Rossi at ChicagoGsb.edu}.
+}
+
+\seealso{ \code{\link{rbprobitGibbs}} }
+\examples{
+##
+## rordprobitGibbs example
+##
+if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=2000} else {R=10}
+
+## simulate data for ordered probit model
+
+ simordprobit=function(X, betas, cutoff){
+ z = X\%*\%betas + rnorm(nobs)
+ y = cut(z, br = cutoff, right=TRUE, include.lowest = TRUE, labels = FALSE)
+ return(list(y = y, X = X, k=(length(cutoff)-1), betas= betas, cutoff=cutoff ))
+ }
+
+ set.seed(66)
+ nobs=300
+ X=cbind(rep(1,nobs),runif(nobs, min=0, max=5),runif(nobs,min=0, max=5))
+ k=5
+ betas=c(0.5, 1, -0.5)
+ cutoff=c(-100, 0, 1.0, 1.8, 3.2, 100)
+ simout=simordprobit(X, betas, cutoff)
+ Data=list(X=simout$X,y=simout$y, k=k)
+
+## set Mcmc for ordered probit model
+
+ Mcmc=list(R=R)
+ out=rordprobitGibbs(Data=Data,Mcmc=Mcmc)
+
+ cat(" ", fill=TRUE)
+ cat("acceptance rate= ",accept=out$accept,fill=TRUE)
+
+## outputs of betadraw and cut-off draws
+
+ cat(" Summary of betadraws",fill=TRUE)
+ summary(out$betadraw,tvalues=betas)
+ cat(" Summary of cut-off draws",fill=TRUE)
+ summary(out$cutdraw,tvalues=cutoff[2:k])
+
+if(0){
+## plotting examples
+plot(out$cutdraw)
+}
+
+}
+\keyword{ models }
diff --git a/man/rscaleUsage.Rd b/man/rscaleUsage.Rd
index f3ff4b0..2d7a030 100755
--- a/man/rscaleUsage.Rd
+++ b/man/rscaleUsage.Rd
@@ -72,13 +72,11 @@ if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=1000} else {R=5}
data(customerSat)
surveydat = list(k=10,x=as.matrix(customerSat))
-mcmc = list(R=R)
+Mcmc1 = list(R=R)
set.seed(66)
-out=rscaleUsage(Data=surveydat,Mcmc=mcmc)
+out=rscaleUsage(Data=surveydat,Mcmc=Mcmc1)
-cat(" mudraws ",fill=TRUE)
-mat=apply(out$mudraw,2,quantile,probs=c(.01,.05,.5,.95,.99))
-print(mat)
+summary(out$mudraw)
}
}
diff --git a/man/rsurGibbs.Rd b/man/rsurGibbs.Rd
index 53e0c70..da0f5a8 100755
--- a/man/rsurGibbs.Rd
+++ b/man/rsurGibbs.Rd
@@ -81,18 +81,20 @@ y2=X2\%*\%beta2+E[,2]
regdata=NULL
regdata[[1]]=list(y=y1,X=X1)
regdata[[2]]=list(y=y2,X=X2)
-Mcmc=list(R=R)
-out=rsurGibbs(Data=list(regdata=regdata),Mcmc=Mcmc)
-##
-## summarize GS output
-if(R < 100) {begin=1} else {begin=100}
-end=Mcmc$R
-cat(" betadraws ",fill=TRUE)
-mat=apply(out$betadraw[begin:end,],2,quantile,probs=c(.01,.05,.5,.95,.99))
-mat=rbind(c(beta1,beta2),mat); rownames(mat)[1]="beta"; print(mat)
-cat(" Sigmadraws ",fill=TRUE)
-mat=apply(out$Sigmadraw[begin:end,],2,quantile,probs=c(.01,.05,.5,.95,.99))
-mat=rbind(as.vector(Sigma),mat); rownames(mat)[1]="sigma"; print(mat)
+
+Mcmc1=list(R=R)
+
+out=rsurGibbs(Data=list(regdata=regdata),Mcmc=Mcmc1)
+
+cat("Summary of beta draws",fill=TRUE)
+summary(out$betadraw,tvalues=c(beta1,beta2))
+cat("Summary of Sigmadraws",fill=TRUE)
+summary(out$Sigmadraw,tvalues=as.vector(Sigma[upper.tri(Sigma,diag=TRUE)]))
+
+if(0){
+plot(out$betadraw,tvalues=c(beta1,beta2))
+}
+
}
\keyword{ regression}
diff --git a/man/runireg.Rd b/man/runireg.Rd
index 9eab58b..67b064b 100755
--- a/man/runireg.Rd
+++ b/man/runireg.Rd
@@ -55,11 +55,15 @@ X=cbind(rep(1,n),runif(n)); beta=c(1,2); sigsq=.25
y=X\%*\%beta+rnorm(n,sd=sqrt(sigsq))
out=runireg(Data=list(y=y,X=X),Mcmc=list(R=R))
-cat(" betadraws ",fill=TRUE)
-mat=apply(out$betadraw,2,quantile,probs=c(.01,.05,.5,.95,.99))
-mat=rbind(beta,mat); rownames(mat)[1]="beta"; print(mat)
-cat(" Sigma-sq draws",fill=TRUE)
-cat(" sigma-sq= ",sigsq,fill=TRUE)
-print(quantile(out$sigmasqdraw,probs=c(.01,.05,.5,.95,.99)))
+
+cat("Summary of beta/sigma-sq draws",fill=TRUE)
+summary(out$betadraw,tvalues=beta)
+summary(out$sigmasqdraw,tvalues=sigsq)
+
+if(0){
+## plotting examples
+plot(out$betadraw)
+}
+
}
\keyword{ regression }
diff --git a/man/runiregGibbs.Rd b/man/runiregGibbs.Rd
index ab04c83..7b4d1c0 100755
--- a/man/runiregGibbs.Rd
+++ b/man/runiregGibbs.Rd
@@ -55,16 +55,18 @@ n=100
X=cbind(rep(1,n),runif(n)); beta=c(1,2); sigsq=.25
y=X\%*\%beta+rnorm(n,sd=sqrt(sigsq))
-A=diag(c(.05,.05)); betabar=c(0,0)
-nu=3; ssq=1.0
+Data1=list(y=y,X=X); Mcmc1=list(R=R)
+
+out=runiregGibbs(Data=Data1,Mcmc=Mcmc1)
+
+cat("Summary of beta and Sigma draws",fill=TRUE)
+summary(out$betadraw,tvalues=beta)
+summary(out$sigmasqdraw,tvalues=sigsq)
+
+if(0){
+## plotting examples
+plot(out$betadraw)
+}
-Data=list(y=y,X=X); Mcmc=list(R=R,keep=1) ; Prior=list(A=A,betabar=betabar,nu=nu,ssq=ssq)
-out=runiregGibbs(Data=Data,Prior=Prior,Mcmc=Mcmc)
-cat(" betadraws ",fill=TRUE)
-mat=apply(out$betadraw,2,quantile,probs=c(.01,.05,.5,.95,.99))
-mat=rbind(beta,mat); rownames(mat)[1]="beta"; print(mat)
-cat(" Sigma-sq draws",fill=TRUE)
-cat(" sigma-sq= ",sigsq,fill=TRUE)
-print(quantile(out$sigmasqdraw,probs=c(.01,.05,.5,.95,.99)))
}
\keyword{ regression }
diff --git a/man/simmnl.Rd b/man/simmnl.Rd
deleted file mode 100755
index 1c097cc..0000000
--- a/man/simmnl.Rd
+++ /dev/null
@@ -1,42 +0,0 @@
-\name{simmnl}
-\alias{simmnl}
-\concept{logit}
-\concept{mnl}
-\title{ Simulate from Multinomial Logit Model }
-\description{
- \code{simmnl} simulates from the MNL model.
-}
-\usage{
-simmnl(p, n, beta)
-}
-\arguments{
- \item{p}{ number choice alternatives }
- \item{n}{ number of observations }
- \item{beta}{ MNL coefficient vector }
-}
-\details{
- \code{simmnl} will simulate two uniformly distributed X vars and add intercepts.
-}
-
-\value{
- \item{y}{ n x 1 vector of multinomial outcomes (1, \ldots, p)}
- \item{X} { n x k design matrix }
- \item{beta}{beta vector }
- \item{prob}{n x j array of choice probabilities }
-}
-\references{ For further discussion, see \emph{Bayesian Statistics and Marketing}
- by Rossi, Allenby, and McCulloch. \cr
- \url{http://faculty.chicagogsb.edu/peter.rossi/research/bsm.html}
-}
-
-\author{ Peter Rossi, Graduate School of Business, University of Chicago,
- \email{Peter.Rossi at ChicagoGsb.edu}.
-}
-\section{Warning}{
- This routine is a utility routine that does \strong{not} check the
- input arguments for proper dimensions and type.
-}
-
-\seealso{ \code{\link{llmnl}}, \code{\link{rmnlIndepMetrop}} }
-}
-\keyword{ models }
diff --git a/man/simmnlwX.Rd b/man/simmnlwX.Rd
deleted file mode 100755
index 07dfc9a..0000000
--- a/man/simmnlwX.Rd
+++ /dev/null
@@ -1,36 +0,0 @@
-\name{simmnlwX}
-\alias{simmnlwX}
-\title{ Simulate from MNL given X Matrix }
-\description{
- \code{simmnlwX} simulates from MNL given X Matrix.
-}
-\usage{
-simmnlwX(n, X, beta)
-}
-\arguments{
- \item{n}{ number of obs }
- \item{X}{ n*p x k Design matrix (p is number of choice alts)}
- \item{beta}{ k x 1 coefficient vector }
-}
-\value{
- a list containing:
- \item{y}{ n x 1 vector of multinomial outcomes (1, \ldots, nrow(X)/n)}
- \item{X}{ Design matrix }
- \item{beta}{ coefficient vector }
- \item{prob}{ n x p array of choice probs }
-}
-\references{ For further discussion, see \emph{Bayesian Statistics and Marketing}
- by Rossi, Allenby and McCulloch. \cr
- \url{http://faculty.chicagogsb.edu/peter.rossi/research/bsm.html}
-}
-
-\author{ Peter Rossi, Graduate School of Business, University of Chicago,
- \email{Peter.Rossi at ChicagoGsb.edu}.
-}
-
-\section{Warning}{
- This routine is a utility routine that does \strong{not} check the
- input arguments for proper dimensions and type.
-}
-\seealso{ \code{\link{simmnl}} }
-\keyword{ models }
diff --git a/man/simmnp.Rd b/man/simmnp.Rd
deleted file mode 100755
index 0245f09..0000000
--- a/man/simmnp.Rd
+++ /dev/null
@@ -1,44 +0,0 @@
-\name{simmnp}
-\alias{simmnp}
-\concept{simulation}
-\concept{probit}
-\concept{mnp}
-\title{ Simulate from Multinomial Probit Model }
-\description{
- \code{simmvp} simulates from the multinomial probit model.
-}
-\usage{
-simmnp(X, p, n, beta, sigma)
-}
-\arguments{
- \item{X}{ n*(p-1) x length(beta) Design matrix }
- \item{p}{ number of choice alternatives }
- \item{n}{ number of observations }
- \item{beta}{ coefficient vector }
- \item{sigma}{ (p-1) x (p-1) covariance matrix }
-}
-\value{
- a list of
- \item{y}{n vector of multinomial (1, \ldots,p) outcomes}
- \item{X}{Design matrix}
- \item{beta}{ coefficients}
- \item{sigma}{ covariance matrix }
-}
-\references{ For further discussion, see \emph{Bayesian Statistics and Marketing}
- by Rossi, Allenby and McCulloch, Chapter 4. \cr
- \url{http://faculty.chicagogsb.edu/peter.rossi/research/bsm.html}
-}
-
-\author{ Peter Rossi, Graduate School of Business, University of Chicago,
- \email{Peter.Rossi at ChicagoGsb.edu}.
-}
-
-\section{Warning}{
- This routine is a utility routine that does \strong{not} check the
- input arguments for proper dimensions and type.
-}
-
-\seealso{ \code{\link{rmnpGibbs}} }
-
-}
-\keyword{ models }
diff --git a/man/simmvp.Rd b/man/simmvp.Rd
deleted file mode 100755
index 8c06042..0000000
--- a/man/simmvp.Rd
+++ /dev/null
@@ -1,43 +0,0 @@
-\name{simmvp}
-\alias{simmvp}
-\concept{probit}
-\concept{mvp}
-\title{ Simulate from Multivariate Probit Model }
-\description{
- \code{simmvp} simulates from the multivariate probit model.
-}
-\usage{
-simmvp(X, p, n, beta, sigma)
-}
-\arguments{
- \item{X}{ n*p x length(beta) Design matrix }
- \item{p}{ dimension of the MVP}
- \item{n}{ number of observations }
- \item{beta}{ coefficient vector }
- \item{sigma}{ p x p covariance matrix }
-}
-\value{
- a list of
- \item{y}{p*n vector of 0/1 binary outcomes}
- \item{X}{Design matrix}
- \item{beta}{ coefficients}
- \item{sigma}{ covariance matrix }
-}
-\references{ For further discussion, see \emph{Bayesian Statistics and Marketing}
- by Rossi, Allenby and McCulloch, Chapter 4. \cr
- \url{http://faculty.chicagogsb.edu/peter.rossi/research/bsm.html}
-}
-
-\author{ Peter Rossi, Graduate School of Business, University of Chicago,
- \email{Peter.Rossi at ChicagoGsb.edu}.
-}
-
-\section{Warning}{
- This routine is a utility routine that does \strong{not} check the
- input arguments for proper dimensions and type.
-}
-
-\seealso{ \code{\link{rmvpGibbs}} }
-
-}
-\keyword{ models }
diff --git a/man/summary.bayesm.mat.Rd b/man/summary.bayesm.mat.Rd
new file mode 100755
index 0000000..c6adc78
--- /dev/null
+++ b/man/summary.bayesm.mat.Rd
@@ -0,0 +1,42 @@
+\name{summary.bayesm.mat}
+\alias{summary.bayesm.mat}
+\title{Summarize Mcmc Parameter Draws }
+\description{
+ \code{summary.bayesm.mat} is an S3 method to summarize marginal distributions given an
+ array of draws
+}
+\usage{
+\method{summary}{bayesm.mat}(object, names, burnin = trunc(0.1 * nrow(X)), tvalues, QUANTILES = TRUE, TRAILER = TRUE,...)
+}
+\arguments{
+ \item{object}{ \code{object} (hereafter \code{X}) is an array of draws, usually an object of class "bayesm.mat" }
+ \item{names}{ optional character vector of names for the columns of \code{X}}
+ \item{burnin}{ number of draws to burn-in, def: .1*nrow(X) }
+ \item{tvalues}{ optional vector of "true" values for use in simulation examples }
+ \item{QUANTILES}{ logical for should quantiles be displayed, def: TRUE }
+ \item{TRAILER}{ logical for should a trailer be displayed, def: TRUE }
+ \item{...}{ optional arguments for generic function }
+}
+\details{
+ Typically, \code{summary.bayesm.nmix} will be invoked by a call to the generic summary function as in
+ \code{summary(object)} where object is of class bayesm.mat. Mean, Std Dev, Numerical Standard error (of
+ estimate of posterior mean), relative numerical efficiency (see \code{numEff}) and effective sample
+ size are displayed. If QUANTILES=TRUE, quantiles of marginal distirbutions in the columns of X are displayed. \cr
+ \cr
+ \code{summary.bayesm.mat} is also exported for direct use as a standard function, as in
+ \code{summary.bayesm.mat(matrix)}.
+}
+\author{ Peter Rossi, Graduate School of Business, University of Chicago,
+ \email{Peter.Rossi at ChicagoGsb.edu}.
+}
+\seealso{ \code{\link{summary.bayesm.var}}, \code{\link{summary.bayesm.nmix}}}
+\examples{
+##
+## not run
+# out=rmnpGibbs(Data,Prior,Mcmc)
+# summary(out$betadraw)
+#
+
+}
+\keyword{ univar }
+
diff --git a/man/summary.bayesm.nmix.Rd b/man/summary.bayesm.nmix.Rd
new file mode 100755
index 0000000..ada7fab
--- /dev/null
+++ b/man/summary.bayesm.nmix.Rd
@@ -0,0 +1,38 @@
+\name{summary.bayesm.nmix}
+\alias{summary.bayesm.nmix}
+\title{Summarize Draws of Normal Mixture Components }
+\description{
+ \code{summary.bayesm.nmix} is an S3 method to display summaries of the distribution implied
+ by draws of Normal Mixture Components. Posterior means and Variance-Covariance matrices are
+ displayed.\cr
+ \cr
+ Note: 1st and 2nd moments may not be very interpretable for mixtures of normals. This summary function
+ can take a minute or so. The current implementation is not efficient.
+}
+\usage{
+ \method{summary}{bayesm.nmix}(object, names,burnin = trunc(0.1 * nrow(probdraw)), ...)
+}
+\arguments{
+ \item{object}{ an object of class "bayesm.nmix" -- a list of lists of draws}
+ \item{names}{ optional character vector of names fo reach dimension of the density}
+ \item{burnin}{ number of draws to burn-in, def: .1*nrow(probdraw)}
+ \item{...}{ parms to send to summary}
+}
+\details{
+ an object of class "bayesm.nmix" is a list of three components:
+ \item{probdraw}{ a matrix of R/keep rows by dim of normal mix of mixture prob draws}
+ \item{second comp}{ not used}
+ \item{compdraw}{ list of list of lists with draws of mixture comp parms}
+}
+\author{ Peter Rossi, Graduate School of Business, University of Chicago,
+ \email{Peter.Rossi at ChicagoGsb.edu}.
+}
+\seealso{ \code{\link{summary.bayesm.mat}}, \code{\link{summary.bayesm.var}}}
+\examples{
+##
+## not run
+# out=rnmix(Data,Prior,Mcmc)
+# summary(out)
+#
+}
+\keyword{ hplot }
diff --git a/man/summary.bayesm.var.Rd b/man/summary.bayesm.var.Rd
new file mode 100755
index 0000000..abd0405
--- /dev/null
+++ b/man/summary.bayesm.var.Rd
@@ -0,0 +1,43 @@
+\name{summary.bayesm.var}
+\alias{summary.bayesm.var}
+\title{Summarize Draws of Var-Cov Matrices}
+\description{
+ \code{summary.bayesm.var} is an S3 method to summarize marginal distributions given an
+ array of draws
+}
+\usage{
+\method{summary}{bayesm.var}(object, names, burnin = trunc(0.1 * nrow(Vard)), tvalues, QUANTILES = FALSE , ...)
+}
+\arguments{
+ \item{object}{ \code{object} (herafter, \code{Vard}) is an array of draws of a covariance matrix }
+ \item{names}{ optional character vector of names for the columns of \code{Vard}}
+ \item{burnin}{ number of draws to burn-in, def: .1*nrow(Vard) }
+ \item{tvalues}{ optional vector of "true" values for use in simulation examples }
+ \item{QUANTILES}{ logical for should quantiles be displayed, def: TRUE }
+ \item{...}{ optional arguments for generic function }
+}
+\details{
+ Typically, \code{summary.bayesm.var} will be invoked by a call to the generic summary function as in
+ summary(object) where object is of class bayesm.var. Mean, Std Dev, Numerical Standard error (of
+ estimate of posterior mean), relative numerical efficiency (see \code{numEff}) and effective sample
+ size are displayed. If QUANTILES=TRUE, quantiles of marginal distirbutions in the columns of Vard
+ are displayed. \cr
+ \cr
+ \code{Vard} is an array of draws of a covariance matrix stored as vectors. Each row is a different draw.
+ \cr
+ The posterior mean of the vector of standard deviations and the correlation matrix are also displayed
+}
+\author{ Peter Rossi, Graduate School of Business, University of Chicago,
+ \email{Peter.Rossi at ChicagoGsb.edu}.
+}
+\seealso{ \code{\link{summary.bayesm.mat}}, \code{\link{summary.bayesm.nmix}}}
+\examples{
+##
+## not run
+# out=rmnpGibbs(Data,Prior,Mcmc)
+# summary(out$sigmadraw)
+#
+
+}
+\keyword{ univar }
+
diff --git a/man/tuna.Rd b/man/tuna.Rd
new file mode 100755
index 0000000..fb8c604
--- /dev/null
+++ b/man/tuna.Rd
@@ -0,0 +1,110 @@
+\name{tuna}
+\alias{tuna}
+\docType{data}
+\title{Sales Data on Canned Tuna Sales}
+\description{
+ Volume of canned tuna sales as well as a measure of display activity, log price and log wholesale price.
+ Weekly data aggregated to the chain level.
+ Brands are seven of the top 10 UPCs in the canned tuna product category.
+}
+\usage{data(tuna)}
+\format{
+ A data frame with 338 observations on the following 30 variables.
+ \describe{
+ \item{\code{WEEK}}{a numeric vector}
+ \item{\code{MOVE1}}{unit sales of Star Kist 6 oz.}
+ \item{\code{MOVE2}}{unit sales of Chicken of the Sea 6 oz.}
+ \item{\code{MOVE3}}{unit sales of Bumble Bee Solid 6.12 oz.}
+ \item{\code{MOVE4}}{unit sales of Bumble Bee Chunk 6.12 oz.}
+ \item{\code{MOVE5}}{unit sales of Geisha 6 oz.}
+ \item{\code{MOVE6}}{unit sales of Bumble Bee Large Cans.}
+ \item{\code{MOVE7}}{unit sales of HH Chunk Lite 6.5 oz.}
+ \item{\code{NSALE1}}{a measure of display activity of Star Kist 6 oz.}
+ \item{\code{NSALE2}}{a measure of display activity of Chicken of the Sea 6 oz.}
+ \item{\code{NSALE3}}{a measure of display activity of Bumble Bee Solid 6.12 oz.}
+ \item{\code{NSALE4}}{a measure of display activity of Bumble Bee Chunk 6.12 oz.}
+ \item{\code{NSALE5}}{a measure of display activity of Geisha 6 oz.}
+ \item{\code{NSALE6}}{a measure of display activity of Bumble Bee Large Cans.}
+ \item{\code{NSALE7}}{a measure of display activity of HH Chunk Lite 6.5 oz.}
+ \item{\code{LPRICE1}}{log of price of Star Kist 6 oz.}
+ \item{\code{LPRICE2}}{log of price of Chicken of the Sea 6 oz.}
+ \item{\code{LPRICE3}}{log of price of Bumble Bee Solid 6.12 oz.}
+ \item{\code{LPRICE4}}{log of price of Bumble Bee Chunk 6.12 oz.}
+ \item{\code{LPRICE5}}{log of price of Geisha 6 oz.}
+ \item{\code{LPRICE6}}{log of price of Bumble Bee Large Cans.}
+ \item{\code{LPRICE7}}{log of price of HH Chunk Lite 6.5 oz.}
+ \item{\code{LWHPRIC1}}{log of wholesale price of Star Kist 6 oz.}
+ \item{\code{LWHPRIC2}}{log of wholesale price of Chicken of the Sea 6 oz.}
+ \item{\code{LWHPRIC3}}{log of wholesale price of Bumble Bee Solid 6.12 oz.}
+ \item{\code{LWHPRIC4}}{log of wholesale price of Bumble Bee Chunk 6.12 oz.}
+ \item{\code{LWHPRIC5}}{log of wholesale price of Geisha 6 oz.}
+ \item{\code{LWHPRIC6}}{log of wholesale price of Bumble Bee Large Cans.}
+ \item{\code{LWHPRIC7}}{log of wholesale price of HH Chunk Lite 6.5 oz.}
+ \item{\code{FULLCUST}}{total customers visits}
+ }
+}
+\source{
+ Chevalier, A. Judith, Anil K. Kashyap and Peter E. Rossi
+ (2003), "Why Don't Prices Rise During Periods of Peak Demand? Evidence from Scanner Data,"
+ \emph{The American Economic Review} , 93(1), 15-37.
+}
+\references{
+ Chapter 7, \emph{Bayesian Statistics and Marketing} by Rossi et al. \cr
+ \url{http://faculty.chicagogsb.edu/peter.rossi/research/bsm.html}
+}
+\examples{
+data(tuna)
+cat(" Quantiles of sales",fill=TRUE)
+mat=apply(as.matrix(tuna[,2:5]),2,quantile)
+print(mat)
+
+##
+## example of processing for use with rivGibbs
+##
+if(0)
+{
+ data(tuna)
+ t = dim(tuna)[1]
+ customers = tuna[,30]
+ sales = tuna[,2:8]
+ lnprice = tuna[,16:22]
+ lnwhPrice= tuna[,23:29]
+ share=sales/mean(customers)
+ shareout=as.vector(1-rowSums(share))
+ lnprob=log(share/shareout)
+
+# create w matrix
+
+ I1=as.matrix(rep(1, t))
+ I0=as.matrix(rep(0, t))
+ intercept=rep(I1, 4)
+ brand1=rbind(I1, I0, I0, I0)
+ brand2=rbind(I0, I1, I0, I0)
+ brand3=rbind(I0, I0, I1, I0)
+ w=cbind(intercept, brand1, brand2, brand3)
+
+## choose brand 1 to 4
+
+ y=as.vector(as.matrix(lnprob[,1:4]))
+ X=as.vector(as.matrix(lnprice[,1:4]))
+ lnwhPrice=as.vector(as.matrix (lnwhPrice[1:4]))
+ z=cbind(w, lnwhPrice)
+
+ Data=list(z=z, w=w, x=X, y=y)
+ Mcmc=list(R=R, keep=1)
+ set.seed(66)
+ out=rivGibbs(Data=Data,Mcmc=Mcmc)
+
+ cat(" betadraws ",fill=TRUE)
+ summary(out$betadraw)
+
+
+if(0){
+## plotting examples
+plot(out$betadraw)
+}
+}
+
+
+}
+\keyword{datasets}
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/r-cran-bayesm.git
More information about the debian-science-commits
mailing list