[r-cran-bayesm] 29/44: Import Upstream version 2.2-0
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 a93d83be696ef1519357ecf630dbabc85578185a
Author: Andreas Tille <tille at debian.org>
Date: Thu Sep 7 13:09:43 2017 +0200
Import Upstream version 2.2-0
---
DESCRIPTION | 7 +-
NAMESPACE | 2 +-
R/plot.bayesm.hcoef.R | 2 +-
R/plot.bayesm.mat.R | 5 +-
R/plot.bayesm.nmix.R | 85 ++++--
R/rDPGibbs.R | 569 ++++++++++++++++++++++++++++++++++++
R/rhierLinearMixture.R | 5 +-
R/rhierLinearModel.R | 7 +-
R/rhierNegbinRw.R | 31 +-
R/rivDP.R | 697 +++++++++++++++++++++++++++++++++++++++++++++
R/rnegbinRw.R | 1 +
R/summary.bayesm.mat.R | 2 +-
R/summary.bayesm.nmix.R | 4 +-
inst/doc/bayesm-manual.pdf | Bin 464284 -> 464797 bytes
man/plot.bayesm.mat.Rd | 3 +-
man/plot.bayesm.nmix.Rd | 18 +-
man/rDPGibbs.Rd | 212 ++++++++++++++
man/rhierNegbinRw.Rd | 4 +-
man/rivDP.Rd | 148 ++++++++++
man/rnmixGibbs.Rd | 2 +-
man/rscaleUsage.Rd | 2 +-
man/tuna.Rd | 5 +-
src/thetadraw.c | 150 ++++++++++
23 files changed, 1891 insertions(+), 70 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index 86448c6..a8dc57f 100755
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
Package: bayesm
-Version: 2.1-2
-Date: 2007-03-15
+Version: 2.2-0
+Date: 2008-03-28
Title:Bayesian Inference for Marketing/Micro-econometrics
Author: Peter Rossi <peter.rossi at ChicagoGsb.edu>,
Rob McCulloch <robert.mcculloch at ChicagoGsb.edu>.
@@ -16,6 +16,7 @@ Description: bayesm covers many important models used
Multivariate Probit,
Negative Binomial (Poisson) Regression,
Multivariate Mixtures of Normals (including clustering),
+ Dirichlet Process Prior Density Estimation with normal base
Hierarchical Linear Models with normal prior and covariates,
Hierarchical Linear Models with a mixture of normals prior and covariates,
Hierarchical Multinomial Logits with a mixture of normals prior
@@ -30,4 +31,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: Thu Mar 15 17:10:45 2007; per
+Packaged: Fri Mar 7 11:04:01 2008; per
diff --git a/NAMESPACE b/NAMESPACE
index f640295..db36bd0 100755
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -8,7 +8,7 @@ rmvpGibbs,rhierLinearModel,rhierMnlRwMixture,rivGibbs,
rmnlIndepMetrop,rscaleUsage,ghkvec,condMom,logMargDenNR,
rhierBinLogit,rnegbinRw,rhierNegbinRw,rbiNormGibbs,clusterMix,rsurGibbs,
mixDenBi,mnpProb,rhierLinearMixture,summary.bayesm.mat,plot.bayesm.mat,
-plot.bayesm.hcoef,plot.bayesm.nmix,rordprobitGibbs)
+plot.bayesm.hcoef,plot.bayesm.nmix,rordprobitGibbs,rivGibbs,rivDP,rDPGibbs)
## register S3 methods
diff --git a/R/plot.bayesm.hcoef.R b/R/plot.bayesm.hcoef.R
index b031a91..a7e6d7d 100755
--- a/R/plot.bayesm.hcoef.R
+++ b/R/plot.bayesm.hcoef.R
@@ -38,6 +38,6 @@ plot.bayesm.hcoef=function(x,burnin=trunc(.1*R),...){
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,...)
+ plot(pmeans,names,TRACEPLOT=FALSE,INT=FALSE,DEN=FALSE,CHECK_NDRAWS=FALSE,...)
invisible()
}
diff --git a/R/plot.bayesm.mat.R b/R/plot.bayesm.mat.R
index d714239..f0f180c 100755
--- a/R/plot.bayesm.mat.R
+++ b/R/plot.bayesm.mat.R
@@ -1,4 +1,5 @@
-plot.bayesm.mat=function(x,names,burnin=trunc(.1*nrow(X)),tvalues,TRACEPLOT=TRUE,DEN=TRUE,INT=TRUE,...){
+plot.bayesm.mat=function(x,names,burnin=trunc(.1*nrow(X)),tvalues,TRACEPLOT=TRUE,DEN=TRUE,INT=TRUE,
+ CHECK_NDRAWS=TRUE,...){
#
# S3 method to print matrices of draws the object X is of class "bayesm.mat"
#
@@ -11,7 +12,7 @@ plot.bayesm.mat=function(x,names,burnin=trunc(.1*nrow(X)),tvalues,TRACEPLOT=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(nrow(X) < 100 & CHECK_NDRAWS) {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
diff --git a/R/plot.bayesm.nmix.R b/R/plot.bayesm.nmix.R
index 6c498bd..be4cb72 100755
--- a/R/plot.bayesm.nmix.R
+++ b/R/plot.bayesm.nmix.R
@@ -1,9 +1,13 @@
-plot.bayesm.nmix=function(x,names,burnin=trunc(.1*nrow(probdraw)),Grid,bi.sel,nstd=2,...){
+plot.bayesm.nmix=function(x,names,burnin=trunc(.1*nrow(probdraw)),Grid,bi.sel,nstd=2,marg=TRUE,
+ Data,ngrid=50,ndraw=200,...){
#
# 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
+# P. Rossi 2/07
+# P. Rossi 3/07 fixed problem with dropping dimensions on probdraw (if ncomp=1)
+# P. Rossi 2/08 added marg flag to plot marginals
+# P. Rossi 3/08 added Data argument to paint histograms on the marginal plots
#
nmixlist=x
if(mode(nmixlist) != "list") stop(" Argument must be a list \n")
@@ -15,53 +19,71 @@ plot.bayesm.nmix=function(x,names,burnin=trunc(.1*nrow(probdraw)),Grid,bi.sel,ns
R=nrow(probdraw)
if(R < 100) {cat(" fewer than 100 draws submitted \n"); return(invisible())}
datad=length(compdraw[[1]][[1]]$mu)
+ OneDimData=(datad==1)
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))
+ ind=as.integer(seq(from=(burnin+1),to=R,length.out=max(ndraw,trunc(.05*R))))
+ if(missing(names)) {names=as.character(1:datad)}
+ if(!missing(Data)){
+ if(!is.matrix(Data)) stop("Data argument must be a matrix \n")
+ if(ncol(Data)!= datad) stop("Data matrix is of wrong dimension \n")
}
+
+ if(missing(Grid)){
+ Grid=matrix(0,nrow=ngrid,ncol=datad)
+ if(!missing(Data))
+ {for(i in 1:datad) Grid[,i]=c(seq(from=range(Data[,i])[1],to=range(Data[,i])[2],length=ngrid))}
+ else
+ {out=momMix(probdraw[ind,,drop=FALSE],compdraw[ind])
+ mu=out$mu
+ sd=out$sd
+ for(i in 1:datad ) Grid[,i]=c(seq(from=(mu[i]-nstd*sd[i]),
+ to=(mu[i]+nstd*sd[i]),length=ngrid))
+ }
+ }
#
# 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(marg){
+ mden=eMixMargDen(Grid,probdraw[ind,,drop=FALSE],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){
+ 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")
+ if(!missing(Data)){
+ deltax=(range(Grid[,index])[2]-range(Grid[,index])[1])/nrow(Grid)
+ hist(Data[,index],xlim=range(Grid[,index]),
+ freq=FALSE,col="yellow",breaks=max(20,.1*nrow(Data)),add=TRUE)
+ lines(Grid[,index],mden[,index]/(sum(mden[,index])*deltax),col="red",lwd=2)}
+ else
+ {lines(Grid[,index],mden[,index],col="black",lwd=2)
+ 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
#
+ if(!OneDimData){
par(ask=dev.interactive())
nsel=length(bi.sel)
- ngrid=50
- den=array(0,dim=c(nsel,ngrid,ngrid))
+ den=array(0,dim=c(ngrid,ngrid,nsel))
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)
+ xi=Grid[,i]
+ xj=Grid[,j]
lstxixj[[sel]]=list(xi,xj)
for(elt in ind){
- den[sel,,]=den[sel,,]+mixDenBi(i,j,xi,xj,probdraw[elt,],compdraw[[elt]])
+ den[,,sel]=den[,,sel]+mixDenBi(i,j,xi,xj,probdraw[elt,,drop=FALSE],compdraw[[elt]])
}
+ den[,,sel]=den[,,sel]/sum(den[,,sel])
}
nx=nsel
par(mfrow=c(1,1))
@@ -69,10 +91,11 @@ plot.bayesm.nmix=function(x,names,burnin=trunc(.1*nrow(probdraw)),Grid,bi.sel,ns
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)
+ xlabtxt=names[bi.sel[[index]][1]]
+ ylabtxt=names[bi.sel[[index]][2]]
+ 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/rDPGibbs.R b/R/rDPGibbs.R
new file mode 100755
index 0000000..9367537
--- /dev/null
+++ b/R/rDPGibbs.R
@@ -0,0 +1,569 @@
+rDPGibbs=
+function(Prior,Data,Mcmc)
+{
+#
+# Revision History:
+# 5/06 add rthetaDP
+# 7/06 include rthetaDP in main body to avoid copy overhead
+# 1/08 add scaling
+# 2/08 add draw of lambda
+# 3/08 changed nu prior support to dim(y) + exp(unif gird on nulim[1],nulim[2])
+#
+# purpose: do Gibbs sampling for density estimation using Dirichlet process model
+#
+# arguments:
+# Data is a list of y which is an n x k matrix of data
+# Prior is a list of (alpha,lambda,Prioralpha)
+# alpha: starting value
+# lambda_hyper: hyperparms of prior on lambda
+# Prioralpha: hyperparms of alpha prior; a list of (Istarmin,Istarmax,power)
+# if elements of the prior don't exist, defaults are assumed
+# Mcmc is a list of (R,keep,maxuniq)
+# R: number of draws
+# keep: thinning parameter
+# maxuniq: the maximum number of unique thetaStar values
+#
+# Output:
+# list with elements
+# alphadraw: vector of length R/keep, [i] is ith draw of alpha
+# Istardraw: vector of length R/keep, [i] is the number of unique theta's drawn from ith iteration
+# adraw
+# nudraw
+# vdraw
+# thetaNp1draws: list, [[i]] is ith draw of theta_{n+1}
+# inddraw: R x n matrix, [,i] is indicators of identity for each obs in ith iteration
+#
+# Model:
+# y_i ~ f(y|thetai)
+# thetai|G ~ G
+# G|lambda,alpha ~ DP(G|G0(lambda),alpha)
+#
+# Priors:
+# alpha: starting value
+#
+# lambda:
+# G0 ~ N(mubar,Sigma (x) Amu^-1)
+# mubar=vec(mubar)
+# Sigma ~ IW(nu,nu*v*I) note: mode(Sigma)=nu/(nu+2)*v*I
+# mubar=0
+# amu is uniform on grid specified by alim
+# nu is log uniform, nu=d-1+exp(Z) z is uniform on seq defined bvy nulim
+# v is uniform on sequence specificd by vlim
+#
+# Prioralpha:
+# alpha ~ (1-(alpha-alphamin)/(alphamax-alphamin))^power
+# alphamin=exp(digamma(Istarmin)-log(gamma+log(N)))
+# alphamax=exp(digamma(Istarmax)-log(gamma+log(N)))
+# gamma= .5772156649015328606
+#
+#
+#
+# define needed functions
+#
+# -----------------------------------------------------------------------------------------
+#
+q0=function(y,lambda,eta){
+#
+# function to compute a vector of int f(y[i]|theta) p(theta|lambda)dlambda
+# here p(theta|lambda) is G0 the base prior
+#
+# implemented for a multivariate normal data density and standard conjugate
+# prior:
+# theta=list(mu,Sigma)
+# f(y|theta,eta) is N(mu,Sigma)
+# lambda=list(mubar,Amu,nu,V)
+# mu|Sigma ~ N(mubar,Sigma (x) Amu^-1)
+# Sigma ~ IW(nu,V)
+#
+# arguments:
+# Y is n x k matrix of observations
+# lambda=list(mubar,Amu,nu,V)
+# eta is not used
+#
+# output:
+# vector of q0 values for each obs (row of Y)
+#
+# p. rossi 12/05
+#
+# here y is matrix of observations (each row is an obs)
+
+mubar=lambda$mubar; nu=lambda$nu ; Amu=lambda$Amu; V=lambda$V
+k=ncol(y)
+R=chol(V)
+logdetR=sum(log(diag(R)))
+if (k > 1)
+ {lnk1k2=(k/2)*log(2)+log((nu-k)/2)+lgamma((nu-k)/2)-lgamma(nu/2)+sum(log(nu/2-(1:(k-1))/2))}
+else
+ {lnk1k2=(k/2)*log(2)+log((nu-k)/2)+lgamma((nu-k)/2)-lgamma(nu/2)}
+constant=-(k/2)*log(2*pi)+(k/2)*log(Amu/(1+Amu)) + lnk1k2 + nu*logdetR
+#
+# note: here we are using the fact that |V + S_i | = |R|^2 (1 + v_i'v_i)
+# where v_i = sqrt(Amu/(1+Amu))*t(R^-1)*(y_i-mubar), R is chol(V)
+#
+# and S_i = Amu/(1+Amu) * (y_i-mubar)(y_i-mubar)'
+#
+mat=sqrt(Amu/(1+Amu))*t(backsolve(R,diag(ncol(y))))%*%(t(y)-mubar)
+vivi=colSums(mat^2)
+
+lnq0v=constant-((nu+1)/2)*(2*logdetR+log(1+vivi))
+
+return(exp(lnq0v))
+}
+# ----------------------------------------------------------------------------------------------
+ rmultinomF=
+ function(p) {
+ return(sum(runif(1) > cumsum(p))+1)
+ }
+# -----------------------------------------------------------------------------------------------
+ alphaD=function(Prioralpha,Istar,gridsize){
+#
+# function to draw alpha using prior, p(alpha)= (1-(alpha-alphamin)/(alphamax-alphamin))**power
+#
+ power=Prioralpha$power
+ alphamin=Prioralpha$alphamin
+ alphamax=Prioralpha$alphamax
+ n=Prioralpha$n
+ alpha=seq(from=alphamin,to=(alphamax-0.000001),len=gridsize)
+ lnprob=Istar*log(alpha) + lgamma(alpha) - lgamma(n+alpha) +
+ power*log(1-(alpha-alphamin)/(alphamax-alphamin))
+ lnprob=lnprob-median(lnprob)
+ probs=exp(lnprob)
+ probs=probs/sum(probs)
+ return(alpha[rmultinomF(probs)])
+}
+
+
+#
+# ------------------------------------------------------------------------------------------
+#
+yden=function(thetaStar,y,eta){
+#
+# function to compute f(y | theta)
+# computes f for all values of theta in theta list of lists
+#
+# arguments:
+# thetaStar is a list of lists. thetaStar[[i]] is a list with components, mu, rooti
+# y |theta[[i]] ~ N(mu,(rooti %*% t(rooti))^-1) rooti is inverse of Chol root of Sigma
+# eta is not used
+#
+# output:
+# length(thetaStar) x n array of values of f(y[j,]|thetaStar[[i]]
+#
+
+nunique=length(thetaStar)
+n=nrow(y)
+ydenmat=matrix(double(n*nunique),ncol=n)
+k=ncol(y)
+for(i in 1:nunique){
+
+ # now compute vectorized version of lndMvn
+ # compute y_i'RIRI'y_i for all i
+ #
+ mu=thetaStar[[i]]$mu; rooti=thetaStar[[i]]$rooti
+ quads=colSums((crossprod(rooti,(t(y)-mu)))^2)
+ ydenmat[i,]=exp(-(k/2)*log(2*pi) + sum(log(diag(rooti))) - .5*quads)
+
+}
+return(ydenmat)
+}
+
+#
+# -----------------------------------------------------------------------------------------
+#
+GD=function(lambda){
+#
+# function to draw from prior for Multivariate Normal Model
+#
+# mu|Sigma ~ N(mubar,Sigma x Amu^-1)
+# Sigma ~ IW(nu,V)
+#
+# note: we must insure that mu is a vector to use most efficient
+# lndMvn routine
+#
+nu=lambda$nu
+V=lambda$V
+mubar=lambda$mubar
+Amu=lambda$Amu
+k=length(mubar)
+Sigma=rwishart(nu,chol2inv(chol(lambda$V)))$IW
+root=chol(Sigma)
+mu=mubar+(1/sqrt(Amu))*t(root)%*%matrix(rnorm(k),ncol=1)
+return(list(mu=as.vector(mu),rooti=backsolve(root,diag(k))))
+}
+
+#
+# -------------------------------------------------------------------------------------------
+#
+thetaD=function(y,lambda,eta){
+#
+# function to draw from posterior of theta given data y and base prior G0(lambda)
+#
+# here y ~ N(mu,Sigma)
+# theta = list(mu=mu,rooti=chol(Sigma)^-1)
+# mu|Sigma ~ N(mubar,Sigma (x) Amu-1)
+# Sigma ~ IW(nu,V)
+#
+# arguments:
+# y is n x k matrix of obs
+# lambda is list(mubar,Amu,nu,V)
+# eta is not used
+# output:
+# one draw of theta, list(mu,rooti)
+# Sigma=inv(rooti)%*%t(inv(rooti))
+#
+# note: we assume that y is a matrix. if there is only one obs, y is a 1 x k matrix
+#
+rout=rmultireg(y,matrix(c(rep(1,nrow(y))),ncol=1),matrix(lambda$mubar,nrow=1),matrix(lambda$Amu,ncol=1),
+ lambda$nu,lambda$V)
+return(list(mu=as.vector(rout$B),rooti=backsolve(chol(rout$Sigma),diag(ncol(y)))))
+}
+
+#
+# --------------------------------------------------------------------------------------------
+# load a faster version of lndMvn
+# note: version of lndMvn below assumes x,mu is a vector!
+lndMvn=function (x, mu, rooti)
+{
+ return(-(length(x)/2) * log(2 * pi) - 0.5 * sum(((x-mu)%*%rooti)**2) + sum(log(diag(rooti))))
+}
+# -----------------------------------------------------------------------------------------
+ lambdaD=function(lambda,thetastar,alim=c(.01,2),nulim=c(.01,2),vlim=c(.1,5),gridsize=20){
+#
+# revision history
+# p. rossi 7/06
+# vectorized 1/07
+# changed 2/08 to paramaterize V matrix of IW prior to nu*v*I; then mode of Sigma=nu/(nu+2)vI
+# this means that we have a reparameterization to v* = nu*v
+#
+# function to draw (nu, v, a) using uniform priors
+#
+# theta_j=(mu_j,Sigma_j) mu_j~N(0,Sigma_j/a) Sigma_j~IW(nu,vI)
+# recall E[Sigma]= vI/(nu-dim-1)
+#
+# define functions needed
+# ----------------------------------------------------------------------------------------------
+ rmultinomF=
+ function(p) {
+ return(sum(runif(1) > cumsum(p))+1)
+ }
+echo=function(lst){return(t(lst[[2]]))}
+rootiz=function(lst){crossprod(lst[[2]],lst[[1]])}
+#
+# ------------------------------------------------------------------------------------------
+
+ d=length(thetastar[[1]]$mu)
+ Istar=length(thetastar)
+ aseq=seq(from=alim[1],to=alim[2],len=gridsize)
+ nuseq=d-1+exp(seq(from=nulim[1],to=nulim[2],len=gridsize)) # log uniform grid
+ vseq=seq(from=vlim[1],to=vlim[2],len=gridsize)
+#
+# "brute" force approach would simply loop over the
+# "observations" (theta_j) and use log of the appropriate densities. To vectorize, we
+# notice that the "data" comes via various statistics:
+# 1. sum of log(diag(rooti_j)
+# 2. sum of tr(V%*%rooti_j%*%t(rooti_j)) where V=vI_d
+# 3. quadratic form t(mu_j-0)%*%rooti%*%t(rooti)%*%(mu_j-0)
+# thus, we will compute these first.
+# for documentation purposes, we leave brute force code in comment fields
+#
+# extract needed info from thetastar list
+#
+ out=double(Istar*d*d)
+ out=sapply(thetastar,echo)
+ dim(out)=c(d,Istar*d) # out has the rootis in form: [t(rooti_1), t(rooti_2), ...,t(rooti_Istar)]
+ sumdiagriri=sum(colSums(out^2)) # sum_j tr(rooti_j%*%t(rooti_j))
+# now get diagonals of rooti
+ ind=cbind(c(1:(d*Istar)),rep((1:d),Istar))
+ out=t(out)
+ sumlogdiag=sum(log(out[ind]))
+ rimu=sapply(thetastar,rootiz) # columns of rimu contain t(rooti_j)%*%mu_j
+ dim(rimu)=c(d,Istar)
+ sumquads=sum(colSums(rimu^2))
+#
+# draw a (conditionally indep of nu,v given theta_j)
+ lnprob=double(length(aseq))
+ #for(i in seq(along=aseq)){
+ #for(j in seq(along=thetastar)){
+ #lnprob[i]=lnprob[i]+lndMvn(thetastar[[j]]$mu,c(rep(0,d)),thetastar[[j]]$rooti*sqrt(aseq[i]))}
+ lnprob=Istar*(-(d/2)*log(2*pi))-.5*aseq*sumquads+Istar*d*log(sqrt(aseq))+sumlogdiag
+ lnprob=lnprob-max(lnprob)+200
+ probs=exp(lnprob)
+ probs=probs/sum(probs)
+ adraw=aseq[rmultinomF(probs)]
+#
+# draw nu given v
+#
+ V=lambda$V
+ lnprob=double(length(nuseq))
+ #for(i in seq(along=nuseq)){
+ #for(j in seq(along=thetastar)){
+ #Sigma_j=crossprod(backsolve(thetastar[[j]]$rooti,diag(d)))
+ #lnprob[i]=lnprob[i]+lndIWishart(nuseq[i],V,Sigma_j)}
+ arg=rep(c(1:d),gridsize)
+ dim(arg)=c(d,gridsize)
+ arg=t(arg)
+ arg=(nuseq+1-arg)/2
+ lnprob=-Istar*log(2)*d/2*nuseq - Istar*rowSums(lgamma(arg)) +
+ Istar*d*log(sqrt(V[1,1]))*nuseq + sumlogdiag*nuseq
+ lnprob=lnprob-max(lnprob)+200
+ probs=exp(lnprob)
+ probs=probs/sum(probs)
+ nudraw=nuseq[rmultinomF(probs)]
+#
+# draw v given nu
+#
+ lnprob=double(length(vseq))
+ #for(i in seq(along=vseq)){
+ #V=vseq[i]*diag(d)
+ #for(j in seq(along=thetastar)){
+ #Sigma_j=crossprod(backsolve(thetastar[[j]]$rooti,diag(d)))
+ #lnprob[i]=lnprob[i]+lndIWishart(nudraw,V,Sigma_j)}
+# lnprob=Istar*nudraw*d*log(sqrt(vseq))-.5*sumdiagriri*vseq
+ lnprob=Istar*nudraw*d*log(sqrt(vseq*nudraw))-.5*sumdiagriri*vseq*nudraw
+ lnprob=lnprob-max(lnprob)+200
+ probs=exp(lnprob)
+ probs=probs/sum(probs)
+ vdraw=vseq[rmultinomF(probs)]
+#
+# put back into lambda
+#
+ return(list(mubar=c(rep(0,d)),Amu=adraw,nu=nudraw,V=nudraw*vdraw*diag(d)))
+}
+# -----------------------------------------------------------------------------------------
+
+# check arguments
+#
+pandterm=function(message) {stop(message,call.=FALSE)}
+if(missing(Data)) {pandterm("Requires Data argument -- list of y")}
+ if(is.null(Data$y)) {pandterm("Requires Data element y")}
+ y=Data$y
+#
+# check data for validity
+#
+if(!is.matrix(y)) {pandterm("y must be a matrix")}
+nobs=nrow(y)
+dimy=ncol(y)
+#
+# check for Prior
+#
+alimdef=c(.01,10)
+nulimdef=c(.01,3)
+vlimdef=c(.1,4)
+if(missing(Prior)) {pandterm("requires Prior argument ")}
+else
+ {
+ if(is.null(Prior$lambda_hyper)) {lambda_hyper=list(alim=alimdef,nulim=nulimdef,vlim=vlimdef)}
+ else {lambda_hyper=Prior$lambda_hyper;
+ if(is.null(lambda_hyper$alim)) {lambda_hyper$alim=alimdef}
+ if(is.null(lambda_hyper$nulim)) {lambda_hyper$nulim=nulimdef}
+ if(is.null(lambda_hyper$vlim)) {lambda_hyper$vlim=vlimdef}
+ }
+ if(is.null(Prior$Prioralpha)) {Prioralpha=list(Istarmin=1,Istarmax=min(50,0.1*nobs),power=0.8)}
+ else {Prioralpha=Prior$Prioralpha;
+ if(is.null(Prioralpha$Istarmin)) {Prioralpha$Istarmin=1} else {Prioralpha$Istarmin=Prioralpha$Istarmin}
+ if(is.null(Prioralpha$Istarmax))
+ {Prioralpha$Istarmax=min(50,0.1*nobs)} else {Prioralpha$Istarmax=Prioralpha$Istarmax}
+ if(is.null(Prioralpha$power)) {Prioralpha$power=0.8}
+ }
+ }
+gamma= .5772156649015328606
+Prioralpha$alphamin=exp(digamma(Prioralpha$Istarmin)-log(gamma+log(nobs)))
+Prioralpha$alphamax=exp(digamma(Prioralpha$Istarmax)-log(gamma+log(nobs)))
+Prioralpha$n=nobs
+#
+# check Prior arguments for valdity
+#
+if(lambda_hyper$alim[1]<0) {pandterm("alim[1] must be >0")}
+if(lambda_hyper$nulim[1]<0) {pandterm("nulim[1] must be >0")}
+if(lambda_hyper$vlim[1]<0) {pandterm("vlim[1] must be >0")}
+if(Prioralpha$Istarmin <1){pandterm("Prioralpha$Istarmin must be >= 1")}
+if(Prioralpha$Istarmax <= Prioralpha$Istarmin){pandterm("Prioralpha$Istarmin must be > Prioralpha$Istarmax")}
+#
+# 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$maxuniq)) {maxuniq=200} else {maxuniq=Mcmc$maxuniq}
+ if(is.null(Mcmc$SCALE)) {SCALE=TRUE} else {SCALE=Mcmc$SCALE}
+ if(is.null(Mcmc$gridsize)) {gridsize=20} else {gridsize=Mcmc$gridsize}
+ }
+
+#
+# print out the problem
+#
+cat(" Starting Gibbs Sampler for Density Estimation Using Dirichlet Process Model",fill=TRUE)
+cat(" ",nobs," observations on ",dimy," dimensional data",fill=TRUE)
+cat(" ",fill=TRUE)
+cat(" SCALE=",SCALE,fill=TRUE)
+cat(" ",fill=TRUE)
+cat(" Prior Parms: ",fill=TRUE)
+cat(" G0 ~ N(mubar,Sigma (x) Amu^-1)",fill=TRUE)
+cat(" mubar = ",0,fill=TRUE)
+cat(" Sigma ~ IW(nu,nu*v*I)",fill=TRUE)
+cat(" Amu ~ uniform[",lambda_hyper$alim[1],",",lambda_hyper$alim[2],"]",fill=TRUE)
+cat(" nu ~ uniform on log grid on [",dimy-1+exp(lambda_hyper$nulim[1]),
+ ",",dimy-1+exp(lambda_hyper$nulim[2]),"]",fill=TRUE)
+cat(" v ~ uniform[",lambda_hyper$vlim[1],",",lambda_hyper$vlim[2],"]",fill=TRUE)
+cat(" ",fill=TRUE)
+cat(" alpha ~ (1-(alpha-alphamin)/(alphamax-alphamin))^power",fill=TRUE)
+cat(" Istarmin = ",Prioralpha$Istarmin,fill=TRUE)
+cat(" Istarmax = ",Prioralpha$Istarmax,fill=TRUE)
+cat(" alphamin = ",Prioralpha$alphamin,fill=TRUE)
+cat(" alphamax = ",Prioralpha$alphamax,fill=TRUE)
+cat(" power = ",Prioralpha$power,fill=TRUE)
+cat(" ",fill=TRUE)
+cat(" Mcmc Parms: R= ",R," keep= ",keep," maxuniq= ",maxuniq," gridsize for lambda hyperparms= ",gridsize,
+ fill=TRUE)
+cat(" ",fill=TRUE)
+
+# initialize theta, thetastar, indic
+
+theta=vector("list",nobs)
+for(i in 1:nobs) {theta[[i]]=list(mu=rep(0,dimy),rooti=diag(dimy))}
+indic=double(nobs)
+thetaStar=unique(theta)
+nunique=length(thetaStar)
+for(j in 1:nunique){
+ indic[which(sapply(theta,identical,thetaStar[[j]]))]=j
+}
+#
+# initialize lambda
+#
+lambda=list(mubar=rep(0,dimy),Amu=.1,nu=dimy+1,V=(dimy+1)*diag(dimy))
+
+#
+# initialize alpha
+#
+alpha=1
+
+alphadraw=double(floor(R/keep))
+Istardraw=double(floor(R/keep))
+adraw=double(floor(R/keep))
+nudraw=double(floor(R/keep))
+vdraw=double(floor(R/keep))
+thetaNp1draw=vector("list",floor(R/keep))
+inddraw=matrix(double((floor(R/keep))*nobs),ncol=nobs)
+
+#
+# do scaling
+#
+if(SCALE){
+ dvec=sqrt(apply(y,2,var))
+ ybar=apply(y,2,mean)
+ y=scale(y,center=ybar,scale=dvec)
+ dvec=1/dvec # R function scale divides by scale
+}
+#
+# note on scaling
+#
+# we model scaled y, z_i=D(y_i-ybar) D=diag(1/sigma1, ..., 1/sigma_dimy)
+#
+# if p_z= 1/R sum(phi(z|mu,Sigma))
+# p_y=1/R sum(phi(y|D^-1mu+ybar,D^-1SigmaD^-1)
+# rooti_y=Drooti_z
+#
+# you might want to use quantiles instead, like median and (10,90)
+#
+
+#
+# start main iteration loop
+#
+itime=proc.time()[3]
+cat("MCMC Iteration (est time to end -min) ",fill=TRUE)
+fsh()
+
+for(rep in 1:R)
+{
+ n = length(theta)
+
+ eta=NULL # note eta is not used
+ thetaNp1=NULL
+ q0v = q0(y,lambda,eta) # now that we draw lambda we need to recompute q0v each time
+
+ p=c(rep(1/(alpha+(n-1)),n-1),alpha/(alpha+(n-1)))
+
+ nunique=length(thetaStar)
+
+ if(nunique > maxuniq ) { pandterm("maximum number of unique thetas exceeded")}
+ ydenmat=matrix(double(maxuniq*n),ncol=n)
+ ydenmat[1:nunique,]=yden(thetaStar,y,eta)
+ # ydenmat is a length(thetaStar) x n array of density values given f(y[j,] | thetaStar[[i]]
+ # note: due to remix step (below) we must recompute ydenmat each time!
+
+ # use .Call to draw theta list
+ out= .Call("thetadraw",y,ydenmat,indic,q0v,p,theta,lambda,eta=eta,
+ thetaD=thetaD,yden=yden,maxuniq,nunique,new.env())
+
+ # theta has been modified by thetadraw so we need to recreate thetaStar
+ thetaStar=unique(theta)
+ nunique=length(thetaStar)
+
+ #thetaNp1 and remix
+ probs=double(nunique+1)
+ for(j in 1:nunique) {
+ ind = which(sapply(theta,identical,thetaStar[[j]]))
+ probs[j]=length(ind)/(alpha+n)
+ new_utheta=thetaD(y[ind,,drop=FALSE],lambda,eta)
+ for(i in seq(along=ind)) {theta[[ind[i]]]=new_utheta}
+ indic[ind]=j
+ thetaStar[[j]]=new_utheta
+ }
+ probs[nunique+1]=alpha/(alpha+n)
+ ind=rmultinomF(probs)
+ if(ind==length(probs)) {
+ thetaNp1=GD(lambda)
+ } else {
+ thetaNp1=thetaStar[[ind]]
+ }
+
+ # draw alpha
+ alpha=alphaD(Prioralpha,nunique,gridsize=gridsize)
+
+ # draw lambda
+ lambda=lambdaD(lambda,thetaStar,alim=lambda_hyper$alim,nulim=lambda_hyper$nulim,
+ vlim=lambda_hyper$vlim,gridsize=gridsize)
+
+ 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
+ alphadraw[mkeep]=alpha
+ Istardraw[mkeep]=nunique
+ adraw[mkeep]=lambda$Amu
+ nudraw[mkeep]=lambda$nu
+ vdraw[mkeep]=lambda$V[1,1]/lambda$nu
+ if(SCALE){
+ thetaNp1[[1]]=thetaNp1[[1]]/dvec+ybar
+ if(ncol(y)>1)
+ {thetaNp1[[2]]=diag(dvec)%*%thetaNp1[[2]]}
+ else
+ {thetaNp1[[2]]=dvec*thetaNp1[[2]]}
+ }
+ thetaNp1draw[[mkeep]]=list(list(mu=thetaNp1[[1]],rooti=thetaNp1[[2]]))
+ # here we put the draws into the list of lists of list format useful for
+ # finite mixture of normals utilities
+ inddraw[mkeep,]=indic
+ }
+}
+ctime=proc.time()[3]
+cat("Total Time Elapsed= ",round((ctime-itime)/60,2),fill=TRUE)
+nmix=list(probdraw=matrix(c(rep(1,nrow(inddraw))),ncol=1),zdraw=inddraw,compdraw=thetaNp1draw)
+attributes(nmix)$class="bayesm.nmix"
+attributes(alphadraw)$class=c("bayesm.mat","mcmc")
+attributes(Istardraw)$class=c("bayesm.mat","mcmc")
+attributes(adraw)$class=c("bayesm.mat","mcmc")
+attributes(nudraw)$class=c("bayesm.mat","mcmc")
+attributes(vdraw)$class=c("bayesm.mat","mcmc")
+return(list(alphadraw=alphadraw,Istardraw=Istardraw,adraw=adraw,nudraw=nudraw,
+ vdraw=vdraw,nmix=nmix))
+}
diff --git a/R/rhierLinearMixture.R b/R/rhierLinearMixture.R
index e4acf77..e024a19 100755
--- a/R/rhierLinearMixture.R
+++ b/R/rhierLinearMixture.R
@@ -57,7 +57,10 @@ function(Data,Prior,Mcmc)
#
append=function(l) { l=c(l,list(XpX=crossprod(l$X),Xpy=crossprod(l$X,l$y)))}
#
-getvar=function(l) { var(l$y)}
+getvar=function(l) {
+ v=var(l$y)
+ if(is.na(v)) return(1)
+ if(v>0) return (v) else return (1)}
#
runiregG=
function(y,X,XpX,Xpy,sigmasq,rooti,betabar,nu,ssq){
diff --git a/R/rhierLinearModel.R b/R/rhierLinearModel.R
index 23d12b0..fb75af2 100755
--- a/R/rhierLinearModel.R
+++ b/R/rhierLinearModel.R
@@ -56,7 +56,10 @@ function(Data,Prior,Mcmc)
#------------------------------------------------------------------------------
append=function(l) { l=c(l,list(XpX=crossprod(l$X),Xpy=crossprod(l$X,l$y)))}
#
-getvar=function(l) { var(l$y)}
+getvar=function(l) {
+ v=var(l$y)
+ if(is.na(v)) return(1)
+ if(v>0) return (v) else return (1)}
#
runiregG=
function(y,X,XpX,Xpy,sigmasq,A,betabar,nu,ssq){
@@ -105,7 +108,7 @@ if(missing(Data)) {pandterm("Requires Data argument -- list of regdata and Z")}
regdata=Data$regdata
nreg=length(regdata)
if(is.null(Data$Z)) { cat("Z not specified -- putting in iota",fill=TRUE); fsh() ; Z=matrix(rep(1,nreg),ncol=1)}
- else {if (nrow(Data$Z) != nreg) {pandterm(paste("Nrow(Z) ",nrow(Z),"ne number regressions ",nlgt))}
+ else {if (nrow(Data$Z) != nreg) {pandterm(paste("Nrow(Z) ",nrow(Z),"ne number regressions ",nreg))}
else {Z=Data$Z}}
nz=ncol(Z)
#
diff --git a/R/rhierNegbinRw.R b/R/rhierNegbinRw.R
index 3761ca9..f9c5e6a 100755
--- a/R/rhierNegbinRw.R
+++ b/R/rhierNegbinRw.R
@@ -6,6 +6,7 @@ function(Data, Prior, Mcmc) {
# P. Rossi 6/05
# fixed error with nobs not specified and changed llnegbinFract 9/05
# 3/07 added classes
+# 3/08 fixed fractional likelihood
#
# Model
# (y_i|lambda_i,alpha) ~ Negative Binomial(Mean = lambda_i, Overdispersion par = alpha)
@@ -37,7 +38,7 @@ function(Data, Prior, Mcmc) {
# keep is thinning parameter (def = 1)
# s_beta - scaling parameter for beta RW (def = 2.93/sqrt(nvar))
# s_alpha - scaling parameter for alpha RW (def = 2.93)
-# c - fractional weighting parameter (def = 2)
+# w - fractional weighting parameter (def = .1)
# Vbeta0, Delta0 - initial guesses for parameters, if not supplied default values are used
#
@@ -59,11 +60,10 @@ function(par,X,y, nvar) {
}
llnegbinFract =
-function(par,X,y,Xpooled, ypooled, power, nvar,lnalpha) {
+function(par,X,y,Xpooled, ypooled, w,wgt, nvar,lnalpha) {
# Computes the fractional log-likelihood at the unit level
-# = l_i * l_bar^power
theta = c(par,lnalpha)
- llnegbin(theta,X,y,nvar) + power*llnegbin(theta,Xpooled,ypooled, nvar)
+ (1-w)*llnegbin(theta,X,y,nvar) + w*wgt*llnegbin(theta,Xpooled,ypooled, nvar)
}
lpostbetai =
@@ -165,12 +165,12 @@ if(sum(dim(Vbeta0) == c(nvar,nvar)) !=2) pandterm("Vbeta0 is not of dimension nv
if(is.null(Mcmc$Delta0)) {Delta0=matrix(rep(0,nz*nvar),nrow=nz)} else {Delta0=Mcmc$Delta0}
if(sum(dim(Delta0) == c(nz,nvar)) !=2) pandterm("Delta0 is not of dimension nvar by nz")
if(is.null(Mcmc$keep)) {keep=1} else {keep=Mcmc$keep}
-if(is.null(Mcmc$s_alpha)) {cat("Using default s_alpha = 2.93",fill=TRUE); s_alpha=2.93}
+if(is.null(Mcmc$s_alpha)) { s_alpha=2.93}
else {s_alpha= Mcmc$s_alpha }
-if(is.null(Mcmc$s_beta)) {cat("Using default s_beta = 2.93/sqrt(nvar)",fill=TRUE); s_beta=2.93/sqrt(nvar)}
+if(is.null(Mcmc$s_beta)) { s_beta=2.93/sqrt(nvar)}
else {s_beta=Mcmc$s_beta }
-if(is.null(Mcmc$c)) {cat("Using default c = 2"); c=2}
- else {c = Mcmc$c}
+if(is.null(Mcmc$w)) { w=.1}
+ else {w = Mcmc$w}
#out = rhierNegbinRw(Data, Prior, Mcmc)
# print out problem
@@ -198,7 +198,7 @@ cat("MCMC Parameters:",fill=TRUE)
cat(R," reps; keeping every ",keep,"th draw",fill=TRUE)
cat("s_alpha = ",s_alpha,fill=TRUE)
cat("s_beta = ",s_beta,fill=TRUE)
-cat("Fractional Scaling c = ",c,fill=TRUE)
+cat("Fractional Likelihood Weight Parameter = ",w,fill=TRUE)
cat(" ",fill=TRUE)
par = rep(0,(nvar+1))
@@ -222,16 +222,21 @@ alphacroot = sqrt(alphacvar)
#fsh()
hess_i=NULL
+if(nobs > 1000){
+ sind=sample(c(1:nobs),size=1000)
+ ypooleds=ypooled[sind]
+ Xpooleds=Xpooled[sind,]
+ }
# Find the individual candidate hessian
for (i in 1:nreg) {
- power = length(regdata[[i]]$y)/(c*nobs)
- mle2 = optim(mle$par[1:nvar],llnegbinFract, X=regdata[[i]]$X, y=regdata[[i]]$y, Xpooled=Xpooled,
- ypooled=ypooled, power=power, nvar=nvar, lnalpha=mle$par[nvar+1],
+ wgt = length(regdata[[i]]$y)/length(ypooleds)
+ mle2 = optim(mle$par[1:nvar],llnegbinFract, X=regdata[[i]]$X, y=regdata[[i]]$y, Xpooled=Xpooleds,
+ ypooled=ypooleds, w=w,wgt=wgt, nvar=nvar, lnalpha=mle$par[nvar+1],
method="BFGS", hessian=TRUE, control=list(fnscale=-1, trace=0))
if (mle2$convergence==0)
hess_i[[i]] = list(hess=-mle2$hessian)
else
- hess_i[[i]] = diag(rep(0,nvar))
+ hess_i[[i]] = diag(rep(1,nvar))
if(i%%50 ==0) cat(" completed unit #",i,fill=TRUE)
fsh()
}
diff --git a/R/rivDP.R b/R/rivDP.R
new file mode 100755
index 0000000..0229b74
--- /dev/null
+++ b/R/rivDP.R
@@ -0,0 +1,697 @@
+rivDP =
+function(Data,Prior,Mcmc)
+{
+#
+# revision history:
+# P. Rossi 1/06
+# added draw of alpha 2/06
+# added automatic scaling 2/06
+# removed reqfun 7/07 -- now functions are in rthetaDP
+#
+# purpose:
+# draw from posterior for linear I.V. model with DP process for errors
+#
+# Arguments:
+# Data -- list of z,w,x,y
+# y is vector of obs on lhs var in structural equation
+# x is "endogenous" var in structural eqn
+# w is matrix of obs on "exogenous" vars in the structural eqn
+# z is matrix of obs on instruments
+# Prior -- list of md,Ad,mbg,Abg,mubar,Amu,nuV
+# md is prior mean of delta
+# Ad is prior prec
+# mbg is prior mean vector for beta,gamma
+# Abg is prior prec of same
+# lamda is a list of prior parms for DP draw
+# mubar is prior mean of means for "errors"
+# Amu is scale precision parm for means
+# nu,V parms for IW on Sigma (idential priors for each normal comp
+# alpha prior parm for DP process (weight on base measure)
+# or starting value if there is a prior on alpha (requires element Prioralpha)
+# Prioralpha list of hyperparms for draw of alpha (alphamin,alphamax,power,n)
+#
+# Mcmc -- list of R,keep,starting values for delta,beta,gamma,theta
+# maxuniq is maximum number of unique theta values
+# R is number of draws
+# keep is thinning parameter
+# SCALE if scale data, def: TRUE
+# gridsize is the gridsize parm for alpha draws
+#
+# Output:
+# list of draws of delta,beta,gamma and thetaNp1 which is used for
+# predictive distribution of errors (density estimation)
+#
+# Model:
+#
+# x=z'delta + e1
+# y=beta*x + w'gamma + e2
+# e1,e2 ~ N(theta_i)
+#
+# Priors
+# delta ~ N(md,Ad^-1)
+# vec(beta,gamma) ~ N(mbg,Abg^-1)
+# theta ~ DPP(alpha|lambda)
+#
+#
+# extract data and check dimensios
+#
+pandterm=function(message) {stop(message,call.=FALSE)}
+if(missing(Data)) {pandterm("Requires Data argument -- list of z,w,x,y")}
+ if(is.null(Data$w)) isgamma=FALSE else isgamma=TRUE
+ if(isgamma) w = Data$w #matrix
+ if(is.null(Data$z)) {pandterm("Requires Data element z")}
+ z=Data$z
+ if(is.null(Data$x)) {pandterm("Requires Data element x")}
+ x=as.vector(Data$x)
+ if(is.null(Data$y)) {pandterm("Requires Data element y")}
+ y=as.vector(Data$y)
+
+#
+# check data for validity
+#
+n=length(y)
+if(isgamma)
+ {if(is.matrix(w)) {pandterm("w is not a matrix")}
+ dimg=ncol(w)
+ if(n != nrow(w) ) {pandterm("length(y) ne nrow(w)")}}
+
+if(!is.matrix(z)) {pandterm("z is not a matrix")}
+dimd=ncol(z)
+if(n != length(x) ) {pandterm("length(y) ne length(x)")}
+if(n != nrow(z) ) {pandterm("length(y) ne nrow(z)")}
+
+
+#
+# extract elements corresponding to the prior
+#
+if(missing(Prior))
+ {
+ md=c(rep(0,dimd))
+ Ad=diag(0.01,dimd)
+ if(isgamma) dimbg=1+dimg else dimbg=1
+ mbg=c(rep(0,dimbg))
+ Abg=diag(0.01,dimbg)
+
+
+ gamma= .5772156649015328606
+ Istarmin=1
+ alphamin=exp(digamma(Istarmin)-log(gamma+log(n)))
+ Istarmax=floor(.1*n)
+ alphamax=exp(digamma(Istarmax)-log(gamma+log(n)))
+ power=.8
+ Prioralpha=list(n=n,alphamin=alphamin,alphamax=alphamax,power=power)
+
+ lambda=list(mubar=c(0,0),Amu=.2,nu=3.4,V=1.7*diag(2))
+ }
+
+else
+ {
+ if(is.null(Prior$md)) md=c(rep(0,dimd)) else md=Prior$md
+ if(is.null(Prior$Ad)) Ad=diag(0.01,dimd) else md=Prior$Ad
+ if(isgamma) dimbg=1+dimg else dimbg=1
+ if(is.null(Prior$mbg)) mbg=c(rep(0,dimbg)) else md=Prior$mbg
+ if(is.null(Prior$Abg)) Abg=diag(0.01,dimbg) else md=Prior$Abg
+
+
+ if(!is.null(Prior$Prioralpha))
+ {Prioralpha=Prior$Prioralpha}
+ else
+ {gamma= .5772156649015328606
+ Istarmin=1
+ alphamin=exp(digamma(Istarmin)-log(gamma+log(n)))
+ Istarmax=floor(.1*n)
+ alphamax=exp(digamma(Istarmax)-log(gamma+log(n)))
+ power=.8
+ Prioralpha=list(n=n,alphamin=alphamin,alphamax=alphamax,power=power)}
+
+ if(!is.null(Prior$lambda))
+ {lambda=Prior$lambda}
+ else
+ {lambda=list(mubar=c(0,0),Amu=.2,nu=3.4,V=1.7*diag(2))}
+ }
+
+#
+# obtain starting values for MCMC
+#
+# we draw need inital values of delta, theta and indic
+#
+
+if(missing(Mcmc)) {pandterm("requires Mcmc argument")}
+theta=NULL
+if(!is.null(Mcmc$delta))
+ {delta = Mcmc$delta}
+else
+ {lmxz = lm(x~z,data.frame(x=x,z=z))
+ delta = lmxz$coef[2:(ncol(z)+1)]}
+if(!is.null(Mcmc$theta))
+ {theta=Mcmc$delta }
+else
+ {onecomp=list(mu=c(0,0),rooti=diag(2))
+ theta=vector("list",length(y))
+ for(i in 1:n) {theta[[i]]=onecomp}
+ }
+dimd = length(delta)
+if(is.null(Mcmc$maxuniq))
+ {maxuniq=200}
+else
+ {maxuniq=Mcmc$maxuniq}
+if(is.null(Mcmc$R)) {pandterm("requres Mcmc argument, R")}
+R = Mcmc$R
+if(is.null(Mcmc$keep))
+ {keep=1}
+else
+ {keep=Mcmc$keep}
+if(is.null(Mcmc$gridsize))
+ {gridsize=20}
+else
+ {gridsize=Mcmc$gridsize}
+if(is.null(Mcmc$SCALE))
+ {SCALE=TRUE}
+else
+ {SCALE=Mcmc$SCALE}
+
+
+#
+# scale and center
+#
+if(SCALE){
+ scaley=sqrt(var(y))
+ scalex=sqrt(var(x))
+ meany=mean(y)
+ meanx=mean(x)
+ meanz=apply(z,2,mean)
+ y=(y-meany)/scaley; x=(x-meanx)/scalex
+ z=scale(z,center=TRUE,scale=FALSE)
+ if(isgamma) {meanw=apply(w,2,mean); w=scale(w,center=TRUE,scale=FALSE)}
+}
+
+#
+# print out model
+#
+cat(" ",fill=TRUE)
+cat("Starting Gibbs Sampler for Linear IV Model With DP Process Errors",fill=TRUE)
+cat(" ",fill=TRUE)
+cat(" nobs= ",n,"; ",ncol(z)," instruments",fill=TRUE)
+cat(" ",fill=TRUE)
+cat("Prior Parms: ",fill=TRUE)
+cat("mean of delta ",fill=TRUE)
+print(md)
+cat(" ",fill=TRUE)
+cat("Adelta",fill=TRUE)
+print(Ad)
+cat(" ",fill=TRUE)
+cat("mean of beta/gamma",fill=TRUE)
+print(mbg)
+cat(" ",fill=TRUE)
+cat("Abeta/gamma",fill=TRUE)
+print(Abg)
+cat(" ",fill=TRUE)
+cat("lambda contains: ", fill=TRUE)
+cat("mu Prior Parms:",fill=TRUE)
+cat("mubar= ",lambda$mubar,fill=TRUE)
+cat("Amu= ",lambda$Amu,fill=TRUE)
+cat(" ",fill=TRUE)
+cat("Sigma Prior Parms:",fill=TRUE)
+cat("nu= ",lambda$nu," V=",fill=TRUE)
+print(lambda$V)
+cat(" ",fill=TRUE)
+cat("Parameters of Prior on Dirichlet Process parm (alpha)",fill=TRUE)
+cat("alphamin= ",Prioralpha$alphamin," alphamax= ",Prioralpha$alphamax," power=",
+ Prioralpha$power,fill=TRUE)
+cat("alpha values correspond to Istarmin = ",Istarmin," Istarmax = ",Istarmax,fill=TRUE)
+cat(" ",fill=TRUE)
+cat("MCMC parms: R= ",R," keep= ",keep,fill=TRUE)
+cat(" maximum number of unique thetas= ",maxuniq,fill=TRUE)
+cat(" gridsize for alpha draws= ",gridsize,fill=TRUE)
+cat(" SCALE data= ",SCALE,fill=TRUE)
+cat(" ",fill=TRUE)
+
+
+#
+# define needed functions
+#
+#
+#
+# --------------------------------------------------------------------------------------------
+#
+#
+get_ytxt=function(y,z,delta,x,w,ncomp,indic,comps){
+yt=NULL; xt=NULL;
+if(missing(w)) isw=FALSE else isw=TRUE
+if(isw) ncolw=ncol(w)
+for (k in 1:ncomp)
+{
+ nobs=sum(indic==k)
+ if(nobs > 0)
+ {
+ if(isw) wk=matrix(w[indic==k,],ncol=ncolw)
+ zk=matrix(z[indic==k,],ncol=length(delta))
+ yk=y[indic==k]
+ xk=matrix(x[indic==k],ncol=1)
+ Sigma=backsolve(comps[[k]][[2]],diag(2))
+ Sigma=crossprod(Sigma)
+ mu=comps[[k]][[1]]
+ e1 = as.vector(xk-zk%*%delta)
+ ee2 = mu[2] +(Sigma[1,2]/Sigma[1,1])*(e1-mu[1])
+ sig = sqrt(Sigma[2,2]-(Sigma[1,2]^2/Sigma[1,1]))
+ yt = c(yt,(yk-ee2)/sig)
+ if(isw)
+ {xt = rbind(xt,(cbind(xk,wk)/sig))}
+ else
+ {xt=rbind(xt,xk/sig)}
+ }
+}
+return(list(xt=xt,yt=yt))
+}
+#
+#
+# --------------------------------------------------------------------------------------------
+#
+#
+get_ytxtd=function(y,z,beta,gamma,x,w,ncomp,indic,comps,dimd){
+yt=NULL; xtd=NULL;
+if(missing(w)) isw=FALSE else isw=TRUE
+if(isw) ncolw=ncol(w)
+C = matrix(c(1,beta,0,1),nrow=2)
+for (k in 1:ncomp)
+ {
+ nobs=sum(indic==k)
+ if(nobs > 0)
+ {
+ xtdk=matrix(nrow=2*nobs,ncol=dimd)
+ ind=seq(1,(2*nobs-1),by=2)
+ if(isw) wk=matrix(w[indic==k,],ncol=ncolw)
+ zk=matrix(z[indic==k,],ncol=dimd)
+ zveck=as.vector(t(zk))
+ yk=y[indic==k]
+ xk=x[indic==k]
+ Sigma=backsolve(comps[[k]][[2]],diag(2))
+ Sigma=crossprod(Sigma)
+ mu=comps[[k]][[1]]
+ B = C%*%Sigma%*%t(C)
+ L = t(chol(B))
+ Li=backsolve(L,diag(2),upper.tri=FALSE)
+ if(isw) {u=as.vector((yk-wk%*%gamma-mu[2]-beta*mu[1]))}
+ else {u=as.vector((yk-mu[2]-beta*mu[1]))}
+ ytk = as.vector(Li %*% rbind((xk-mu[1]),u))
+
+ z2=rbind(zveck,beta*zveck)
+ z2=Li%*%z2
+ zt1=z2[1,]
+ zt2=z2[2,]
+
+ dim(zt1)=c(dimd,nobs)
+ zt1=t(zt1)
+ dim(zt2)=c(dimd,nobs)
+ zt2=t(zt2)
+
+ xtdk[ind,]=zt1
+ xtdk[-ind,]=zt2
+
+ yt=c(yt,ytk)
+ xtd=rbind(xtd,xtdk)
+ }
+ }
+return(list(yt=yt,xtd=xtd))
+}
+#
+#
+# --------------------------------------------------------------------------------------------
+#
+#
+rthetaDP= function(maxuniq,alpha,lambda,Prioralpha,theta,thetaStar,indic,q0v,y,gridsize){
+#
+# function to make one draw from DP process
+#
+# P. Rossi 1/06
+# added draw of alpha 2/06
+# removed lambdaD,etaD and function arguments 5/06
+# removed thetaStar argument to .Call and creation of newthetaStar 7/06
+# removed q0 computations as eta is not drawn 7/06
+# changed for new version of thetadraw and removed calculation of thetaStar before
+# .Call 7/07
+#
+# y(i) ~ f(y|theta[[i]],eta)
+# theta ~ DP(alpha,G(lambda))
+# note: eta is not used
+#output:
+# list with components:
+# thetaDraws: list, [[i]] is a list of the ith draw of the n theta's
+# where n is the length of the input theta and nrow(y)
+# thetaNp1Draws: list, [[i]] is ith draw of theta_{n+1}
+#args:
+# maxuniq: the maximum number of unique thetaStar values -- an error will be raised
+# if this is exceeded
+# alpha,lambda: starting values (or fixed DP prior values if not drawn).
+# Prioralpha: list of hyperparms of alpha prior
+# theta: list of starting value for theta's
+# thetaStar: list of unique values of theta, thetaStar[[i]]
+# indic: n vector of indicator for which unique theta (in thetaStar)
+# y: is a matrix nxk
+# thetaStar: list of unique values of theta, thetaStar[[i]]
+# q0v:a double vector with the same number of rows as y, giving \Int f(y(i)|theta,eta) dG_{lambda}(theta).
+#
+# define needed functions for rthetaDP
+# -----------------------------------------------------------------------------------------------
+ pandterm = function(message) {
+ stop(message, call. = FALSE) }
+# ----------------------------------------------------------------------------------------------
+ rmultinomF=
+ function(p) {
+ return(sum(runif(1) > cumsum(p))+1)
+ }
+# -----------------------------------------------------------------------------------------------
+ alphaD=function(Prioralpha,Istar,gridsize){
+#
+# function to draw alpha using prior, p(alpha)= (1-(alpha-alphamin)/(alphamax-alphamin))**power
+#
+ power=Prioralpha$power
+ alphamin=Prioralpha$alphamin
+ alphamax=Prioralpha$alphamax
+ n=Prioralpha$n
+ alpha=seq(from=alphamin,to=(alphamax-0.000001),len=gridsize)
+ lnprob=Istar*log(alpha) + lgamma(alpha) - lgamma(n+alpha) +
+ power*log(1-(alpha-alphamin)/(alphamax-alphamin))
+ lnprob=lnprob-median(lnprob)
+ probs=exp(lnprob)
+ probs=probs/sum(probs)
+ return(alpha[rmultinomF(probs)])
+}
+# -----------------------------------------------------------------------------------------------
+#
+yden=function(thetaStar,y,eta){
+#
+# function to compute f(y | theta)
+# computes f for all values of theta in theta list of lists
+#
+# arguments:
+# thetaStar is a list of lists. thetaStar[[i]] is a list with components, mu, rooti
+# y |theta[[i]] ~ N(mu,(rooti %*% t(rooti))^-1) rooti is inverse of Chol root of Sigma
+# eta is not used
+#
+# output:
+# length(thetaStar) x n array of values of f(y[j,]|thetaStar[[i]]
+#
+
+nunique=length(thetaStar)
+n=nrow(y)
+ydenmat=matrix(double(n*nunique),ncol=n)
+k=ncol(y)
+for(i in 1:nunique){
+
+ # now compute vectorized version of lndMvn
+ # compute y_i'RIRI'y_i for all i
+ #
+ mu=thetaStar[[i]]$mu; rooti=thetaStar[[i]]$rooti
+ quads=colSums((crossprod(rooti,(t(y)-mu)))^2)
+ ydenmat[i,]=exp(-(k/2)*log(2*pi) + sum(log(diag(rooti))) - .5*quads)
+
+}
+return(ydenmat)
+}
+#
+#
+# -----------------------------------------------------------------------------------------
+#
+#
+GD=function(lambda){
+#
+# function to draw from prior for Multivariate Normal Model
+#
+# mu|Sigma ~ N(mubar,Sigma x Amu^-1)
+# Sigma ~ IW(nu,V)
+#
+#
+nu=lambda$nu
+V=lambda$V
+mubar=lambda$mubar
+Amu=lambda$Amu
+k=length(mubar)
+Sigma=rwishart(nu,chol2inv(chol(lambda$V)))$IW
+root=chol(Sigma)
+mu=mubar+(1/sqrt(Amu))*t(root)%*%matrix(rnorm(k),ncol=1)
+return(list(mu=as.vector(mu),rooti=backsolve(root,diag(k))))
+}
+#
+#
+# -------------------------------------------------------------------------------------------
+#
+#
+thetaD=function(y,lambda,eta){
+#
+# function to draw from posterior of theta given data y and base prior G0(lambda)
+#
+# here y ~ N(mu,Sigma)
+# theta = list(mu=mu,rooti=chol(Sigma)^-1)
+# mu|Sigma ~ N(mubar,Sigma (x) Amu-1)
+# Sigma ~ IW(nu,V)
+#
+# arguments:
+# y is n x k matrix of obs
+# lambda is list(mubar,Amu,nu,V)
+# eta is not used
+# output:
+# one draw of theta, list(mu,rooti)
+# Sigma=inv(rooti)%*%t(inv(rooti))
+#
+# note: we assume that y is a matrix. if there is only one obs, y is a 1 x k matrix
+#
+rout=rmultireg(y,matrix(c(rep(1,nrow(y))),ncol=1),matrix(lambda$mubar,nrow=1),matrix(lambda$Amu,ncol=1),
+ lambda$nu,lambda$V)
+return(list(mu=as.vector(rout$B),rooti=backsolve(chol(rout$Sigma),diag(ncol(y)))))
+}
+#
+# END OF REQUIRED FUNCTIONS AREA
+# --------------------------------------------------------------------------------------------
+#
+
+ n = length(theta)
+
+ eta=NULL # note eta is not used
+ thetaNp1=NULL
+
+ p=c(rep(1/(alpha+(n-1)),n-1),alpha/(alpha+(n-1)))
+
+ nunique=length(thetaStar)
+
+ if(nunique > maxuniq ) { pandterm("maximum number of unique thetas exceeded")}
+ ydenmat=matrix(double(maxuniq*n),ncol=n)
+ ydenmat[1:nunique,]=yden(thetaStar,y,eta)
+ # ydenmat is a length(thetaStar) x n array of density values given f(y[j,] | thetaStar[[i]]
+ # note: due to remix step (below) we must recompute ydenmat each time!
+
+ # use .Call to draw theta list
+ out= .Call("thetadraw",y,ydenmat,indic,q0v,p,theta,lambda,eta=eta,
+ thetaD=thetaD,yden=yden,maxuniq,nunique,new.env())
+
+ # theta has been modified by thetadraw so we need to recreate thetaStar
+ thetaStar=unique(theta)
+ nunique=length(thetaStar)
+
+ #thetaNp1 and remix
+ probs=double(nunique+1)
+ for(j in 1:nunique) {
+ ind = which(sapply(theta,identical,thetaStar[[j]]))
+ probs[j]=length(ind)/(alpha+n)
+ new_utheta=thetaD(y[ind,,drop=FALSE],lambda,eta)
+ for(i in seq(along=ind)) {theta[[ind[i]]]=new_utheta}
+ indic[ind]=j
+ thetaStar[[j]]=new_utheta
+ }
+ probs[nunique+1]=alpha/(alpha+n)
+ ind=rmultinomF(probs)
+ if(ind==length(probs)) {
+ thetaNp1=GD(lambda)
+ } else {
+ thetaNp1=thetaStar[[ind]]
+ }
+
+ #alpha
+ alpha=alphaD(Prioralpha,nunique,gridsize)
+
+ return(list(theta=theta,indic=indic,thetaStar=thetaStar,
+ thetaNp1=thetaNp1,alpha=alpha,Istar=nunique))
+}
+#
+#
+# -----------------------------------------------------------------------------------------
+#
+#
+q0=function(y,lambda,eta){
+#
+# function to compute a vector of int f(y[i]|theta) p(theta|lambda)dlambda
+# here p(theta|lambda) is G0 the base prior
+#
+# implemented for a multivariate normal data density and standard conjugate
+# prior:
+# theta=list(mu,Sigma)
+# f(y|theta) is N(mu,Sigma)
+# lambda=list(mubar,Amu,nu,V)
+# mu|Sigma ~ N(mubar,Sigma (x) Amu^-1)
+# Sigma ~ IW(nu,V)
+#
+# arguments:
+# Y is n x k matrix of observations
+# eta is not used
+# lambda=list(mubar,Amu,nu,V)
+#
+# output:
+# vector of q0 values for each obs (row of Y)
+#
+# p. rossi 12/05
+#
+# here y is matrix of observations (each row is an obs)
+
+mubar=lambda$mubar; nu=lambda$nu ; Amu=lambda$Amu; V=lambda$V
+k=ncol(y)
+R=chol(V)
+logdetR=sum(log(diag(R)))
+if (k > 1)
+ {lnk1k2=(k/2)*log(2)+log((nu-k)/2)+lgamma((nu-k)/2)-lgamma(nu/2)+sum(log(nu/2-(1:(k-1))/2))}
+else
+ {lnk1k2=(k/2)*log(2)+log((nu-k)/2)+lgamma((nu-k)/2)-lgamma(nu/2)}
+constant=-(k/2)*log(2*pi)+(k/2)*log(Amu/(1+Amu)) + lnk1k2 + nu*logdetR
+#
+# note: here we are using the fact that |V + S_i | = |R|^2 (1 + v_i'v_i)
+# where v_i = sqrt(Amu/(1+Amu))*t(R^-1)*(y_i-mubar), R is chol(V)
+#
+# and S_i = Amu/(1+Amu) * (y_i-mubar)(y_i-mubar)'
+#
+mat=sqrt(Amu/(1+Amu))*t(backsolve(R,diag(ncol(y))))%*%(t(y)-mubar)
+vivi=colSums(mat^2)
+
+lnq0v=constant-((nu+1)/2)*(2*logdetR+log(1+vivi))
+
+return(exp(lnq0v))
+}
+#
+#
+# --------------------------------------------------------------------------------------------
+#
+#
+# END OF REQUIRED FUNCTIONS AREA
+#
+#
+#initialize comps,indic,ncomp
+comps=unique(theta)
+ncomp=length(comps)
+indic=double(n)
+for(j in 1:ncomp){
+ indic[which(sapply(theta,identical,comps[[j]]))]=j
+ }
+# initialize eta
+eta=NULL
+#
+# initialize alpha
+alpha=1
+
+# reserve space for draws
+#
+deltadraw = matrix(double(floor(R/keep)*dimd),ncol=dimd)
+betadraw = rep(0.0,floor(R/keep))
+alphadraw=double(floor(R/keep))
+Istardraw=double(floor(R/keep))
+if(isgamma) gammadraw = matrix(double(floor(R/keep)*dimg),ncol=dimg)
+thetaNp1draw=vector("list",R)
+
+#
+# start main iteration loop
+#
+itime=proc.time()[3]
+cat("MCMC Iteration (est time to end -min) ",fill=TRUE)
+fsh()
+
+for(rep in 1:R) {
+
+ # draw beta and gamma
+ if(isgamma)
+ {out=get_ytxt(y=y,z=z,delta=delta,x=x,w=w,
+ ncomp=ncomp,indic=indic,comps=comps)}
+ else
+ {out=get_ytxt(y=y,z=z,delta=delta,x=x,
+ ncomp=ncomp,indic=indic,comps=comps)}
+
+ bg = breg(out$yt,out$xt,mbg,Abg)
+ beta = bg[1]
+ if(isgamma) gamma = bg[2:length(bg)]
+
+ # draw delta
+ if(isgamma)
+ {out=get_ytxtd(y=y,z=z,beta=beta,gamma=gamma,
+ x=x,w=w,ncomp=ncomp,indic=indic,comps=comps,dimd=dimd)}
+ else
+ {out=get_ytxtd(y=y,z=z,beta=beta,
+ x=x,ncomp=ncomp,indic=indic,comps=comps,dimd=dimd)}
+
+ delta = breg(out$yt,out$xtd,md,Ad)
+
+ # DP process stuff- theta | lambda
+ if(isgamma) {Err = cbind(x-z%*%delta,y-beta*x-w%*%gamma)}
+ else {Err = cbind(x-z%*%delta,y-beta*x)}
+ q0v = q0(Err,lambda,eta)
+ DPout=rthetaDP(maxuniq=maxuniq,alpha=alpha,lambda=lambda,Prioralpha=Prioralpha,theta=theta,
+ thetaStar=comps,indic=indic,q0v=q0v,y=Err,gridsize=gridsize)
+ indic=DPout$indic
+ theta=DPout$theta
+ comps=DPout$thetaStar
+ alpha=DPout$alpha
+ Istar=DPout$Istar
+ ncomp=length(comps)
+
+ 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
+ deltadraw[mkeep,]=delta
+ betadraw[mkeep]=beta
+ alphadraw[mkeep]=alpha
+ Istardraw[mkeep]=Istar
+ if(isgamma) gammadraw[mkeep,]=gamma
+ thetaNp1draw[[mkeep]]=list(DPout$thetaNp1)
+ }
+}
+#
+# rescale
+#
+if(SCALE){
+ deltadraw=deltadraw*scalex
+ betadraw=betadraw*scaley/scalex
+ if(isgamma) {gammadraw=gammadraw*scaley}
+}
+ctime = proc.time()[3]
+cat(' Total Time Elapsed: ',round((ctime-itime)/60,2),'\n')
+
+densitymix=list(matrix(c(rep(1,length(thetaNp1draw))),ncol=1),list("gunk"),thetaNp1draw)
+#
+# densitymix is in the format to be used with the generic mixture of normals plotting
+# methods (plot.bayesm.nmix)
+#
+attributes(densitymix)$class=c("bayesm.nmix")
+
+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(alphadraw)$class=c("bayesm.mat","mcmc")
+attributes(alphadraw)$mcpar=c(1,R,keep)
+attributes(Istardraw)$class=c("bayesm.mat","mcmc")
+attributes(Istardraw)$mcpar=c(1,R,keep)
+if(isgamma){
+ attributes(gammadraw)$class=c("bayesm.mat","mcmc")
+ attributes(gammadraw)$mcpar=c(1,R,keep)}
+
+if(isgamma)
+ { return(list(deltadraw=deltadraw,betadraw=betadraw,alphadraw=alphadraw,Istardraw=Istardraw,
+ gammadraw=gammadraw,densitymix=densitymix))}
+ else
+ { return(list(deltadraw=deltadraw,betadraw=betadraw,alphadraw=alphadraw,Istardraw=Istardraw,
+ densitymix=densitymix))}
+}
+
+
diff --git a/R/rnegbinRw.R b/R/rnegbinRw.R
index 79cf0a1..1ebe6c4 100755
--- a/R/rnegbinRw.R
+++ b/R/rnegbinRw.R
@@ -75,6 +75,7 @@ if(is.null(Data$y)) {pandterm("Requires Data element y")} else {y=Data$y}
nvar = ncol(X)
if (length(y) != nrow(X)) {pandterm("Mismatch in the number of observations in X and y")}
+nobs=length(y)
#
# check for prior elements
diff --git a/R/summary.bayesm.mat.R b/R/summary.bayesm.mat.R
index c09b1a4..3bef1b6 100755
--- a/R/summary.bayesm.mat.R
+++ b/R/summary.bayesm.mat.R
@@ -27,7 +27,7 @@ summary.bayesm.mat=function(object,names,burnin=trunc(.1*nrow(X)),tvalues,QUANTI
rownames(mat)[2]="std dev"
rownames(mat)[3]="num se"
rownames(mat)[4]="rel eff"
- rownames(mat)[5]="s size"
+ rownames(mat)[5]="sam 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")
diff --git a/R/summary.bayesm.nmix.R b/R/summary.bayesm.nmix.R
index bbb2a6d..569ed9f 100755
--- a/R/summary.bayesm.nmix.R
+++ b/R/summary.bayesm.nmix.R
@@ -25,9 +25,9 @@ for(i in (burnin+1):R){
cat("\nNormal Mixture Moments\n Mean\n")
attributes(mumat)$class="bayesm.mat"
attributes(sigmat)$class="bayesm.var"
-summary(mumat,names,QUANTILES=FALSE,TRAILER=FALSE)
+summary(mumat,names,burnin=burnin,QUANTILES=FALSE,TRAILER=FALSE)
cat(" \n")
-summary(sigmat)
+summary(sigmat,burnin=burnin)
cat("note: 1st and 2nd Moments for a Normal Mixture \n")
cat(" may not be interpretable, consider plots\n")
invisible()
diff --git a/inst/doc/bayesm-manual.pdf b/inst/doc/bayesm-manual.pdf
index f19cf6c..502dbc0 100755
Binary files a/inst/doc/bayesm-manual.pdf and b/inst/doc/bayesm-manual.pdf differ
diff --git a/man/plot.bayesm.mat.Rd b/man/plot.bayesm.mat.Rd
index 259f490..b9181ee 100755
--- a/man/plot.bayesm.mat.Rd
+++ b/man/plot.bayesm.mat.Rd
@@ -9,7 +9,7 @@
array correspond to parameters and the rows to MCMC draws.
}
\usage{
-\method{plot}{bayesm.mat}(x,names,burnin,tvalues,TRACEPLOT,DEN,INT, ...)
+\method{plot}{bayesm.mat}(x,names,burnin,tvalues,TRACEPLOT,DEN,INT,CHECK_NDRAWS, ...)
}
\arguments{
\item{x}{ An object of either S3 class, bayesm.mat, or S3 class, mcmc }
@@ -19,6 +19,7 @@
\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{CHECK_NDRAWS}{ logical, TRUE check that there are at least 100 draws, def: TRUE }
\item{...}{ standard graphics parameters }
}
\details{
diff --git a/man/plot.bayesm.nmix.Rd b/man/plot.bayesm.nmix.Rd
index 489377d..757a7f7 100755
--- a/man/plot.bayesm.nmix.Rd
+++ b/man/plot.bayesm.nmix.Rd
@@ -10,15 +10,20 @@
are produced.
}
\usage{
-\method{plot}{bayesm.nmix}(x,names,burnin,Grid,bi.sel,nstd, ...)
+\method{plot}{bayesm.nmix}(x,names,burnin,Grid,bi.sel,nstd,marg,Data,ngrid,ndraw, ...)
}
\arguments{
- \item{x}{ An object of either S3 class, bayesm.nmix, or S3 class, mcmc }
+ \item{x}{ An object of S3 class bayesm.nmix }
\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{Grid}{matrix of grid points for densities, def: mean +/- nstd std deviations (if Data no supplied),
+ range of Data if supplied)}
\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{marg}{logical, if TRUE display marginals, def: TRUE}
+ \item{Data}{matrix of data points, used to paint histograms on marginals and for grid }
+ \item{ngrid}{number of grid points for density estimates, def:50}
+ \item{ndraw}{number of draws to average Mcmc estimates over, def:200}
\item{...}{ standard graphics parameters }
}
\details{
@@ -33,13 +38,14 @@
\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}}}
+\seealso{ \code{\link{rnmixGibbs}}, \code{\link{rhierMnlRwMixture}}, \code{\link{rhierLinearMixture}},
+ \code{\link{rDPGibbs}}}
\examples{
##
## not run
-# out=rnmix(Data,Prior,Mcmc)
+# out=rnmixGibbs(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
+# # plot bivariate distributions for dimension 1,2; 3,4; and 1,3
#
}
diff --git a/man/rDPGibbs.Rd b/man/rDPGibbs.Rd
new file mode 100755
index 0000000..0b4e33c
--- /dev/null
+++ b/man/rDPGibbs.Rd
@@ -0,0 +1,212 @@
+\name{rDPGibbs}
+\alias{rDPGibbs}
+\concept{bayes}
+\concept{MCMC}
+\concept{normal mixtures}
+\concept{Dirichlet Process}
+\concept{Gibbs Sampling}
+\title{ Density Estimation with Dirichlet Process Prior and Normal Base }
+\description{
+ \code{rDPGibbs} implements a Gibbs Sampler to draw from the posterior for a normal mixture problem
+ with a Dirichlet Process prior. The base distribution is a multivariate normal distribution. A natural
+ conjugate base prior is used with priors on the hyper parameters of this distribution. One interpretation
+ of this model is as a normal mixture with a random number of components.
+}
+
+\usage{
+rDPGibbs(Prior, Data, Mcmc)
+}
+
+\arguments{
+ \item{Prior}{ list(Prioralpha,lambda\_hyper) }
+ \item{Data}{ list(y) }
+ \item{Mcmc}{ list(R,keep,maxuniq,SCALE,gridsize) }
+}
+
+\details{
+
+Model: \cr
+ \eqn{y_i} \eqn{\sim}{~} \eqn{N(mu_i,Sigma_i)}. \cr
+
+Priors:\cr
+ \eqn{theta_i=(mu_i,Sigma_i)} \eqn{\sim}{~} \eqn{DP(G_0(lambda),alpha)}\cr
+ \eqn{G_0(lambda):}\cr
+ \eqn{mu_i given Sigma_i} \eqn{\sim}{~} \eqn{N(0,Sigma_i (x) a^{-1})}\cr
+ \eqn{Sigma_i} \eqn{\sim}{~} \eqn{IW(nu,nu*v*I)}
+
+ \eqn{lambda(a,nu,v):}\cr
+ \eqn{a} \eqn{sim}{~} \eqn{uniform on grid[alim[1],alimb[2]]}\cr
+ \eqn{nu} \eqn{sim}{~} \eqn{uniform on grid[dim(data)-1 + exp(nulim[1]),dim(data)-1 +exp(nulim[2])]}\cr
+ \eqn{v} \eqn{sim}{~} \eqn{uniform on grid[vlim[1],vlim[2]]}
+
+ \eqn{alpha} \eqn{sim}{~} \eqn{(1-(alpha-alphamin)/(alphamax-alphamin))^power} \cr
+ alpha= alphamin then expected number of compoents = Istarmin \cr
+ alpha= alphamax then expected number of compoents = Istarmax \cr
+
+list arguments
+
+Data:\cr
+ \itemize{
+ \item{y}{N x k matrix of observations on k dimensional data}
+ }
+
+Prioralpha:\cr
+ \itemize{
+ \item{Istarmin}{expected number of components at lower bound of support of alpha}
+ \item{Istarmax}{expected number of components at upper bound of support of alpha}
+ \item{power}{power parameter for alpha prior}
+ }
+
+lambda\_hyper:\cr
+ \itemize{
+ \item{alim}{defines support of a distribution,def:c(.01,10) }
+ \item{nulim}{defines support of nu distribution, def:c(.01,3)}
+ \item{vlim}{defines support of v distribution, def:c(.1,4)}
+ }
+Mcmc:\cr
+ \itemize{
+ \item{R}{number of mcmc draws}
+ \item{keep}{thinning parm, keep every keepth draw}
+ \item{maxuniq}{storage constraint on the number of unique components}
+ \item{SCALE}{should data be scaled by mean,std deviation before posterior draws, def: TRUE}
+ \item{gridsize}{number of discrete points for hyperparameter priors,def: 20}
+ }
+
+output:\cr
+
+the basic output are draws from the predictive distribution of the data in the object, \code{nmix}.
+The average of these draws is the Bayesian analogue of a density estimate.
+
+nmix:\cr
+ \itemize{
+ \item{probdraw}{R/keep x 1 matrix of 1s}
+ \item{zdraw}{R/keep x N matrix of draws of indicators of which component each obs is assigned to}
+ \item{compdraw}{R/keep list of draws of normals}
+ }
+ Output of the components is in the form of a list of lists. \cr
+ compdraw[[i]] is ith draw -- list of lists. \cr
+ compdraw[[i]][[1]] is list of parms for normal component. \cr
+ compdraw[[i]][1]][[1]] is the mean vector.
+ \eqn{Sigma} = t(R)\%*\%R, \eqn{R^{-1}} = compdraw[[i]][[1]][[2]].
+}
+
+
+\note{
+ we parameterize the prior on \eqn{Sigma_i} such that \eqn{mode(Sigma)= nu/(nu+2) vI}.
+ The support of nu enforces valid IW density; \eqn{nulim[1] > 0}
+
+ We use the structure for \code{nmix} that is compatible with the \code{bayesm} routines for finite mixtures of normals.
+ This allows us to use the same summary and plotting methods.
+
+ The default choices of alim,nulim, and vlim determine the location and approximate size of candidate
+ "atoms" or possible normal components. The defaults are sensible given that we scale the data. Without
+ scaling, you want to insure that alim is set for a wide enough range of values (remember a is a precision
+ parameter) and the v is big enough to propose Sigma matrices wide enough to cover the data range.
+
+ A careful analyst should look at the posterior distribution of a, nu, v to make sure that the support is
+ set correctly in alim, nulim, vlim. In other words, if we see the posterior bunched up at one end of these
+ support ranges, we should widen the range and rerun.
+
+ If you want to force the procedure to use many small atoms, then set nulim to consider only large values and
+ set vlim to consider only small scaling constants. Set Istarmax to a large number. This will create a very
+ "lumpy" density estimate somewhat like the classical Kernel density estimates. Of course, this is not advised
+ if you have a prior belief that densities are relatively smooth.
+
+}
+
+
+\value{
+ \item{nmix}{a list containing: probdraw,zdraw,compdraw}
+ \item{alphadraw}{vector of draws of DP process tightness parameter}
+ \item{nudraw}{vector of draws of base prior hyperparameter}
+ \item{adraw}{vector of draws of base prior hyperparameter}
+ \item{vdraw}{vector of draws of base prior hyperparameter}
+}
+
+\author{ Peter Rossi, Graduate School of Business, University of Chicago,
+ \email{Peter.Rossi at ChicagoGsb.edu}.
+}
+
+\seealso{ \code{\link{rnmixGibbs}},\code{\link{rmixture}}, \code{\link{rmixGibbs}} ,
+ \code{\link{eMixMargDen}}, \code{\link{momMix}}, \code{\link{mixDen}}, \code{\link{mixDenBi}}}
+
+\examples{
+## simulate univariate data from Chi-Sq
+
+if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=2000} else {R=10}
+
+set.seed(66)
+N=200
+chisqdf=8; y1=as.matrix(rchisq(N,df=chisqdf))
+
+## set arguments for rDPGibbs
+
+Data1=list(y=y1)
+Prioralpha=list(Istarmin=1,Istarmax=10,power=.8)
+Prior1=list(alpha=1,Prioralpha=Prioralpha)
+
+Mcmc=list(R=R,keep=1,maxuniq=200)
+
+out1=rDPGibbs(Prior=Prior1,Data=Data1,Mcmc)
+
+if(0){
+## plotting examples
+rgi=c(0,20); grid=matrix(seq(from=rgi[1],to=rgi[2],length.out=50),ncol=1)
+deltax=(rgi[2]-rgi[1])/nrow(grid)
+plot(out1$nmix,Grid=grid,Data=y1)
+## plot true density with historgram
+plot(range(grid[,1]),1.5*range(dchisq(grid[,1],df=chisqdf)),type="n",xlab=paste("Chisq ; ",N," obs",sep=""), ylab="")
+hist(y1,xlim=rgi,freq=FALSE,col="yellow",breaks=20,add=TRUE)
+lines(grid[,1],dchisq(grid[,1],df=chisqdf)/(sum(dchisq(grid[,1],df=chisqdf))*deltax),col="blue",lwd=2)
+}
+
+
+## simulate bivariate data from the "Banana" distribution (Meng and Barnard)
+banana=function(A,B,C1,C2,N,keep=10,init=10)
+{ R=init*keep+N*keep
+ x1=x2=0
+ bimat=matrix(double(2*N),ncol=2)
+ for (r in 1:R)
+ { x1=rnorm(1,mean=(B*x2+C1)/(A*(x2^2)+1),sd=sqrt(1/(A*(x2^2)+1)))
+ x2=rnorm(1,mean=(B*x2+C2)/(A*(x1^2)+1),sd=sqrt(1/(A*(x1^2)+1)))
+ if (r>init*keep && r\%\%keep==0) {mkeep=r/keep; bimat[mkeep-init,]=c(x1,x2)} }
+return(bimat)
+}
+
+
+set.seed(66)
+nvar2=2
+A=0.5; B=0; C1=C2=3
+y2=banana(A=A,B=B,C1=C1,C2=C2,1000)
+
+Data2=list(y=y2)
+Prioralpha=list(Istarmin=1,Istarmax=10,power=.8)
+Prior2=list(alpha=1,Prioralpha=Prioralpha)
+R=2000
+Mcmc=list(R=R,keep=1,maxuniq=200)
+
+out2=rDPGibbs(Prior=Prior2,Data=Data2,Mcmc)
+
+
+if(0){
+## plotting examples
+
+rx1=range(y2[,1]); rx2=range(y2[,2])
+x1=seq(from=rx1[1],to=rx1[2],length.out=50)
+x2=seq(from=rx2[1],to=rx2[2],length.out=50)
+grid=cbind(x1,x2)
+
+plot(out2$nmix,Grid=grid,Data=y2)
+
+## plot true bivariate density
+tden=matrix(double(50*50),ncol=50)
+for (i in 1:50){ for (j in 1:50)
+ {tden[i,j]=exp(-0.5*(A*(x1[i]^2)*(x2[j]^2)+(x1[i]^2)+(x2[j]^2)-2*B*x1[i]*x2[j]-2*C1*x1[i]-2*C2*x2[j]))}
+}
+tden=tden/sum(tden)
+image(x1,x2,tden,col=terrain.colors(100),xlab="",ylab="")
+contour(x1,x2,tden,add=TRUE,drawlabels=FALSE)
+title("True Density")
+}
+}
+\keyword{ multivariate }
diff --git a/man/rhierNegbinRw.Rd b/man/rhierNegbinRw.Rd
index 3f09a34..8cc5f04 100755
--- a/man/rhierNegbinRw.Rd
+++ b/man/rhierNegbinRw.Rd
@@ -103,9 +103,9 @@ nz=2 # Number of Z variables
# Construct the Z matrix
Z = cbind(rep(1,nreg),rnorm(nreg,mean=1,sd=0.125))
-Delta = cbind(c(0.4,0.2), c(0.1,0.05))
+Delta = cbind(c(4,2), c(0.1,-1))
alpha = 5
-Vbeta = rbind(c(0.1,0),c(0,0.1))
+Vbeta = rbind(c(2,1),c(1,2))
# Construct the regdata (containing X)
simnegbindata = NULL
diff --git a/man/rivDP.Rd b/man/rivDP.Rd
new file mode 100755
index 0000000..775323d
--- /dev/null
+++ b/man/rivDP.Rd
@@ -0,0 +1,148 @@
+\name{rivDP}
+\alias{rivDP}
+\concept{Instrumental Variables}
+\concept{Gibbs Sampler}
+\concept{Dirichlet Process}
+\concept{bayes}
+\concept{endogeneity}
+\concept{simultaneity}
+\concept{MCMC}
+
+\title{ Linear "IV" Model with DP Process Prior for Errors}
+\description{
+ \code{rivDP} is a Gibbs Sampler for a linear structural equation with an arbitrary number of instruments.
+ \code{rivDP} uses a mixture of normals for the structural and reduced form equation implemented with a
+ Dirichlet Process Prior.
+}
+\usage{
+rivDP(Data, Prior, Mcmc)
+}
+\arguments{
+ \item{Data}{ list(z,w,x,y) }
+ \item{Prior}{ list(md,Ad,mbg,Abg,lambda,Prioralpha) (optional) }
+ \item{Mcmc}{ list(R,keep,SCALE) (R required) }
+}
+\details{
+ Model:\cr
+ \eqn{x=z'delta + e1}. \cr
+ \eqn{y=beta*x + w'gamma + e2}. \cr
+ \eqn{e1,e2} \eqn{\sim}{~} \eqn{N(theta_{i})}. \eqn{theta_{i}} represents \eqn{mu_{i},Sigma_{i}}
+
+ Note: Error terms have non-zero means. DO NOT include intercepts in the z or w matrices. This is different
+ from \code{rivGibbs} which requires intercepts to be included explicitly.
+
+ Priors:\cr
+ \eqn{delta} \eqn{\sim}{~} \eqn{N(md,Ad^{-1})}. \eqn{vec(beta,gamma)} \eqn{\sim}{~} \eqn{N(mbg,Abg^{-1})} \cr
+
+ \eqn{theta_{i}\sim{~}G} \cr
+
+ \eqn{G} \eqn{\sim}{~} \eqn{DP(alpha,G_{0})} \cr
+
+ \eqn{G_{0}} is the natural conjugate prior for \eqn{(mu,Sigma)}: \cr
+ \eqn{Sigma} \eqn{\sim}{~} \eqn{IW(nu,vI)} and \eqn{mu | Sigma} \eqn{\sim}{~} \eqn{N(0,1/amu Sigma)} \cr
+ These parameters are collected together in the list \code{lambda}. It is highly
+ recommended that you use the default settings for these hyper-parameters.\cr
+
+ \eqn{alpha} \eqn{\sim}{~} \eqn{(1-(alpha-alpha_{min})/(alpha_{max}-alpha{min}))^{omega}} \cr
+ where \eqn{alpha_{min}} and \eqn{alpha_{max}} are set using the arguments in the reference
+ below. It is highly recommended that you use the default values for the hyperparameters
+ of the prior on alpha
+
+ List arguments contain:
+ \itemize{
+ \item{\code{z}}{ matrix of obs on instruments}
+ \item{\code{y}}{ vector of obs on lhs var in structural equation}
+ \item{\code{x}}{ "endogenous" var in structural eqn}
+ \item{\code{w}}{ matrix of obs on "exogenous" vars in the structural eqn}
+ \item{\code{md}}{ prior mean of delta (def: 0)}
+ \item{\code{Ad}}{ pds prior prec for prior on delta (def: .01I)}
+ \item{\code{mbg}}{ prior mean vector for prior on beta,gamma (def: 0)}
+ \item{\code{Abg}}{ pds prior prec for prior on beta,gamma (def: .01I)}
+ \item{\code{lambda}}{ list of hyperparameters for theta prior- use default settings }
+ \item{\code{Prioralpha}}{ list of hyperparameters for theta prior- use default settings }
+ \item{\code{R}}{ number of MCMC draws}
+ \item{\code{keep}}{ MCMC thinning parm: keep every keepth draw (def: 1)}
+ \item{\code{SCALE}}{ scale data, def: TRUE}
+ \item{\code{gridsize}}{ gridsize parm for alpha draws (def: 20)}
+ }
+}
+\value{
+ a list containing:
+ \item{deltadraw}{R/keep x dim(delta) array of delta draws}
+ \item{betadraw}{R/keep x 1 vector of beta draws}
+ \item{gammadraw}{R/keep x dim(gamma) array of gamma draws }
+ \item{Istardraw}{R/keep x 1 array of drawsi of the number of unique normal components}
+ \item{alphadraw}{R/keep x 1 array of draws of Dirichlet Process tightness parameter}
+ \item{densitymix}{R/keep x list of draws for predictive distribution of errors}
+}
+\references{ For further discussion, see "A Semi-Parametric Bayesian Approach to the Instrumental
+ Variable Problem," by Conley, Hansen, McCulloch and Rossi, Journal of Econometrics (2008).\cr
+}
+\seealso{\code{rivGibbs}}
+
+\author{ Peter Rossi, Graduate School of Business, University of Chicago,
+ \email{Peter.Rossi at ChicagoGsb.edu}.
+}
+\examples{
+##
+if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=2000} else {R=10}
+
+##
+## simulate scaled log-normal errors and run
+##
+set.seed(66)
+k=10
+delta=1.5
+Sigma=matrix(c(1,.6,.6,1),ncol=2)
+N=1000
+tbeta=4
+set.seed(66)
+scalefactor=.6
+root=chol(scalefactor*Sigma)
+mu=c(1,1)
+##
+## compute interquartile ranges
+##
+ninterq=qnorm(.75)-qnorm(.25)
+error=matrix(rnorm(100000*2),ncol=2)%*%root
+error=t(t(error)+mu)
+Err=t(t(exp(error))-exp(mu+.5*scalefactor*diag(Sigma)))
+lnNinterq=quantile(Err[,1],prob=.75)-quantile(Err[,1],prob=.25)
+##
+## simulate data
+##
+error=matrix(rnorm(N*2),ncol=2)\%*\%root
+error=t(t(error)+mu)
+Err=t(t(exp(error))-exp(mu+.5*scalefactor*diag(Sigma)))
+#
+# scale appropriately
+Err[,1]=Err[,1]*ninterq/lnNinterq
+Err[,2]=Err[,2]*ninterq/lnNinterq
+z=matrix(runif(k*N),ncol=k)
+x=z\%*\%(delta*c(rep(1,k)))+Err[,1]
+y=x*tbeta+Err[,2]
+
+# set intial values for MCMC
+Data = list(); Mcmc=list()
+Data$z = z; Data$x=x; Data$y=y
+
+# start MCMC and keep results
+Mcmc$maxuniq=100
+Mcmc$R=R
+end=Mcmc$R
+begin=100
+
+out=rivDP(Data=Data,Mcmc=Mcmc)
+
+cat("Summary of Beta draws",fill=TRUE)
+summary(out$betadraw,tvalues=tbeta)
+
+if(0){
+## plotting examples
+plot(out$betadraw,tvalues=tbeta)
+plot(out$densitymix) ## plot "fitted" density of the errors
+##
+
+}
+}
+\keyword{ models }
diff --git a/man/rnmixGibbs.Rd b/man/rnmixGibbs.Rd
index 5f906da..ad9776b 100755
--- a/man/rnmixGibbs.Rd
+++ b/man/rnmixGibbs.Rd
@@ -109,7 +109,7 @@ if(0){
##
## plotting examples
##
-plot(out)
+plot(out,Data=dm$x)
}
}
diff --git a/man/rscaleUsage.Rd b/man/rscaleUsage.Rd
index 2d7a030..2a3930f 100755
--- a/man/rscaleUsage.Rd
+++ b/man/rscaleUsage.Rd
@@ -67,7 +67,7 @@ rscaleUsage(Data,Prior, Mcmc)
\examples{
##
-if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=1000} else {R=5}
+if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=1000} else {R=1}
{
data(customerSat)
surveydat = list(k=10,x=as.matrix(customerSat))
diff --git a/man/tuna.Rd b/man/tuna.Rd
index fb8c604..71f98ba 100755
--- a/man/tuna.Rd
+++ b/man/tuna.Rd
@@ -1,10 +1,11 @@
\name{tuna}
\alias{tuna}
\docType{data}
-\title{Sales Data on Canned Tuna Sales}
+\title{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.
+ Weekly data aggregated to the chain level. This data is extracted from the Dominick's Finer Foods database
+ maintained by the University of Chicago \url{http://http://research.chicagogsb.edu/marketing/databases/dominicks/dataset.aspx}.
Brands are seven of the top 10 UPCs in the canned tuna product category.
}
\usage{data(tuna)}
diff --git a/src/thetadraw.c b/src/thetadraw.c
new file mode 100755
index 0000000..7a2b417
--- /dev/null
+++ b/src/thetadraw.c
@@ -0,0 +1,150 @@
+#include <R.h>
+#include <Rmath.h>
+#include <math.h>
+#include <Rinternals.h>
+#include <Rdefines.h>
+
+/* modified by rossi 7/06 to remove thetaStar argument and copy */
+/* modified by rossi 7/06 to remove theta copy and modify directly */
+
+/* function to make multinomial draws */
+
+int rmultin(double *probs, int nprob )
+{
+ double cumprob,rnd;
+ int i;
+ GetRNGstate();
+ rnd=unif_rand();
+ cumprob=0.0;
+ for (i=0; i < nprob; i++){
+ if(rnd > cumprob && rnd <= cumprob+probs[i])
+ { break; }
+ else
+ {cumprob=cumprob+probs[i];}
+ }
+ PutRNGstate();
+ return i+1 ;
+}
+
+
+/* gets a row of a matrix and returns this as a matrix by setting dim */
+
+SEXP getrow(SEXP mat, int row, int nrow, int ncol){
+ int i,ind;
+ SEXP ans, ndim;
+ PROTECT(ans=NEW_NUMERIC(ncol));
+ PROTECT(ndim=NEW_INTEGER(2));
+ for(i =0; i < ncol; i++){
+ ind=i*nrow+row;
+ NUMERIC_POINTER(ans)[i]=NUMERIC_POINTER(mat)[ind];
+ }
+ INTEGER_POINTER(ndim)[0]=1;
+ INTEGER_POINTER(ndim)[1]=ncol;
+ SET_DIM(ans,ndim);
+ UNPROTECT(2);
+return(ans);
+}
+
+
+/* theta draw routine to be used with .Call */
+
+SEXP thetadraw( SEXP y, SEXP ydenmatO, SEXP indicO, SEXP q0v, SEXP p,
+ SEXP theta, SEXP lambda, SEXP eta,
+ SEXP thetaD, SEXP yden,
+ SEXP maxuniqS,SEXP nuniqueS,
+ SEXP rho) {
+ int nunique,n,ncol,j,i,maxuniq,inc,index,ii,jj,ind ;
+ SEXP R_fc_thetaD, R_fc_yden, yrow, ydim, onetheta, lofone, newrow,
+ ydenmat, ydendim ;
+ double *probs;
+ int *indmi;
+ int *indic;
+ double sprob;
+
+ nunique=INTEGER_VALUE(nuniqueS);
+ n=length(theta);
+ maxuniq=INTEGER_VALUE(maxuniqS);
+
+ /* create new lists for use and output */
+ PROTECT(lofone=NEW_LIST(1));
+
+ /* create R function call object, lang4 creates a pairwise (linked) list with
+ 4 values -- function, first arg, sec arg, third arg. R_NilValue is a placeholder until
+ we associate first argument (which varies in our case) */
+ PROTECT(R_fc_thetaD=lang4(thetaD,R_NilValue,lambda,eta));
+ PROTECT(R_fc_yden=lang4(yden,R_NilValue,y,eta));
+
+ PROTECT(ydim=GET_DIM(y));
+ ncol=INTEGER_POINTER(ydim)[1];
+ PROTECT(yrow=NEW_NUMERIC(ncol));
+ PROTECT(newrow=NEW_NUMERIC(ncol));
+ PROTECT(ydenmat=NEW_NUMERIC(maxuniq*n));
+ PROTECT(ydendim=NEW_INTEGER(2));
+ INTEGER_POINTER(ydendim)[0]=maxuniq;
+ INTEGER_POINTER(ydendim)[1]=n;
+
+ /* copy iformation from R objects that will be modified
+ note that we must access elements in the lists (generic vectors) by using VECTOR_ELT
+ we can't use the pointer and deferencing directly like we can for numeric and integer
+ vectors */
+ for(j=0;j < maxuniq*n; j++){NUMERIC_POINTER(ydenmat)[j]=NUMERIC_POINTER(ydenmatO)[j];}
+ SET_DIM(ydenmat,ydendim);
+
+ /* allocate space for local vectors */
+ probs=(double *)R_alloc(n,sizeof(double));
+ indmi=(int *)R_alloc((n-1),sizeof(int));
+ indic=(int *)R_alloc(n,sizeof(int));
+
+ /* copy information from R object indicO to indic */
+ for(j=0;j < n; j++) {indic[j]=NUMERIC_POINTER(indicO)[j];}
+
+ /* start loop over observations */
+
+ for(i=0;i < n; i++){
+ probs[n-1]=NUMERIC_POINTER(q0v)[i]*NUMERIC_POINTER(p)[n-1];
+
+ /* make up indmi -- vector of length n-1 consisting of -i as in R notation --
+ 1, ...,i-1, ,i+1,...,n */
+ inc=0;
+ for(j=0;j < (n-1); j++){
+ if(j==i) {inc=inc+1;};
+ indmi[j]=inc;
+ inc=inc+1;
+ }
+ for(j=0;j < (n-1); j++){
+ ii=indic[indmi[j]]; jj=i; /* find element ydenmat(ii,jj+1) */
+ index=jj*maxuniq+(ii-1);
+ probs[j]=NUMERIC_POINTER(p)[j]*NUMERIC_POINTER(ydenmat)[index];
+ }
+ sprob=0.0;
+ for(j=0;j<n;j++){sprob=sprob+probs[j];}
+ for(j=0;j<n;j++){probs[j]=probs[j]/sprob;}
+ ind=rmultin(probs,n);
+
+ if(ind == n){
+ yrow=getrow(y,i,n,ncol);
+ SETCADR(R_fc_thetaD,yrow); /* set the second argument to yrow -- head of the tail */
+ onetheta=eval(R_fc_thetaD,rho);
+ SET_ELEMENT(theta,i,onetheta);
+ if((nunique) > (maxuniq-1)) {error("max number of unique thetas exceeded");}
+ /* check to make sure we don't exceed max number of unique theta */
+ SET_ELEMENT(lofone,0,onetheta);
+ SETCADR(R_fc_yden,lofone);
+ newrow=eval(R_fc_yden,rho);
+ for(j=0;j<n; j++)
+ { NUMERIC_POINTER(ydenmat)[j*maxuniq+nunique]=NUMERIC_POINTER(newrow)[j];}
+ indic[i]=nunique+1;
+ nunique=nunique+1;
+ }
+ else {
+ onetheta=VECTOR_ELT(theta,indmi[ind-1]);
+ SET_ELEMENT(theta,i,onetheta);
+ indic[i]=indic[indmi[ind-1]];
+ }
+ }
+
+ UNPROTECT(8);
+ return(nuniqueS); /* returns argument -- function now is called for its effect on theta */
+}
+
+
--
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