R : Copyright 2005, The R Foundation for Statistical Computing Version 2.1.1 (2005-06-20), ISBN 3-900051-07-0 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for a HTML browser interface to help. Type 'q()' to quit R. > ### *
> ### > attach(NULL, name = "CheckExEnv") > assign(".CheckExEnv", as.environment(2), pos = length(search())) # base > ## add some hooks to label plot pages for base and grid graphics > setHook("plot.new", ".newplot.hook") > setHook("persp", ".newplot.hook") > setHook("grid.newpage", ".gridplot.hook") > > assign("cleanEx", + function(env = .GlobalEnv) { + rm(list = ls(envir = env, all.names = TRUE), envir = env) + RNGkind("default", "default") + set.seed(1) + options(warn = 1) + delayedAssign("T", stop("T used instead of TRUE"), + assign.env = .CheckExEnv) + delayedAssign("F", stop("F used instead of FALSE"), + assign.env = .CheckExEnv) + sch <- search() + newitems <- sch[! sch %in% .oldSearch] + for(item in rev(newitems)) + eval(substitute(detach(item), list(item=item))) + missitems <- .oldSearch[! .oldSearch %in% sch] + if(length(missitems)) + warning("items ", paste(missitems, collapse=", "), + " have been removed from the search path") + }, + env = .CheckExEnv) > assign("..nameEx", "__{must remake R-ex/*.R}__", env = .CheckExEnv) # for now > assign("ptime", proc.time(), env = .CheckExEnv) > grDevices::postscript("bayesm-Examples.ps") > assign("par.postscript", graphics::par(no.readonly = TRUE), env = .CheckExEnv) > options(contrasts = c(unordered = "contr.treatment", ordered = "contr.poly")) > options(warn = 1) > library('bayesm') > > assign(".oldSearch", search(), env = .CheckExEnv) > assign(".oldNS", loadedNamespaces(), env = .CheckExEnv) > cleanEx(); ..nameEx <- "Scotch" > > ### * Scotch > > flush(stderr()); flush(stdout()) > > ### Name: Scotch > ### Title: Survey Data on Brands of Scotch Consumed > ### Aliases: Scotch > ### Keywords: datasets > > ### ** Examples > > data(Scotch) > cat(" Frequencies of Brands", fill=TRUE) Frequencies of Brands > mat=apply(as.matrix(Scotch),2,mean) > print(mat) Chivas.Regal Dewar.s.White.Label 0.36339044 0.23309288 Johnnie.Walker.Black.Label J...B 0.22633003 0.20649234 Johnnie.Walker.Red.Label Other.Brands 0.19116321 0.18665464 Glenlivet Cutty.Sark 0.15960325 0.15284040 Glenfiddich Pinch..Haig. 0.15058611 0.05275023 Clan.MacGregor Ballantine 0.04643823 0.04463481 Macallan Passport 0.04283138 0.03697024 Black...White Scoresby.Rare 0.03651939 0.03561767 Grants Ushers 0.03336339 0.03020739 White.Horse Knockando 0.02795311 0.02119026 the.Singleton 0.01397656 > ## > ## use Scotch data to run Multivariate Probit Model > ## > if(nchar(Sys.getenv("LONG_TEST")) != 0){ + ## + + y=as.matrix(Scotch) + p=ncol(y); n=nrow(y) + dimnames(y)=NULL + y=as.vector(t(y)) + y=as.integer(y) + I_p=diag(p) + X=rep(I_p,n) + X=matrix(X,nrow=p) + X=t(X) + + R=2000 + Data=list(p=p,X=X,y=y) + Mcmc=list(R=R) + set.seed(66) + 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) + rdraw=matrix(double((R)*p*p),ncol=p*p) + rdraw=t(apply(out$sigmadraw,1,nmat)) + 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)) + + } > > > > > cleanEx(); ..nameEx <- "bank" > > ### * bank > > flush(stderr()); flush(stdout()) > > ### Name: bank > ### Title: Bank Card Conjoint Data of Allenby and Ginter (1995) > ### Aliases: bank > ### Keywords: datasets > > ### ** Examples > > data(bank) > cat(" table of Binary Dep Var", fill=TRUE) table of Binary Dep Var > print(table(bank$choiceAtt[,2])) 0 1 8326 6473 > cat(" table of Attribute Variables",fill=TRUE) table of Attribute Variables > mat=apply(as.matrix(bank$choiceAtt[,3:16]),2,table) > print(mat) $Med_FInt -1 0 1 2381 10409 2009 $Low_FInt -1 0 1 1181 12185 1433 $Med_VInt 0 1 13252 1547 $Rewrd_2 -1 0 1 2129 10427 2243 $Rewrd_3 -1 0 1 2243 10436 2120 $Rewrd_4 -1 0 1 1281 12799 719 $Med_Fee -1 0 1 1788 10983 2028 $Low_Fee -1 0 1 1304 12771 724 $Bank_B -1 0 1 2688 10129 1982 $Out_State -1 0 1 1497 12819 483 $Med_Rebate -1 0 1 1902 10875 2022 $High_Rebate -1 0 1 1776 12777 246 $High_CredLine -1 0 1 2392 12287 120 $Long_Grace -1 0 2133 12666 > cat(" means of Demographic Variables",fill=TRUE) means of Demographic Variables > mat=apply(as.matrix(bank$demo[,2:3]),2,mean) > print(mat) age income 46.3055 46.1945 > > ## example of processing for use with rhierBinLogit > ## > if(nchar(Sys.getenv("LONG_TEST")) != 0) + { + choiceAtt=bank$choiceAtt + Z=bank$demo + + ## center demo data so that mean of random-effects + ## distribution can be interpretted as the average respondents + + Z[,1]=rep(1,nrow(Z)) + Z[,2]=Z[,2]-mean(Z[,2]) + Z[,3]=Z[,3]-mean(Z[,3]) + Z[,4]=Z[,4]-mean(Z[,4]) + Z=as.matrix(Z) + + hh=levels(factor(choiceAtt$id)) + nhh=length(hh) + lgtdata=NULL + for (i in 1:nhh) { + y=choiceAtt[choiceAtt[,1]==hh[i],2] + nobs=length(y) + X=as.matrix(choiceAtt[choiceAtt[,1]==hh[i],c(3:16)]) + lgtdata[[i]]=list(y=y,X=X) + } + + cat("Finished Reading data",fill=TRUE) + fsh() + + Data=list(lgtdata=lgtdata,Z=Z) + 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) + } > > > > > cleanEx(); ..nameEx <- "breg" > > ### * breg > > flush(stderr()); flush(stdout()) > > ### Name: breg > ### Title: Posterior Draws from a Univariate Regression with Unit Error > ### Variance > ### Aliases: breg > ### Keywords: models regression distribution > > ### ** Examples > > ## > > if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=1000} else {R=10} > > ## simulate data > set.seed(66) > n=100 > X=cbind(rep(1,n),runif(n)); beta=c(1,2) > y=X%*%beta+rnorm(n) > ## > ## set prior > A=diag(c(.05,.05)); betabar=c(0,0) > ## > ## make draws from posterior > betadraw=matrix(double(R*2),ncol=2) > for (rep in 1:R) {betadraw[rep,]=breg(y,X,betabar,A)} > ## > ## summarize draws > mat=apply(betadraw,2,quantile,probs=c(.01,.05,.5,.95,.99)) > mat=rbind(beta,mat); rownames(mat)[1]="beta"; print(mat) [,1] [,2] beta 1.0000000 2.000000 1% 0.5651933 2.010858 5% 0.6446456 2.204045 50% 0.8463473 2.571649 95% 0.9939671 2.784626 99% 1.0491400 2.795733 > > > > cleanEx(); ..nameEx <- "cgetC" > > ### * cgetC > > flush(stderr()); flush(stdout()) > > ### Name: cgetC > ### Title: Obtain A List of Cut-offs for Scale Usage Problems > ### Aliases: cgetC > ### Keywords: utilities > > ### ** Examples > > ## > cgetC(.1,10) [1] -1000.000000 2.537352 2.811348 3.285343 3.959338 [6] 4.833333 5.907329 7.181324 8.655319 10.329314 [11] 1000.000000 > > > > cleanEx(); ..nameEx <- "cheese" > > ### * cheese > > flush(stderr()); flush(stdout()) > > ### Name: cheese > ### Title: Sliced Cheese Data > ### Aliases: cheese > ### Keywords: datasets > > ### ** Examples > > data(cheese) > cat(" Quantiles of the Variables ",fill=TRUE) Quantiles of the Variables > mat=apply(as.matrix(cheese[,2:4]),2,quantile) > print(mat) VOLUME DISP PRICE 0% 231.0 0.00000000 1.319907 25% 1989.5 0.00000000 2.457262 50% 3408.0 0.04736842 2.703250 75% 5519.5 0.16600000 3.203279 100% 148109.0 1.00000000 4.641757 > > ## > ## example of processing for use with rhierLinearModel > ## > if(nchar(Sys.getenv("LONG_TEST")) != 0) + { + + retailer=levels(cheese$RETAILER) + nreg=length(retailer) + nvar=3 + regdata=NULL + for (reg in 1:nreg) { + y=log(cheese$VOLUME[cheese$RETAILER==retailer[reg]]) + iota=c(rep(1,length(y))) + X=cbind(iota,cheese$DISP[cheese$RETAILER==retailer[reg]], + log(cheese$PRICE[cheese$RETAILER==retailer[reg]])) + regdata[[reg]]=list(y=y,X=X) + } + Z=matrix(c(rep(1,nreg)),ncol=1) + nz=ncol(Z) + ## + ## run each individual regression and store results + ## + lscoef=matrix(double(nreg*nvar),ncol=nvar) + for (reg in 1:nreg) { + coef=lsfit(regdata[[reg]]$X,regdata[[reg]]$y,intercept=FALSE)$coef + if (var(regdata[[reg]]$X[,2])==0) { lscoef[reg,1]=coef[1]; lscoef[reg,3]=coef[2]} + else {lscoef[reg,]=coef } + } + + 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) + + 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)) + } + } + + } > > > > graphics::par(get("par.postscript", env = .CheckExEnv)) Warning: calling par(new=) with no plot > cleanEx(); ..nameEx <- "condMom" > > ### * condMom > > flush(stderr()); flush(stdout()) > > ### Name: condMom > ### Title: Computes Conditional Mean/Var of One Element of MVN given All > ### Others > ### Aliases: condMom > ### Keywords: distribution > > ### ** Examples > > ## > sig=matrix(c(1,.5,.5,.5,1,.5,.5,.5,1),ncol=3) > sigi=chol2inv(chol(sig)) > mu=c(1,2,3) > x=c(1,1,1) > condMom(x,mu,sigi,2) $cmean [1] 1.333333 $cvar [1] 0.6666667 > > > > > cleanEx(); ..nameEx <- "createX" > > ### * createX > > flush(stderr()); flush(stdout()) > > ### Name: createX > ### Title: Create X Matrix for Use in Multinomial Logit and Probit Routines > ### Aliases: createX > ### Keywords: array utilities > > ### ** Examples > > na=2; nd=1; p=3 > vec=c(1,1.5,.5,2,3,1,3,4.5,1.5) > Xa=matrix(vec,byrow=TRUE,ncol=3) > Xa=cbind(Xa,-Xa) > Xd=matrix(c(-1,-2,-3),ncol=1) > createX(p=p,na=na,nd=nd,Xa=Xa,Xd=Xd) [,1] [,2] [,3] [,4] [,5] [,6] [1,] 1 0 -1 0 1.0 -1.0 [2,] 0 1 0 -1 1.5 -1.5 [3,] 0 0 0 0 0.5 -0.5 [4,] 1 0 -2 0 2.0 -2.0 [5,] 0 1 0 -2 3.0 -3.0 [6,] 0 0 0 0 1.0 -1.0 [7,] 1 0 -3 0 3.0 -3.0 [8,] 0 1 0 -3 4.5 -4.5 [9,] 0 0 0 0 1.5 -1.5 > createX(p=p,na=na,nd=nd,Xa=Xa,Xd=Xd,base=1) [,1] [,2] [,3] [,4] [,5] [,6] [1,] 0 0 0 0 1.0 -1.0 [2,] 1 0 -1 0 1.5 -1.5 [3,] 0 1 0 -1 0.5 -0.5 [4,] 0 0 0 0 2.0 -2.0 [5,] 1 0 -2 0 3.0 -3.0 [6,] 0 1 0 -2 1.0 -1.0 [7,] 0 0 0 0 3.0 -3.0 [8,] 1 0 -3 0 4.5 -4.5 [9,] 0 1 0 -3 1.5 -1.5 > createX(p=p,na=na,nd=nd,Xa=Xa,Xd=Xd,DIFF=TRUE) [,1] [,2] [,3] [,4] [,5] [,6] [1,] 1 0 -1 0 0.5 -0.5 [2,] 0 1 0 -1 1.0 -1.0 [3,] 1 0 -2 0 1.0 -1.0 [4,] 0 1 0 -2 2.0 -2.0 [5,] 1 0 -3 0 1.5 -1.5 [6,] 0 1 0 -3 3.0 -3.0 > createX(p=p,na=na,nd=nd,Xa=Xa,Xd=Xd,DIFF=TRUE,base=2) [,1] [,2] [,3] [,4] [,5] [,6] [1,] 1 0 -1 0 -0.5 0.5 [2,] 0 1 0 -1 -1.0 1.0 [3,] 1 0 -2 0 -1.0 1.0 [4,] 0 1 0 -2 -2.0 2.0 [5,] 1 0 -3 0 -1.5 1.5 [6,] 0 1 0 -3 -3.0 3.0 > createX(p=p,na=na,nd=NULL,Xa=Xa,Xd=NULL) [,1] [,2] [,3] [,4] [1,] 1 0 1.0 -1.0 [2,] 0 1 1.5 -1.5 [3,] 0 0 0.5 -0.5 [4,] 1 0 2.0 -2.0 [5,] 0 1 3.0 -3.0 [6,] 0 0 1.0 -1.0 [7,] 1 0 3.0 -3.0 [8,] 0 1 4.5 -4.5 [9,] 0 0 1.5 -1.5 > createX(p=p,na=NULL,nd=nd,Xa=NULL,Xd=Xd) [,1] [,2] [,3] [,4] [1,] 1 0 -1 0 [2,] 0 1 0 -1 [3,] 0 0 0 0 [4,] 1 0 -2 0 [5,] 0 1 0 -2 [6,] 0 0 0 0 [7,] 1 0 -3 0 [8,] 0 1 0 -3 [9,] 0 0 0 0 > > > > cleanEx(); ..nameEx <- "customerSat" > > ### * customerSat > > flush(stderr()); flush(stdout()) > > ### Name: customerSat > ### Title: Customer Satifaction Data > ### Aliases: customerSat > ### Keywords: datasets > > ### ** Examples > > data(customerSat) > apply(as.matrix(customerSat),2,table) q1 q2 q3 q4 q5 q6 q7 q8 q9 q10 1 99 163 140 180 129 100 50 47 34 36 2 115 89 73 135 91 107 49 34 20 26 3 111 111 88 141 97 120 62 41 28 34 4 127 110 104 140 100 116 66 33 33 46 5 328 368 324 345 343 346 240 237 182 177 6 174 172 151 151 175 181 131 132 92 120 7 222 201 200 183 225 227 213 232 201 194 8 327 304 350 273 309 305 375 413 418 403 9 118 102 139 90 116 106 197 192 245 233 10 190 191 242 173 226 203 428 450 558 542 > > > > cleanEx(); ..nameEx <- "detailing" > > ### * detailing > > flush(stderr()); flush(stdout()) > > ### Name: detailing > ### Title: Physician Detailing Data from Manchanda et al (2004) > ### Aliases: detailing > ### Keywords: datasets > > ### ** Examples > > data(detailing) > cat(" table of Counts Dep Var", fill=TRUE) table of Counts Dep Var > print(table(detailing$counts[,2])) 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 5558 2822 2620 2130 1798 1387 1154 867 654 540 490 360 290 255 235 193 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 161 146 147 111 98 94 69 82 57 51 59 49 51 42 38 36 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 32 25 19 22 30 17 11 18 22 12 10 6 6 10 18 9 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 2 6 9 11 7 5 4 4 2 1 1 3 1 4 1 4 64 65 66 67 71 72 74 76 78 79 88 93 96 6 1 2 1 2 2 1 2 2 2 1 1 1 > cat(" means of Demographic Variables",fill=TRUE) means of Demographic Variables > mat=apply(as.matrix(detailing$demo[,2:4]),2,mean) > print(mat) generalphys specialist mean_samples 0.6010000 0.1850000 0.5639783 > > ## > ## example of processing for use with rhierNegbinRw > ## > if(nchar(Sys.getenv("LONG_TEST")) != 0) + { + + data(detailing) + counts = detailing$counts + Z = detailing$demo + + # Construct the Z matrix + Z[,1] = 1 + Z[,2]=Z[,2]-mean(Z[,2]) + Z[,3]=Z[,3]-mean(Z[,3]) + Z[,4]=Z[,4]-mean(Z[,4]) + Z=as.matrix(Z) + id=levels(factor(counts$id)) + nreg=length(id) + nobs = nrow(counts$id) + + regdata=NULL + for (i in 1:nreg) { + X = counts[counts[,1] == id[i],c(3:4)] + X = cbind(rep(1,nrow(X)),X) + y = counts[counts[,1] == id[i],2] + X = as.matrix(X) + regdata[[i]]=list(X=X, y=y) + } + nvar=ncol(X) # Number of X variables + nz=ncol(Z) # Number of Z variables + rm(detailing,counts) + cat("Finished Reading data",fill=TRUE) + fsh() + + Data = list(regdata=regdata, Z=Z) + deltabar = matrix(rep(0,nvar*nz),nrow=nz) + Vdelta = 0.01 * diag(nz) + 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) + + 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) + + # 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(" 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(" alphadraws ",fill=TRUE) + mat=apply(matrix(out$alphadraw),2,quantile,probs=c(.01,.05,.5,.95,.99)) + print(mat) + } > > > > cleanEx(); ..nameEx <- "ghkvec" > > ### * ghkvec > > flush(stderr()); flush(stdout()) > > ### Name: ghkvec > ### Title: Compute GHK approximation to Multivariate Normal Integrals > ### Aliases: ghkvec > ### Keywords: distribution > > ### ** Examples > > ## > > Sigma=matrix(c(1,.5,.5,1),ncol=2) > L=t(chol(Sigma)) > trunpt=c(0,0,1,1) > above=c(1,1) > ghkvec(L,trunpt,above,100) [1] 0.3257788 0.7497757 > > > > cleanEx(); ..nameEx <- "llmnl" > > ### * llmnl > > flush(stderr()); flush(stdout()) > > ### Name: llmnl > ### Title: Evaluate Log Likelihood for Multinomial Logit Model > ### Aliases: llmnl > ### Keywords: models > > ### ** Examples > > ## > ## Not run: ll=llmnl(y,X,beta) > > > > cleanEx(); ..nameEx <- "llmnp" > > ### * llmnp > > flush(stderr()); flush(stdout()) > > ### Name: llmnp > ### Title: Evaluate Log Likelihood for Multinomial Probit Model > ### Aliases: llmnp > ### Keywords: models > > ### ** Examples > > ## > ## Not run: ll=llmnp(X,y,beta,Sigma,r) > > > > cleanEx(); ..nameEx <- "llnhlogit" > > ### * llnhlogit > > flush(stderr()); flush(stdout()) > > ### Name: llnhlogit > ### Title: Evaluate Log Likelihood for non-homothetic Logit Model > ### Aliases: llnhlogit > ### Keywords: models > > ### ** Examples > > ## > ## Not run: ll=llnhlogit(theta,choice,lnprices,Xexpend) > > > > cleanEx(); ..nameEx <- "lndIChisq" > > ### * lndIChisq > > flush(stderr()); flush(stdout()) > > ### Name: lndIChisq > ### Title: Compute Log of Inverted Chi-Squared Density > ### Aliases: lndIChisq > ### Keywords: distribution > > ### ** Examples > > ## > lndIChisq(3,1,2) [1] -1.753888 > > > > cleanEx(); ..nameEx <- "lndIWishart" > > ### * lndIWishart > > flush(stderr()); flush(stdout()) > > ### Name: lndIWishart > ### Title: Compute Log of Inverted Wishart Density > ### Aliases: lndIWishart > ### Keywords: distribution > > ### ** Examples > > ## > lndIWishart(5,diag(3),(diag(3)+.5)) [1] -12.40291 > > > > cleanEx(); ..nameEx <- "lndMvn" > > ### * lndMvn > > flush(stderr()); flush(stdout()) > > ### Name: lndMvn > ### Title: Compute Log of Multivariate Normal Density > ### Aliases: lndMvn > ### Keywords: distribution > > ### ** Examples > > ## > Sigma=matrix(c(1,.5,.5,1),ncol=2) > lndMvn(x=c(rep(0,2)),mu=c(rep(0,2)),rooti=backsolve(chol(Sigma),diag(2))) [,1] [1,] 0.1438410 > > > > cleanEx(); ..nameEx <- "lndMvst" > > ### * lndMvst > > flush(stderr()); flush(stdout()) > > ### Name: lndMvst > ### Title: Compute Log of Multivariate Student-t Density > ### Aliases: lndMvst > ### Keywords: distribution > > ### ** Examples > > ## > Sigma=matrix(c(1,.5,.5,1),ncol=2) > lndMvst(x=c(rep(0,2)),nu=4,mu=c(rep(0,2)),rooti=backsolve(chol(Sigma),diag(2))) [,1] [1,] -4.015042 > > > > cleanEx(); ..nameEx <- "margarine" > > ### * margarine > > flush(stderr()); flush(stdout()) > > ### Name: margarine > ### Title: Household Panel Data on Margarine Purchases > ### Aliases: margarine > ### Keywords: datasets > > ### ** Examples > > data(margarine) > cat(" Table of Choice Variable ",fill=TRUE) Table of Choice Variable > print(table(margarine$choicePrice[,2])) 1 2 3 4 5 6 7 8 9 10 1766 699 243 593 315 74 319 203 225 33 > cat(" Means of Prices",fill=TRUE) Means of Prices > mat=apply(as.matrix(margarine$choicePrice[,3:12]),2,mean) > print(mat) PPk_Stk PBB_Stk PFl_Stk PHse_Stk PGen_Stk PImp_Stk PSS_Tub PPk_Tub 0.5184362 0.5432103 1.0150201 0.4371477 0.3452819 0.7807785 0.8250895 1.0774094 PFl_Tub PHse_Tub 1.1893758 0.5686734 > cat(" Quantiles of Demographic Variables",fill=TRUE) Quantiles of Demographic Variables > mat=apply(as.matrix(margarine$demos[,2:8]),2,quantile) > print(mat) Income Fs3_4 Fs5. Fam_Size college whtcollar retired 0% 2.5 0 0 1 0 0 0 25% 17.5 0 0 2 0 0 0 50% 22.5 0 0 3 0 1 0 75% 32.5 1 0 4 1 1 0 100% 130.0 1 1 8 1 1 1 > > ## > ## example of processing for use with rhierMnlRwMixture > ## > if(nchar(Sys.getenv("LONG_TEST")) != 0) + { + select= c(1:5,7) ## select brands + chPr=as.matrix(margarine$choicePrice) + ## make sure to log prices + chPr=cbind(chPr[,1],chPr[,2],log(chPr[,2+select])) + demos=as.matrix(margarine$demos[,c(1,2,5)]) + + ## remove obs for other alts + chPr=chPr[chPr[,2] <= 7,] + chPr=chPr[chPr[,2] != 6,] + + ## recode choice + chPr[chPr[,2] == 7,2]=6 + + hhidl=levels(as.factor(chPr[,1])) + lgtdata=NULL + nlgt=length(hhidl) + p=length(select) ## number of choice alts + ind=1 + for (i in 1:nlgt) { + nobs=sum(chPr[,1]==hhidl[i]) + if(nobs >=5) { + data=chPr[chPr[,1]==hhidl[i],] + y=data[,2] + names(y)=NULL + X=createX(p=p,na=1,Xa=data[,3:8],nd=NULL,Xd=NULL,INT=TRUE,base=1) + lgtdata[[ind]]=list(y=y,X=X,hhid=hhidl[i]); ind=ind+1 + } + } + nlgt=length(lgtdata) + ## + ## now extract demos corresponding to hhs in lgtdata + ## + Z=NULL + nlgt=length(lgtdata) + for(i in 1:nlgt){ + Z=rbind(Z,demos[demos[,1]==lgtdata[[i]]$hhid,2:3]) + } + ## + ## take log of income and family size and demean + ## + Z=log(Z) + Z[,1]=Z[,1]-mean(Z[,1]) + Z[,2]=Z[,2]-mean(Z[,2]) + + keep=5 + 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) + + } > > > > cleanEx(); ..nameEx <- "mnlHess" > > ### * mnlHess > > flush(stderr()); flush(stdout()) > > ### Name: mnlHess > ### Title: Computes -Expected Hessian for Multinomial Logit > ### Aliases: mnlHess > ### Keywords: models > > ### ** Examples > > ## > ## Not run: mnlHess(y,X,beta) > > > > cleanEx(); ..nameEx <- "nmat" > > ### * nmat > > flush(stderr()); flush(stdout()) > > ### Name: nmat > ### Title: Convert Covariance Matrix to a Correlation Matrix > ### Aliases: nmat > ### Keywords: utilities array > > ### ** Examples > > ## > set.seed(66) > X=matrix(rnorm(200,4),ncol=2) > Varmat=var(X) > nmat(as.vector(Varmat)) [1] 1.00000000 -0.03206572 -0.03206572 1.00000000 > > > > cleanEx(); ..nameEx <- "numEff" > > ### * numEff > > flush(stderr()); flush(stdout()) > > ### Name: numEff > ### Title: Compute Numerical Standard Error and Relative Numerical > ### Efficiency > ### Aliases: numEff > ### Keywords: ts utilities > > ### ** Examples > > numEff(rnorm(1000),m=20) $stderr [,1] [1,] 0.02889749 $f [,1] [1,] 0.7796689 $m [1] 20 > numEff(rnorm(1000)) $stderr [,1] [1,] 0.03698488 $f [,1] [1,] 1.264728 $m [1] 44 > > > > cleanEx(); ..nameEx <- "rbprobitGibbs" > > ### * rbprobitGibbs > > flush(stderr()); flush(stdout()) > > ### Name: rbprobitGibbs > ### Title: Gibbs Sampler (Albert and Chib) for Binary Probit > ### Aliases: rbprobitGibbs > ### Keywords: models > > ### ** Examples > > ## > ## rbprobitGibbs example > ## > if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=2000} else {R=10} > > set.seed(66) > simbprobit= + function(X,beta) { + ## function to simulate from binary probit including x variable + y=ifelse((X%*%beta+rnorm(nrow(X)))<0,0,1) + list(X=X,y=y,beta=beta) + } > > nobs=200 > X=cbind(rep(1,nobs),runif(nobs),runif(nobs)) > 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) > > out=rbprobitGibbs(Data=Data,Mcmc=Mcmc) Starting Gibbs Sampler for Binary Probit Model Prior Parms: betabar [1] 0 0 0 A [,1] [,2] [,3] [1,] 0.01 0.00 0.00 [2,] 0.00 0.01 0.00 [3,] 0.00 0.00 0.01 MCMC parms: R= 10 keep= 1 MCMC Iteration (est time to end - min) Total Time Elapsed: 0 > > cat(" Betadraws ",fill=TRUE) Betadraws > mat=apply(out$betadraw,2,quantile,probs=c(.01,.05,.5,.95,.99)) > mat=rbind(beta,mat); rownames(mat)[1]="beta"; print(mat) [,1] [,2] [,3] beta 0.00000000 1.0000000 -1.0000000 1% 0.04489536 0.5207186 -1.8111980 5% 0.07030921 0.5302811 -1.6453382 50% 0.37172881 0.6874956 -0.9858278 95% 0.70089020 0.8834252 -0.7156283 99% 0.85166245 0.9304226 -0.7022984 > > > > cleanEx(); ..nameEx <- "rdirichlet" > > ### * rdirichlet > > flush(stderr()); flush(stdout()) > > ### Name: rdirichlet > ### Title: Draw From Dirichlet Distribution > ### Aliases: rdirichlet > ### Keywords: distribution > > ### ** Examples > > ## > set.seed(66) > rdirichlet(c(rep(3,5))) [1] 0.3603397 0.1367107 0.1534761 0.1056937 0.2437799 > > > > cleanEx(); ..nameEx <- "rhierBinLogit" > > ### * rhierBinLogit > > flush(stderr()); flush(stdout()) > > ### Name: rhierBinLogit > ### Title: MCMC Algorithm for Hierachical Binary Logit > ### Aliases: rhierBinLogit > ### Keywords: models > > ### ** Examples > > ## > if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=10000} else {R=10} > > set.seed(66) > nvar=5 ## number of coefficients > nlgt=1000 ## number of cross-sectional units > nobs=10 ## number of observations per unit > nz=2 ## number of regressors in mixing distribution > > ## set hyper-parameters > ## B=ZDelta + U > > Z=matrix(c(rep(1,nlgt),runif(nlgt,min=-1,max=1)),nrow=nlgt,ncol=nz) > Delta=matrix(c(-2,-1,0,1,2,-1,1,-.5,.5,0),nrow=nz,ncol=nvar) > iota=matrix(1,nrow=nvar,ncol=1) > Vbeta=diag(nvar)+.5*iota%*%t(iota) > > ## simulate data > lgtdata=NULL > > for (i in 1:nlgt) + { beta=t(Delta)%*%Z[i,]+as.vector(t(chol(Vbeta))%*%rnorm(nvar)) + X=matrix(runif(nobs*nvar),nrow=nobs,ncol=nvar) + prob=exp(X%*%beta)/(1+exp(X%*%beta)) + unif=runif(nobs,0,1) + y=ifelse(unif > Data=list(Dat=lgtdata,Demo=Z) > out=rhierBinLogit(Data=list(lgtdata=lgtdata,Z=Z),Mcmc=list(R=R)) Table of Y values pooled over all units ypooled 0 1 4097 5903 Attempting MCMC Inference for Hierarchical Binary Logit: 5 variables in X 2 variables in Z for 1000 cross-sectional units Prior Parms: nu = 8 V [,1] [,2] [,3] [,4] [,5] [1,] 8 0 0 0 0 [2,] 0 8 0 0 0 [3,] 0 0 8 0 0 [4,] 0 0 0 8 0 [5,] 0 0 0 0 8 Deltabar [,1] [,2] [,3] [,4] [,5] [1,] 0 0 0 0 0 [2,] 0 0 0 0 0 ADelta [,1] [,2] [1,] 0.01 0.00 [2,] 0.00 0.01 MCMC Parms: sbeta= 0.2 R= 10 keep= 1 MCMC Iteration (est time to end - min) Total Time Elapsed: 0.08 > > cat(" Deltadraws ",fill=TRUE) Deltadraws > 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) [,1] [,2] [,3] [,4] [,5] [,6] delta -2.000000000 -1.00000000 0.000000000 1.00000000 2.000000000 -1.000000000 1% 0.002070074 -0.06636228 0.007852401 -0.03617905 0.009220025 -0.056725649 5% 0.004659883 -0.06020435 0.008781428 -0.03270958 0.012528509 -0.055186793 50% 0.027635284 -0.02702937 0.021459189 -0.01064371 0.028847655 -0.021631955 95% 0.043856470 0.01473469 0.041457829 0.02009469 0.048052871 0.007830176 99% 0.046143741 0.02224859 0.041874301 0.02593557 0.049029208 0.012580049 [,7] [,8] [,9] [,10] delta 1.000000000 -0.500000000 0.500000000 0.000000000 1% 0.005939521 -0.065119882 -0.001054080 -0.064228871 5% 0.008839744 -0.063971604 0.001906039 -0.059556319 50% 0.029647060 -0.042678671 0.024451108 -0.028568448 95% 0.064944083 -0.005928613 0.050872125 -0.007310795 99% 0.068957594 -0.000206114 0.054797323 -0.002162184 > cat(" Vbetadraws ",fill=TRUE) Vbetadraws > 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) [,1] [,2] [,3] [,4] [,5] Vbeta 1.5000000 0.500000000 0.500000000 0.5000000000 0.500000000 1% 0.1067367 -0.009289902 0.001014991 -0.0139946484 -0.006065785 5% 0.1084473 -0.008479434 0.001304096 -0.0123409794 -0.005396924 50% 0.1276727 -0.001398834 0.007783606 -0.0056337836 0.002961030 95% 0.1367680 0.008868526 0.018851885 -0.0000958734 0.007024520 99% 0.1371433 0.009677800 0.019909834 0.0022777784 0.007190185 [,6] [,7] [,8] [,9] [,10] [,11] Vbeta 0.500000000 1.5000000 0.500000000 0.5000000000 0.500000000 0.500000000 1% -0.009289902 0.1015031 -0.001356852 -0.0131672385 0.001825738 0.001014991 5% -0.008479434 0.1019631 -0.001008748 -0.0117350614 0.002384750 0.001304096 50% -0.001398834 0.1058310 0.003634585 -0.0063620617 0.005044709 0.007783606 95% 0.008868526 0.1181146 0.010698996 -0.0010235724 0.009873027 0.018851885 99% 0.009677800 0.1214456 0.011262689 -0.0009940816 0.011032195 0.019909834 [,12] [,13] [,14] [,15] [,16] Vbeta 0.500000000 1.5000000 0.500000000 0.500000000 0.5000000000 1% -0.001356852 0.1029966 0.001827070 0.002523644 -0.0139946484 5% -0.001008748 0.1035354 0.003239515 0.002862334 -0.0123409794 50% 0.003634585 0.1135662 0.007285149 0.006237820 -0.0056337836 95% 0.010698996 0.1235704 0.018672874 0.013261132 -0.0000958734 99% 0.011262689 0.1249173 0.019941812 0.013547746 0.0022777784 [,17] [,18] [,19] [,20] [,21] [,22] Vbeta 0.5000000000 0.500000000 1.5000000 0.500000000 0.500000000 0.500000000 1% -0.0131672385 0.001827070 0.1076981 0.002369389 -0.006065785 0.001825738 5% -0.0117350614 0.003239515 0.1088122 0.004004557 -0.005396924 0.002384750 50% -0.0063620617 0.007285149 0.1343052 0.011601002 0.002961030 0.005044709 95% -0.0010235724 0.018672874 0.1557306 0.021620214 0.007024520 0.009873027 99% -0.0009940816 0.019941812 0.1586929 0.022659191 0.007190185 0.011032195 [,23] [,24] [,25] Vbeta 0.500000000 0.500000000 1.5000000 1% 0.002523644 0.002369389 0.1018685 5% 0.002862334 0.004004557 0.1031904 50% 0.006237820 0.011601002 0.1171665 95% 0.013261132 0.021620214 0.1375157 99% 0.013547746 0.022659191 0.1378044 > > 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) + } > > > > > graphics::par(get("par.postscript", env = .CheckExEnv)) Warning: calling par(new=) with no plot > cleanEx(); ..nameEx <- "rhierLinearModel" > > ### * rhierLinearModel > > flush(stderr()); flush(stdout()) > > ### Name: rhierLinearModel > ### Title: Gibbs Sampler for Hierarchical Linear Model > ### Aliases: rhierLinearModel > ### Keywords: regression > > ### ** Examples > > ## > if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=2000} else {R=10} > > nreg=100; nobs=100; nvar=3 > Vbeta=matrix(c(1,.5,0,.5,2,.7,0,.7,1),ncol=3) > Z=cbind(c(rep(1,nreg)),3*runif(nreg)); Z[,2]=Z[,2]-mean(Z[,2]) > nz=ncol(Z) > Delta=matrix(c(1,-1,2,0,1,0),ncol=2) > Delta=t(Delta) # first row of Delta is means of betas > Beta=matrix(rnorm(nreg*nvar),nrow=nreg)%*%chol(Vbeta)+Z%*%Delta > tau=.1 > iota=c(rep(1,nobs)) > 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) > > 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) Starting Gibbs Sampler for Linear Hierarchical Model 100 Regressions 2 Variables in Z (if 1, then only intercept) Prior Parms: Deltabar [,1] [,2] [,3] [1,] 0 0 0 [2,] 0 0 0 A [,1] [,2] [1,] 1 0 [2,] 0 1 nu.e (d.f. parm for regression error variances)= 3 Vbeta ~ IW(nu,V) nu = 6 V [,1] [,2] [,3] [1,] 6 0 0 [2,] 0 6 0 [3,] 0 0 6 MCMC parms: R= 10 keep= 1 MCMC Iteration (est time to end - min) Total Time Elapsed: 0.02 > > cat(" Deltadraws ",fill=TRUE) Deltadraws > 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) [,1] [,2] [,3] [,4] [,5] [,6] delta 1.0000000 0.00000000 -1.0000000 1.0000000 2.000000 0.00000000 1% 0.9049336 -0.08587196 -1.2947333 0.7355300 1.929383 -0.34261218 5% 0.9104958 -0.03508422 -1.2071392 0.7401179 1.967624 -0.33977007 50% 0.9701896 0.13157680 -0.9558319 0.9356783 2.080607 -0.18673750 95% 1.0301048 0.41086563 -0.6876690 1.2454277 2.177079 -0.06597003 99% 1.0337341 0.43523270 -0.6514537 1.2502983 2.181942 -0.05042893 > cat(" Vbetadraws ",fill=TRUE) Vbetadraws > 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) [,1] [,2] [,3] [,4] [,5] [,6] [,7] Vbeta 1.0000000 0.5000000 0.00000000 0.5000000 2.000000 0.7000000 0.00000000 1% 0.7286528 0.3215465 -0.09218703 0.3215465 1.533121 0.3957224 -0.09218703 5% 0.7626521 0.3442439 -0.06880559 0.3442439 1.786221 0.4840415 -0.06880559 50% 0.9501390 0.4907790 0.04771368 0.4907790 2.252869 0.7542128 0.04771368 95% 1.3012112 0.8826273 0.22043954 0.8826273 2.812445 0.9024790 0.22043954 99% 1.4628746 1.0464061 0.27552917 1.0464061 2.944918 0.9264810 0.27552917 [,8] [,9] Vbeta 0.7000000 1.0000000 1% 0.3957224 0.8195369 5% 0.4840415 0.8419077 50% 0.7542128 0.9657227 95% 0.9024790 1.1108132 99% 0.9264810 1.1255840 > > > > > cleanEx(); ..nameEx <- "rhierMnlRwMixture" > > ### * rhierMnlRwMixture > > flush(stderr()); flush(stdout()) > > ### Name: rhierMnlRwMixture > ### Title: MCMC Algorithm for Hierarchical Multinomial Logit with Mixture > ### of Normals Heterogeneity > ### Aliases: rhierMnlRwMixture > ### Keywords: models > > ### ** Examples > > ## > if(nchar(Sys.getenv("LONG_TEST")) != 0) + { + + set.seed(66) + p=3 # num of choice alterns + ncoef=3 + nlgt=300 # num of cross sectional units + nz=2 + Z=matrix(runif(nz*nlgt),ncol=nz) + Z=t(t(Z)-apply(Z,2,mean)) # demean Z + ncomp=3 # no of mixture components + Delta=matrix(c(1,0,1,0,1,2),ncol=2) + comps=NULL + comps[[1]]=list(mu=c(0,-1,-2),rooti=diag(rep(1,3))) + 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) + + ## 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) } + outa=simmnlwX(ni[i],X,betai) + simlgtdata[[i]]=list(y=outa$y,X=X,beta=betai) + } + + ## plot betas + if(0){ + ## set if(1) above to produce plots + bmat=matrix(0,nlgt,ncoef) + for(i in 1:nlgt) {bmat[i,]=simlgtdata[[i]]$beta} + par(mfrow=c(ncoef,1)) + 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,ncoef) + + R=10000 + keep=5 + c=2 + s=2.93/sqrt(ncoef) + 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,c=c,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)) + + 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")} + } + } > > > > > graphics::par(get("par.postscript", env = .CheckExEnv)) Warning: calling par(new=) with no plot > cleanEx(); ..nameEx <- "rhierNegbinRw" > > ### * rhierNegbinRw > > flush(stderr()); flush(stdout()) > > ### Name: rhierNegbinRw > ### Title: MCMC Algorithm for Negative Binomial Regression > ### Aliases: rhierNegbinRw > ### Keywords: models > > ### ** Examples > > ## > if(nchar(Sys.getenv("LONG_TEST")) != 0) + { + ## + set.seed(66) + simnegbin = + function(X, beta, alpha) { + # Simulate from the Negative Binomial Regression + lambda = exp(X %*% beta) + y=NULL + for (j in 1:length(lambda)) + y = c(y,rnbinom(1,mu = lambda[j],size = alpha)) + return(y) + } + + nreg = 100 # Number of cross sectional units + T = 50 # Number of observations per unit + nobs = nreg*T + nvar=2 # Number of X variables + 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)) + alpha = 5 + Vbeta = rbind(c(0.1,0),c(0,0.1)) + + # Construct the regdata (containing X) + simnegbindata = NULL + for (i in 1:nreg) { + betai = as.vector(Z[i,]%*%Delta) + chol(Vbeta)%*%rnorm(nvar) + X = cbind(rep(1,T),rnorm(T,mean=2,sd=0.25)) + simnegbindata[[i]] = list(y=simnegbin(X,betai,alpha), X=X,beta=betai) + } + + 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) + + 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) + + # 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(" 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) + } > > > > cleanEx(); ..nameEx <- "rivGibbs" > > ### * rivGibbs > > flush(stderr()); flush(stdout()) > > ### Name: rivGibbs > ### Title: Gibbs Sampler for Linear "IV" Model > ### Aliases: rivGibbs > ### Keywords: models > > ### ** Examples > > ## > if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=2000} else {R=10} > > set.seed(66) > simIV = function(delta,beta,Sigma,n,z,w,gamma) { + eps = matrix(rnorm(2*n),ncol=2) %*% chol(Sigma) + x = z %*% delta + eps[,1]; y = beta*x + eps[,2] + w%*%gamma + list(x=as.vector(x),y=as.vector(y)) } > n = 200 ; p=1 # number of instruments > z = cbind(rep(1,n),matrix(runif(n*p),ncol=p)) > w = matrix(1,n,1) > rho=.8 > 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) Starting Gibbs Sampler for Linear IV Model nobs= 200 ; 2 instruments; 1 included exog vars Prior Parms: mean of delta [1] 0 0 Adelta [,1] [,2] [1,] 0.01 0.00 [2,] 0.00 0.01 mean of beta/gamma [1] 0 0 Abeta/gamma [,1] [,2] [1,] 0.01 0.00 [2,] 0.00 0.01 Sigma Prior Parms nu= 3 V= [,1] [,2] [1,] 3 0 [2,] 0 3 MCMC parms: R= 10 keep= 1 MCMC Iteration (est time to end -min) > > cat(" deltadraws ",fill=TRUE) deltadraws > mat=apply(out$deltadraw,2,quantile,probs=c(.01,.05,.5,.95,.99)) > mat=rbind(delta,mat); rownames(mat)[1]="delta"; print(mat) [,1] [,2] delta 1.0000000 4.000000 1% 0.8284758 3.869874 5% 0.8392815 3.968293 50% 0.9489632 4.328563 95% 1.1376847 4.495194 99% 1.2062319 4.575587 > cat(" betadraws ",fill=TRUE) betadraws > 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) [,1] beta 0.5000000 1% 0.5127719 5% 0.5197926 50% 0.5836201 95% 0.7459027 99% 0.8165350 > cat(" Sigma draws",fill=TRUE) Sigma draws > 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) [,1] [,2] [,3] [,4] Sigma 1.0000000 0.8000000 0.8000000 1.0000000 1% 0.8575873 0.5040537 0.5040537 0.7118259 5% 0.8616333 0.5502728 0.5502728 0.7732516 50% 0.9597415 0.7582458 0.7582458 0.9547994 95% 1.1301939 0.8922998 0.8922998 1.0709312 99% 1.1640956 0.8932660 0.8932660 1.0935975 > > > > cleanEx(); ..nameEx <- "rmnlIndepMetrop" > > ### * rmnlIndepMetrop > > flush(stderr()); flush(stdout()) > > ### Name: rmnlIndepMetrop > ### Title: MCMC Algorithm for Multinomial Logit Model > ### Aliases: rmnlIndepMetrop > ### Keywords: models > > ### ** Examples > > ## > > 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) > 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) Starting Independence Metropolis Sampler for Multinomial Logit Model with 3 alternatives Prior Parms: betabar [1] 0 0 0 0 A [,1] [,2] [,3] [,4] [1,] 0.01 0.00 0.00 0.00 [2,] 0.00 0.01 0.00 0.00 [3,] 0.00 0.00 0.01 0.00 [4,] 0.00 0.00 0.00 0.01 MCMC parms: R= 10 keep= 1 nu (df for st candidates) = 6 MCMC Iteration (est time to end - min) Total Time Elapsed: 0 > cat(" Betadraws ",fill=TRUE) Betadraws > mat=apply(out$betadraw,2,quantile,probs=c(.01,.05,.5,.95,.99)) > mat=rbind(beta,mat); rownames(mat)[1]="beta"; print(mat) [,1] [,2] [,3] [,4] beta 1.0000000 -1.0000000 1.500000 0.5000000 1% 0.7957876 -1.6230397 1.040821 0.1729344 5% 0.8007537 -1.4267747 1.161966 0.1830847 50% 0.8774211 -0.9454719 1.570660 0.3347438 95% 1.3847157 -0.7088397 1.830373 0.5830628 99% 1.5128462 -0.6862219 1.837606 0.6593664 > > > > > cleanEx(); ..nameEx <- "rmnpGibbs" > > ### * rmnpGibbs > > flush(stderr()); flush(stdout()) > > ### Name: rmnpGibbs > ### Title: Gibbs Sampler for Multinomial Probit > ### Aliases: rmnpGibbs > ### Keywords: models > > ### ** Examples > > ## > if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=2000} else {R=10} > > set.seed(66) > p=3 > 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) > simout=simmnp(X,p,500,beta,Sigma) > > Data=list(p=p,y=simout$y,X=simout$X) > Mcmc=list(R=R,keep=1) > > out=rmnpGibbs(Mcmc=Mcmc,Data=Data) Table of y values y 1 2 3 58 321 121 Starting Gibbs Sampler for MNP 500 obs; 3 choice alternatives; 4 indep vars (including intercepts) 10 reps; keeping every 1 th draw Prior Parms: betabar [1] 0 0 0 0 A [,1] [,2] [,3] [,4] [1,] 0.01 0.00 0.00 0.00 [2,] 0.00 0.01 0.00 0.00 [3,] 0.00 0.00 0.01 0.00 [4,] 0.00 0.00 0.00 0.01 nu [1] 5 V [,1] [,2] [1,] 5 0 [2,] 0 5 MCMC Parms: R= 10 initial beta= 0 0 0 0 initial sigma= 1 0 0 1 MCMC Iteration (est time to end - min) > > cat(" Betadraws ",fill=TRUE) Betadraws > 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) [,1] [,2] [,3] [,4] beta -1.0000000 1.0000000 1.0000000 2.0000000 1% -1.2559544 0.2900262 0.3908498 0.8050299 5% -1.2498188 0.3057026 0.4312908 0.9638476 50% -1.0661117 0.6737744 0.8511478 1.8451683 95% -0.7603244 0.9211566 1.3277862 2.1507614 99% -0.7095687 0.9618609 1.3922750 2.1757275 > cat(" Sigmadraws ",fill=TRUE) Sigmadraws > 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) [,1] [,2] [,3] [,4] sigma 1 0.500000000 0.500000000 1.000000 1% 1 -0.004993161 -0.004993161 1.019450 5% 1 0.074007401 0.074007401 1.122840 50% 1 0.283975627 0.283975627 1.349785 95% 1 0.324487979 0.324487979 1.545948 99% 1 0.337346298 0.337346298 1.546357 > > > > > cleanEx(); ..nameEx <- "rmultireg" > > ### * rmultireg > > flush(stderr()); flush(stdout()) > > ### Name: rmultireg > ### Title: Draw from the Posterior of a Multivariate Regression > ### Aliases: rmultireg > ### Keywords: regression > > ### ** Examples > > ## > if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=2000} else {R=10} > > set.seed(66) > n=200 > m=2 > X=cbind(rep(1,n),runif(n)) > k=ncol(X) > B=matrix(c(1,2,-1,3),ncol=m) > Sigma=matrix(c(1,.5,.5,1),ncol=m); RSigma=chol(Sigma) > Y=X%*%B+matrix(rnorm(m*n),ncol=m)%*%RSigma > > betabar=rep(0,k*m);Bbar=matrix(betabar,ncol=m) > A=diag(rep(.01,k)) > nu=3; V=nu*diag(m) > > betadraw=matrix(double(R*k*m),ncol=k*m) > Sigmadraw=matrix(double(R*m*m),ncol=m*m) > for (rep in 1:R) + {out=rmultireg(Y,X,Bbar,A,nu,V);betadraw[rep,]=out$B + Sigmadraw[rep,]=out$Sigma} > > cat(" Betadraws ",fill=TRUE) Betadraws > mat=apply(betadraw,2,quantile,probs=c(.01,.05,.5,.95,.99)) > mat=rbind(as.vector(B),mat); rownames(mat)[1]="beta" > print(mat) [,1] [,2] [,3] [,4] beta 1.0000000 2.000000 -1.0000000 3.000000 1% 0.8935809 1.661419 -1.2027068 2.629172 5% 0.9188952 1.685108 -1.1443588 2.665079 50% 1.0258584 2.071904 -0.9345797 3.167549 95% 1.2922225 2.450720 -0.6758421 3.539084 99% 1.3162992 2.475153 -0.6736775 3.685967 > cat(" Sigma draws",fill=TRUE) Sigma draws > mat=apply(Sigmadraw,2,quantile,probs=c(.01,.05,.5,.95,.99)) > mat=rbind(as.vector(Sigma),mat); rownames(mat)[1]="Sigma" > print(mat) [,1] [,2] [,3] [,4] Sigma 1.0000000 0.5000000 0.5000000 1.0000000 1% 0.8671183 0.4611315 0.4611315 0.8108963 5% 0.8694217 0.4665847 0.4665847 0.8346087 50% 0.9615998 0.5242548 0.5242548 0.9472345 95% 1.0561278 0.6926645 0.6926645 1.2994628 99% 1.0700665 0.6974689 0.6974689 1.3800294 > > > > cleanEx(); ..nameEx <- "rmvpGibbs" > > ### * rmvpGibbs > > flush(stderr()); flush(stdout()) > > ### Name: rmvpGibbs > ### Title: Gibbs Sampler for Multivariate Probit > ### Aliases: rmvpGibbs > ### Keywords: models multivariate > > ### ** Examples > > ## > if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=2000} else {R=10} > > set.seed(66) > p=3 > n=500 > beta=c(-2,0,2) > 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 > 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) Table of y values y 0 1 752 748 Starting Gibbs Sampler for MVP 500 obs of 3 binary indicators; 3 indep vars (including intercepts) 10 reps; keeping every 1 th draw Prior Parms: betabar [1] 0 0 0 A [,1] [,2] [,3] [1,] 0.01 0.00 0.00 [2,] 0.00 0.01 0.00 [3,] 0.00 0.00 0.01 nu [1] 6 V [,1] [,2] [,3] [1,] 6 0 0 [2,] 0 6 0 [3,] 0 0 6 MCMC Parms: R= 10 initial beta= 0 0 0 initial sigma= [,1] [,2] [,3] [1,] 1 0 0 [2,] 0 1 0 [3,] 0 0 1 MCMC Iteration (est time to end - min) > > ind=seq(from=0,by=p,length=k) > inda=1:3 > ind=ind+inda > cat(" Betadraws ",fill=TRUE) Betadraws > 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) [,1] [,2] [,3] beta -2.000000 0.00000000 2.000000 1% -1.833383 -0.12208477 1.210622 5% -1.828588 -0.10633056 1.318589 50% -1.765531 -0.04514932 1.777502 95% -1.261021 0.03576498 2.010827 99% -1.125576 0.05506411 2.028375 > rdraw=matrix(double((R)*p*p),ncol=p*p) > rdraw=t(apply(out$sigmadraw,1,nmat)) > cat(" Draws of Correlation Matrix ",fill=TRUE) Draws of Correlation Matrix > mat=apply(rdraw,2,quantile,probs=c(.01,.05,.5,.95,.99)) > mat=rbind(as.vector(Sigma),mat); rownames(mat)[1]="Sigma"; print(mat) [,1] [,2] [,3] [,4] [,5] [,6] [,7] Sigma 1 0.50000000 0.500000000 0.50000000 1 0.5000000000 0.500000000 1% 1 0.04614673 -0.018068193 0.04614673 1 -0.0061935989 -0.018068193 5% 1 0.04722861 0.001076157 0.04722861 1 0.0006719889 0.001076157 50% 1 0.12833087 0.077087682 0.12833087 1 0.0639324283 0.077087682 95% 1 0.24912851 0.123232430 0.24912851 1 0.1684707015 0.123232430 99% 1 0.25951470 0.124050114 0.25951470 1 0.1728003932 0.124050114 [,8] [,9] Sigma 0.5000000000 1 1% -0.0061935989 1 5% 0.0006719889 1 50% 0.0639324283 1 95% 0.1684707015 1 99% 0.1728003932 1 > > > > > cleanEx(); ..nameEx <- "rmvst" > > ### * rmvst > > flush(stderr()); flush(stdout()) > > ### Name: rmvst > ### Title: Draw from Multivariate Student-t > ### Aliases: rmvst > ### Keywords: distribution > > ### ** Examples > > ## > set.seed(66) > rmvst(nu=5,mu=c(rep(0,2)),root=chol(matrix(c(2,1,1,2),ncol=2))) [,1] [1,] 3.201216 [2,] 1.859446 > > > > cleanEx(); ..nameEx <- "rnegbinRw" > > ### * rnegbinRw > > flush(stderr()); flush(stdout()) > > ### Name: rnegbinRw > ### Title: MCMC Algorithm for Negative Binomial Regression > ### Aliases: rnegbinRw > ### Keywords: models > > ### ** Examples > > ## > if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=1000} else {R=10} > > set.seed(66) > simnegbin = + function(X, beta, alpha) { + # Simulate from the Negative Binomial Regression + lambda = exp(X %*% beta) + y=NULL + for (j in 1:length(lambda)) + y = c(y,rnbinom(1,mu = lambda[j],size = alpha)) + return(y) + } > > nobs = 500 > nvar=2 # Number of X variables > alpha = 5 > Vbeta = diag(nvar)*0.01 > > # Construct the regdata (containing X) > 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) Starting Random Walk Metropolis Sampler for Negative Binomial Regression 500 obs; 2 covariates (including intercept); Prior Parameters: betabar [1] 0 0 A [,1] [,2] [1,] 0.01 0.00 [2,] 0.00 0.01 a [1] 0.5 b [1] 0.1 MCMC Parms: 10 reps; keeping every 1 th draw s_alpha = 2.93 s_beta = 2.071823 Initializing RW Increment Covariance Matrix... beta_mle = 0.7069472 0.1525201 alpha_mle = 6.272214 MCMC Iteration (est time to end - min) Total Time Elapsed: 0 > > cat(" betadraws ",fill=TRUE) betadraws > mat=apply(out$betadraw,2,quantile,probs=c(.01,.05,.5,.95,.99)) > mat=rbind(beta,mat); rownames(mat)[1]="beta"; print(mat) [,1] [,2] beta 0.6000000 0.2000000 1% 0.0000000 -0.2413808 5% 0.0000000 -0.2413322 50% 0.3808203 -0.1214455 95% 0.6746438 0.0000000 99% 0.6933401 0.0000000 > cat(" alphadraws ",fill=TRUE) alphadraws > 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) [,1] alpha 5.000000 1% 2.950495 5% 2.965392 50% 3.634174 95% 4.306654 99% 4.332895 > > > > cleanEx(); ..nameEx <- "rnmixGibbs" > > ### * rnmixGibbs > > flush(stderr()); flush(stdout()) > > ### Name: rnmixGibbs > ### Title: Gibbs Sampler for Normal Mixtures > ### Aliases: rnmixGibbs > ### Keywords: multivariate > > ### ** Examples > > ## > if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=2000} else {R=10} > > set.seed(66) > dim=5; k=3 # dimension of simulated data and number of "true" components > sigma = matrix(rep(0.5,dim^2),nrow=dim);diag(sigma)=1 > sigfac = c(1,1,1);mufac=c(1,2,3); compsmv=list() > for(i in 1:k) compsmv[[i]] = list(mu=mufac[i]*1:dim,sigma=sigfac[i]*sigma) > comps = list() # change to "rooti" scale > for(i in 1:k) comps[[i]] = list(mu=compsmv[[i]][[1]],rooti=solve(chol(compsmv[[i]][[2]]))) > pvec=(1:k)/sum(1:k) > > nobs=5000 > dm = rmixture(nobs,pvec,comps) > > Data=list(y=dm$x) > ncomp=9 > Prior=list(ncomp=ncomp,a=c(rep(1,ncomp))) > Mcmc=list(R=R,keep=1) > out=rnmixGibbs(Data=Data,Prior=Prior,Mcmc=Mcmc) Starting Gibbs Sampler for Mixture of Normals 5000 observations on 5 dimensional data using 9 mixture components Prior Parms: mu_j ~ N(mubar,Sigma (x) A^-1) mubar = [,1] [,2] [,3] [,4] [,5] [1,] 0 0 0 0 0 precision parm for prior variance of mu vectors (A)= 0.01 Sigma_j ~ IW(nu,V) nu= 7 V = [,1] [,2] [,3] [,4] [,5] [1,] 7 0 0 0 0 [2,] 0 7 0 0 0 [3,] 0 0 7 0 0 [4,] 0 0 0 7 0 [5,] 0 0 0 0 7 Dirichlet parameters [1] 1 1 1 1 1 1 1 1 1 Mcmc Parms: R= 10 keep= 1 starting value for z z 1 2 3 4 5 6 7 8 9 555 555 555 555 555 555 555 555 560 MCMC Iteration (est time to end -min) > > 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) mu and posterior expectation of mu > print(mat) [,1] [,2] [,3] [,4] [,5] true 2.333333 4.666667 7.000000 9.333333 11.66667 post expect 2.331035 4.662917 7.030559 9.371743 11.73343 > mat=rbind(tmom$sd,pmom$sd) > rownames(mat)=c("true","post expect") > cat(" std dev and posterior expectation of sd",fill=TRUE) std dev and posterior expectation of sd > print(mat) [,1] [,2] [,3] [,4] [,5] true 1.000000 1.000000 1.000000 1.000000 1.000000 post expect 1.221569 1.801891 2.485522 3.235971 3.853351 > 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) corr and posterior expectation of corr > print(t(mat)) true post expect [1,] 1.0 1.0000000 [2,] 0.5 0.7373174 [3,] 0.5 0.7220982 [4,] 0.5 0.6832841 [5,] 0.5 0.6858473 [6,] 0.5 0.7373174 [7,] 1.0 1.0000000 [8,] 0.5 0.8809139 [9,] 0.5 0.8789592 [10,] 0.5 0.8831872 [11,] 0.5 0.7220982 [12,] 0.5 0.8809139 [13,] 1.0 1.0000000 [14,] 0.5 0.9441418 [15,] 0.5 0.9413499 [16,] 0.5 0.6832841 [17,] 0.5 0.8789592 [18,] 0.5 0.9441418 [19,] 1.0 1.0000000 [20,] 0.5 0.9624447 [21,] 0.5 0.6858473 [22,] 0.5 0.8831872 [23,] 0.5 0.9413499 [24,] 0.5 0.9624447 [25,] 1.0 1.0000000 > > > > > cleanEx(); ..nameEx <- "rscaleUsage" > > ### * rscaleUsage > > flush(stderr()); flush(stdout()) > > ### Name: rscaleUsage > ### Title: MCMC Algorithm for Multivariate Ordinal Data with Scale Usage > ### Heterogeneity. > ### Aliases: rscaleUsage > ### Keywords: models > > ### ** Examples > > ## > if(nchar(Sys.getenv("LONG_TEST")) != 0) + { + data(customerSat) + surveydat = list(k=10,x=as.matrix(customerSat)) + + R=1000 + mcmc = list(R=R) + set.seed(66) + out=rscaleUsage(Data=surveydat,Mcmc=mcmc) + + cat(" mudraws ",fill=TRUE) + mat=apply(out$mudraw,2,quantile,probs=c(.01,.05,.5,.95,.99)) + print(mat) + + } > > > > cleanEx(); ..nameEx <- "rtrun" > > ### * rtrun > > flush(stderr()); flush(stdout()) > > ### Name: rtrun > ### Title: Draw from Truncated Univariate Normal > ### Aliases: rtrun > ### Keywords: distribution > > ### ** Examples > > ## > set.seed(66) > rtrun(mu=c(rep(0,10)),sigma=c(rep(1,10)),a=c(rep(0,10)),b=c(rep(2,10))) [1] 1.9180411 1.2589217 0.7708917 0.5216277 0.9001720 0.4705582 0.5315487 [8] 1.2935472 0.4667398 1.5143196 > > > > cleanEx(); ..nameEx <- "runireg" > > ### * runireg > > flush(stderr()); flush(stdout()) > > ### Name: runireg > ### Title: Draw from Posterior for Univariate Regression > ### Aliases: runireg > ### Keywords: regression > > ### ** Examples > > if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=1000} else {R=10} > set.seed(66) > 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 > > betadraw=matrix(double(R*2),ncol=2) > sigsqdraw=double(R) > for (rep in 1:R) + {out=runireg(y,X,betabar,A,nu,ssq);betadraw[rep,]=out$beta + sigsqdraw[rep]=out$sigmasq} > > cat(" Betadraws ",fill=TRUE) Betadraws > mat=apply(betadraw,2,quantile,probs=c(.01,.05,.5,.95,.99)) > mat=rbind(beta,mat); rownames(mat)[1]="beta"; print(mat) [,1] [,2] beta 1.0000000 2.000000 1% 0.7471191 2.193845 5% 0.7545305 2.205619 50% 0.8313063 2.393986 95% 0.9356706 2.615551 99% 0.9466812 2.659226 > cat(" Sigma-sq draws",fill=TRUE) Sigma-sq draws > cat(" sigma-sq= ",sigsq,fill=TRUE) sigma-sq= 0.25 > print(quantile(sigsqdraw,probs=c(.01,.05,.5,.95,.99))) 1% 5% 50% 95% 99% 0.2061553 0.2182414 0.2542022 0.2785026 0.2904749 > > > > cleanEx(); ..nameEx <- "runiregGibbs" > > ### * runiregGibbs > > flush(stderr()); flush(stdout()) > > ### Name: runiregGibbs > ### Title: Gibbs Sampler for Univariate Regression > ### Aliases: runiregGibbs > ### Keywords: regression > > ### ** Examples > > if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=2000} else {R=10} > set.seed(66) > 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 > > 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) Starting Gibbs Sampler for Univariate Regression Model with 100 observations Prior Parms: betabar [1] 0 0 A [,1] [,2] [1,] 0.05 0.00 [2,] 0.00 0.05 nu = 3 ssq= 1 MCMC parms: R= 10 keep= 1 MCMC Iteration (est time to end - min) Total Time Elapsed: 0 > cat(" Betadraws ",fill=TRUE) Betadraws > mat=apply(out$betadraw,2,quantile,probs=c(.01,.05,.5,.95,.99)) > mat=rbind(beta,mat); rownames(mat)[1]="beta"; print(mat) [,1] [,2] beta 1.0000000 2.000000 1% 0.8021015 2.241618 5% 0.8066362 2.249716 50% 0.8946017 2.345063 95% 0.9550541 2.519727 99% 0.9739075 2.527192 > cat(" Sigma-sq draws",fill=TRUE) Sigma-sq draws > cat(" sigma-sq= ",sigsq,fill=TRUE) sigma-sq= 0.25 > print(quantile(out$sigmasqdraw,probs=c(.01,.05,.5,.95,.99))) 1% 5% 50% 95% 99% 0.2203772 0.2217524 0.2452490 0.2768560 0.2865918 > > > > cleanEx(); ..nameEx <- "rwishart" > > ### * rwishart > > flush(stderr()); flush(stdout()) > > ### Name: rwishart > ### Title: Draw from Wishart and Inverted Wishart Distribution > ### Aliases: rwishart > ### Keywords: multivariate > > ### ** Examples > > ## > set.seed(66) > rwishart(5,diag(3))$IW [,1] [,2] [,3] [1,] 0.07984024 0.0270185 0.03571638 [2,] 0.02701850 0.3216721 0.11751868 [3,] 0.03571638 0.1175187 0.34201739 > > > > ### *