[r-cran-mcmcpack] 46/90: Imported Upstream version 1.0-6
Andreas Tille
tille at debian.org
Fri Dec 16 09:07:39 UTC 2016
This is an automated email from the git hooks/post-receive script.
tille pushed a commit to branch master
in repository r-cran-mcmcpack.
commit b1c6797e1e78b6f989cd070f49615fa692439d3c
Author: Andreas Tille <tille at debian.org>
Date: Fri Dec 16 08:07:17 2016 +0100
Imported Upstream version 1.0-6
---
DESCRIPTION | 8 +-
HISTORY | 8 +
R/MCMCoprobit.R | 78 +++++-----
R/MCMCquantreg.R | 37 ++---
R/SSVSquantreg.R | 90 +++++++++++
R/SSVSquantregsummary.R | 82 ++++++++++
R/distn.R | 2 +-
man/MCMCirtKdHet.Rd | 4 +-
man/MCMCquantreg.Rd | 74 ++++-----
man/QRSSVSplot.Rd | 50 ++++++
man/QRSSVSsummary.Rd | 42 ++++++
man/SSVSquantreg.Rd | 136 +++++++++++++++++
man/mptable.Rd | 33 ++++
man/topmodels.Rd | 42 ++++++
src/MCMCdynamicIRT1d-b.cc | 10 +-
src/MCMCdynamicIRT1d.cc | 10 +-
src/MCMCfcds.h | 377 ++++++++++++++++++++++++++++++++++------------
src/MCMCirtHier1d.cc | 13 +-
src/MCMCmedreg.cc | 138 -----------------
src/MCMCquantreg.cc | 27 ++--
src/SSVSquantreg.cc | 304 +++++++++++++++++++++++++++++++++++++
21 files changed, 1185 insertions(+), 380 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index cbe1645..fc000d0 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
Package: MCMCpack
-Version: 1.0-5
-Date: 2009-11-28
+Version: 1.0-6
+Date: 2009-5-11
Title: Markov chain Monte Carlo (MCMC) Package
Author: Andrew D. Martin <admartin at wustl.edu>, Kevin M. Quinn
<kevin_quinn at harvard.edu>, Jong Hee Park <jhp at uchicago.edu>
@@ -19,6 +19,6 @@ Description: This package contains functions to perform Bayesian
License: GPL-2
SystemRequirements: gcc (>= 4.0)
URL: http://mcmcpack.wustl.edu
-Packaged: 2009-11-28 16:57:19 UTC; adm
+Packaged: 2010-05-11 16:18:44 UTC; adm
Repository: CRAN
-Date/Publication: 2009-11-29 17:22:58
+Date/Publication: 2010-05-11 18:43:21
diff --git a/HISTORY b/HISTORY
index b875b3b..5150247 100644
--- a/HISTORY
+++ b/HISTORY
@@ -2,6 +2,14 @@
// Changes and Bug Fixes
//
+
+1.0-5 to 1.0-6
+ * added SSVSquantreg [written by Craig Reed]
+ * added functions to handle SSVSquantreg output [written by Craig Reed]
+ * revisions to the MCMCquantreg function [written by Craig Reed]
+ * revisions to the MCMCquantreg function [written by Craig Reed]
+ * fixed parameterization issue with rinvgamma [thanks to Julian Stander]
+
1.0-4 to 1.0-5
* added heteroskedastic IRT model [written by Ben Lauderdale]
* minor changes to error status in hierarchical IRT [written by Mike Malecki]
diff --git a/R/MCMCoprobit.R b/R/MCMCoprobit.R
index ec58c9e..ef7551b 100644
--- a/R/MCMCoprobit.R
+++ b/R/MCMCoprobit.R
@@ -32,9 +32,9 @@
lecuyer.stream <- seeds[[3]]
## extract X, Y, and variable names from the model formula and frame
- call <- match.call()
- mt <- terms(formula, data=data)
- if(missing(data)) data <- sys.frame(sys.parent())
+ call <- match.call()
+ mt <- terms(formula, data=data)
+ if(missing(data)) data <- sys.frame(sys.parent())
mf <- match.call(expand.dots = FALSE)
mf$burnin <- mf$mcmc <- mf$b0 <- mf$B0 <- NULL
mf$thin <- mf$... <- mf$tune <- mf$tdf <- mf$verbose <- mf$seed <- NULL
@@ -111,68 +111,68 @@
stop("Please respecify and call MCMCoprobit() again.\n")
}
- ## prior for alpha error checking
- if(is.null(dim(a0))) {
- a0 <- a0 * matrix(1, ncat-1, 1)
- }
- if((dim(a0)[1] != ncat-1) || (dim(a0)[2] != 1)) {
- cat("N(a0,A0) prior a0 not conformable.\n")
- stop("Please respecify and call MCMCoprobit() again.\n")
- }
- if(is.null(dim(A0))) {
- A0 <- A0 + diag(ncat - 1)
- }
- if((dim(A0)[1] != ncat - 1) || (dim(A0)[2] != ncat - 1)) {
- cat("N(a0, A0) prior A0 not conformable.\n")
- stop("Please respecify and call MCMCoprobit() again.\n")
- }
-
+ ## prior for alpha error checking
+ if(is.null(dim(a0))) {
+ a0 <- a0 * matrix(1, ncat-1, 1)
+ }
+ if((dim(a0)[1] != ncat-1) || (dim(a0)[2] != 1)) {
+ cat("N(a0,A0) prior a0 not conformable.\n")
+ stop("Please respecify and call MCMCoprobit() again.\n")
+ }
+ if(is.null(dim(A0))) {
+ A0 <- A0 + diag(ncat - 1)
+ }
+ if((dim(A0)[1] != ncat - 1) || (dim(A0)[2] != ncat - 1)) {
+ cat("N(a0, A0) prior A0 not conformable.\n")
+ stop("Please respecify and call MCMCoprobit() again.\n")
+ }
+
## form gamma starting values (note: not changeable)
gamma <- matrix(NA,ncat+1,1)
gamma[1] <- -300
gamma[2] <- 0
gamma[3:ncat] <- (polr.out$zeta[2:(ncat-1)] - polr.out$zeta[1])*.588
gamma[ncat+1] <- 300
-
+
## posterior sample
sample <- matrix(data=0, mcmc/thin, K + ncat + 1)
-
+
## call C++ code to draw sample
- nY <- as.matrix(as.numeric(Y))
-
- ## mcmc.method
- cowles <- as.integer(1)
- if(mcmc.method[1]=="AC") {cowles <- as.integer(0)}
-
- ## form the tuning parameter
- tune <- vector.tune(tune, ncat-1)
+ nY <- as.matrix(as.numeric(Y))
+
+ ## mcmc.method
+ cowles <- as.integer(1)
+ if(mcmc.method[1]=="AC") {cowles <- as.integer(0)}
+
+ ## form the tuning parameter
+ tune <- vector.tune(tune, ncat-1)
posterior <- NULL
-
+
auto.Scythe.call(output.object="posterior", cc.fun.name="MCMCoprobit",
sample.nonconst=sample, Y=as.integer(Y), nY=nY, X=X,
burnin=as.integer(burnin),
mcmc=as.integer(mcmc), thin=as.integer(thin),
tune=tune, tdf=as.double(tdf),
- lecuyer=as.integer(lecuyer),
+ lecuyer=as.integer(lecuyer),
seedarray=as.integer(seed.array),
lecuyerstream=as.integer(lecuyer.stream),
verbose=as.integer(verbose), beta=beta.start,
gamma=gamma, b0=b0, B0=B0, a0=a0, A0=A0,
- cowles=as.integer(cowles))
+ cowles=as.integer(cowles))
## put together matrix and build MCMC object to return
sample <- matrix(posterior$sampledata, posterior$samplerow,
posterior$samplecol, byrow=FALSE)
- if(mcmc.method[1]=="AC"){
- sample[ , 1] <- sample[, 1] - sample[, K+2] ## post-MCMC normalization
- sample[ , (K+2):(K+ncat)] <- sample[ , (K+2):(K+ncat)] - sample[, K+2] ## post-MCMC normalization
- }
- sample <- sample[, c(1:K, (K+3):(K+ncat))]
-
+ if(mcmc.method[1]=="AC"){
+ sample[ , 1] <- sample[, 1] - sample[, K+2] ## post-MCMC normalization
+ sample[ , (K+2):(K+ncat)] <- sample[ , (K+2):(K+ncat)] - sample[, K+2] ## post-MCMC normalization
+ }
+ sample <- sample[, c(1:K, (K+3):(K+ncat))]
+
output <- mcmc(data=sample, start=burnin+1, end=burnin+mcmc, thin=thin)
xnames <- c(X.names, paste("gamma", 2:(ncat-1), sep=""))
varnames(output) <- xnames
attr(output, "title") <- "MCMCoprobit Posterior Sample"
-
+
return(output)
}
diff --git a/R/MCMCquantreg.R b/R/MCMCquantreg.R
index c6f1d66..8bf66f9 100644
--- a/R/MCMCquantreg.R
+++ b/R/MCMCquantreg.R
@@ -1,7 +1,7 @@
"MCMCquantreg" <-
- function(formula, tau=0.5, data=NULL, burnin = 1000, mcmc = 10000,
- thin=1, verbose = 0, seed = NA, beta.start = NA,
- b0 = 0, B0 = 0, c0 = 0.001, d0 = 0.001,
+ function(formula, data=NULL, tau=0.5, burnin = 1000, mcmc = 10000,
+ thin=1, verbose = 0, seed = sample(1:1000000,1), beta.start = NA,
+ b0 = 0, B0 = 0,
...) {
## checks
@@ -10,7 +10,7 @@
cl <- match.call()
if (tau<=0 || tau>=1){
- stop("tau must be in (0,1).\n Please respecify and call again.\n")
+ stop("tau must be in (0,1). \nPlease respecify and call again.\n")
}
## seeds
@@ -27,15 +27,14 @@
K <- ncol(X) # number of covariates
## starting values and priors
- ols.fit <- lm(formula)
+ ols.fit <- lm(formula, data=data)
defaults <- matrix(coef(ols.fit),K,1)
defaults[1] <- defaults[1]+summary(ols.fit)$sigma*qnorm(tau)
beta.start <- coef.start(beta.start, K, formula, family=gaussian, data, defaults=defaults)
mvn.prior <- form.mvn.prior(b0, B0, K)
b0 <- mvn.prior[[1]]
B0 <- mvn.prior[[2]]
- check.ig.prior(c0, d0)
-
+
B0.eigenvalues <- eigen(B0)$values
if (min(B0.eigenvalues) < 0){
@@ -45,24 +44,10 @@
## define holder for posterior sample
- sample <- matrix(data=0, mcmc/thin, K+1)
+ sample <- matrix(data=0, mcmc/thin, K)
posterior <- NULL
- if (tau==0.5) {
- ## call C++/Scythe function "MCMCmedreg" to draw samples
- auto.Scythe.call(output.object="posterior", cc.fun.name="MCMCmedreg",
- sample.nonconst=sample, Y=Y, X=X, burnin=as.integer(burnin),
- mcmc=as.integer(mcmc), thin=as.integer(thin),
- lecuyer=as.integer(lecuyer),
- seedarray=as.integer(seed.array),
- lecuyerstream=as.integer(lecuyer.stream),
- verbose=as.integer(verbose), betastart=beta.start,
- b0=b0, B0=B0, c0=as.double(c0), d0=as.double(d0), package="MCMCpack")
-}
-
- else{
-
- ## call C++/Scythe function "MCMCquantreg" to draw samples
+ ## call C++ code to draw samples
auto.Scythe.call(output.object="posterior", cc.fun.name="MCMCquantreg",
sample.nonconst=sample, tau=as.double(tau), Y=Y, X=X,
burnin=as.integer(burnin),
@@ -71,12 +56,12 @@
seedarray=as.integer(seed.array),
lecuyerstream=as.integer(lecuyer.stream),
verbose=as.integer(verbose), betastart=beta.start,
- b0=b0, B0=B0, c0=as.double(c0), d0=as.double(d0), package="MCMCpack")
-}
+ b0=b0, B0=B0, package="MCMCpack")
+
## pull together matrix and build MCMC object to return
output <- form.mcmc.object(posterior,
- names=c(xnames,"sigma"),
+ names=xnames,
title="MCMCquantreg Posterior Sample",
y=Y, call=cl)
diff --git a/R/SSVSquantreg.R b/R/SSVSquantreg.R
new file mode 100644
index 0000000..1a1e470
--- /dev/null
+++ b/R/SSVSquantreg.R
@@ -0,0 +1,90 @@
+"SSVSquantreg" <-
+ function(formula, data=NULL, tau=0.5, include=NULL, burnin=1000, mcmc = 10000, thin=1,
+ verbose = 0, seed = sample(1:1000000,1), pi0a0=1, pi0b0=1,
+ ...) {
+
+ ## checks
+ check.offset(list(...))
+ if (pi0a0<=0 || pi0b0<=0){
+ stop("Parameters pi0a0 and pi0b0 must be positive. \nPlease respecify and call again.\n")
+ }
+ cl <- match.call()
+ if (tau<=0 || tau>=1){
+ stop("tau must be in (0,1). \nPlease respecify and call again.\n")
+ }
+
+
+ ## seeds
+ seeds <- form.seeds(seed)
+ lecuyer <- seeds[[1]]
+ seed.array <- seeds[[2]]
+ lecuyer.stream <- seeds[[3]]
+
+ ## form response and model matrices
+ holder <- parse.formula(formula, data=data)
+ Y <- holder[[1]]
+ X <- holder[[2]]
+ xnames <- holder[[3]]
+ K <- ncol(X) # number of covariates
+ q <- length(include) #number of covariates that are pre-specified to appear in the model
+
+ if (!is.null(include)){
+ if(is.character(include)){
+ if(!all(include%in%xnames)){
+ include.positions<-NA
+ }
+ else{
+ include.positions<-match(include, xnames)
+ }
+ }
+ else{
+ if(max(include)>length(xnames)){
+ include.positions<-NA
+ }
+ else{
+ include.positions<-include
+ }
+ }
+ if (any(is.na(include.positions))){
+ stop("One or more covariates to be included are not present in the design matrix\n or one or more indices are out of range. Please respecify and call again.\n")
+ }
+
+ ## Bring those covariates that are pre-specified to appear in the model to the first positions in the X matrix
+ X <- cbind(X[,include.positions], X[,-include.positions])
+ xnames <- c(xnames[include.positions], xnames[-include.positions])
+ }
+
+ ## define holder for posterior sample
+ sample <- matrix(data=0, mcmc/thin, 2*K)
+ posterior <- NULL
+
+ ## call C++ code to draw sample
+ auto.Scythe.call(output.object="posterior", cc.fun.name="SSVSquantreg",
+ sample.nonconst=sample, tau=as.double(tau), Y=Y, X=X, q=as.integer(q),
+ burnin=as.integer(burnin), mcmc=as.integer(mcmc), thin=as.integer(thin),
+ lecuyer=as.integer(lecuyer),
+ seedarray=as.integer(seed.array),
+ lecuyerstream=as.integer(lecuyer.stream),
+ verbose=as.integer(verbose),
+ pi0a0 = as.double(pi0a0), pi0b0=as.double(pi0b0),
+ package="MCMCpack")
+
+
+output <- form.mcmc.object(posterior, names=rep(xnames, times=2),
+ title="SSVSquantreg Posterior Sample",
+ y=Y, call=cl
+ )
+
+gammasample<-output[,1:K]
+attr(gammasample, "tau")<-tau
+attr(gammasample, "xnames")<-xnames
+attr(gammasample, "class")<-"qrssvs"
+betasample<-output[,-(1:K)]
+attr(betasample,"tau")<-tau
+attr(gammasample, "call") <- attr(betasample, "call") <- cl
+return(list(gamma=gammasample,beta=betasample))
+}
+
+
+
+
diff --git a/R/SSVSquantregsummary.R b/R/SSVSquantregsummary.R
new file mode 100644
index 0000000..b9d227a
--- /dev/null
+++ b/R/SSVSquantregsummary.R
@@ -0,0 +1,82 @@
+"print.qrssvs"<-function(x, ...){
+x.orig<-x
+cat("Quantile regression stochastic search \nvariable selection (QR-SSVS) output:\nStart = ",
+ attr(x,"mcpar")[1], "\nEnd = ",
+ attr(x,"mcpar")[2], "\nThinning interval = ",
+ attr(x,"mcpar")[3], "\n")
+
+attr(x, "mcpar") <- NULL
+attr(x, "class") <- NULL
+NextMethod("print", ...)
+invisible(x.orig)
+}
+
+"mptable"<-function(qrssvs){
+if (!is(qrssvs, "qrssvs")){
+stop("Can only be used on objects of class qrssvs.\n")
+}
+ssvs.start <- attr(qrssvs, "mcpar")[1]
+ssvs.end <- attr(qrssvs, "mcpar")[2]
+ssvs.thin <- attr(qrssvs, "mcpar")[3]
+nstore <- (ssvs.end-ssvs.start)/ssvs.thin + 1
+probs<-apply(qrssvs,2,function(z){length(which(z==1))})/nstore
+return(data.frame(Probability=probs))
+}
+
+"topmodels"<-function(qrssvs, nmodels=5, abbreviate=FALSE, minlength=3){
+if (!is(qrssvs, "qrssvs")){
+stop("Can only be used on objects of class qrssvs.\n")
+}
+ssvs.start <- attr(qrssvs, "mcpar")[1]
+ssvs.end <- attr(qrssvs, "mcpar")[2]
+ssvs.thin <- attr(qrssvs, "mcpar")[3]
+nstore <- (ssvs.end-ssvs.start)/ssvs.thin + 1
+xnames <- attr(qrssvs, "xnames")
+if (abbreviate){
+xnames <- abbreviate(xnames, minlength)
+}
+model.list<-apply(qrssvs,1,function(z)xnames[which(z==1)])
+model.vector<-sapply(model.list, function(z)paste(z, collapse=","))
+model.count<-sort(table(model.vector), decreasing=T)/nstore
+if (nmodels>length(model.count)){
+warning("Number of models requested exceeds total number of models visited.\n")
+}
+if (rownames(model.count)[1]==""){
+rownames(model.count)[1]<-"Null model"
+}
+return(data.frame(Probability=model.count[1:(min(nmodels, length(model.count)))]))
+}
+
+"plot.qrssvs"<-function(x, ...){
+probs<-mptable(x)
+dotplot(as.matrix(probs),
+ panel=function(x, y, ...){
+ panel.abline(v=0.5, lty=3)
+ panel.dotplot(x, y, ...)
+ },
+ origin=0, type=c("p","h"), pch=16,
+ xlim=c(-0.05,1.05),
+ scales = list(x = list(at = c(0,0.2,0.4,0.5,0.6,0.8,1))),
+ xlab="Marginal inclusion probability", ...)
+}
+
+"summary.qrssvs"<-function(object, ...){
+covnames <- attr(object, "xnames")
+probs<-mptable(object)
+median.model<-covnames[probs>=0.5]
+results<-probs
+attr(results, "median.model")<-median.model
+attr(results, "tau")<-attr(object, "tau")
+class(results)<-"summary.qrssvs"
+return(results)
+}
+
+print.summary.qrssvs<-function(x, digits=max(3, .Options$digits-3), ...){
+attr(x, "class")<-"data.frame"
+cat("\nMarginal inclusion probability of each predictor:\n\n")
+print(x, digits=digits, ...)
+cat("\nFor tau = ", attr(x,"tau"), ", the median probability model \nincludes the following predictors:\n\n",
+paste(attr(x, "median.model"), collapse=", "), ".\n\n", sep="")
+invisible(x)
+}
+
diff --git a/R/distn.R b/R/distn.R
index 00ff600..49e5ba8 100644
--- a/R/distn.R
+++ b/R/distn.R
@@ -221,7 +221,7 @@
"rinvgamma" <-
function(n, shape, scale = 1) {
- return(1 / rgamma(n, shape, scale))
+ return(1 / rgamma(n, shape, scale=scale))
}
##
diff --git a/man/MCMCirtKdHet.Rd b/man/MCMCirtKdHet.Rd
index 8acf1d6..1ec62a0 100644
--- a/man/MCMCirtKdHet.Rd
+++ b/man/MCMCirtKdHet.Rd
@@ -156,7 +156,7 @@ Modified from \code{\link[MCMCpack]{MCMCirtKd}} and \code{\link[MCMCpack]{MCMCor
}
\examples{
-
+\dontrun{
data(Senate)
Y <- as.matrix(Senate[,6:677])
@@ -184,7 +184,7 @@ legend(x="topleft",legend=c("Point sizes proportional to estimated legislator",
"variance under heteroskedastic model.","Some legislators with large variance have",
"more extreme estimated ideal points under the","heteroskedastic model because their",
"deviations from the party line are attributable","to idiosyncrasy rather than moderation."),cex=0.5)
-
+}
}
\keyword{models}
diff --git a/man/MCMCquantreg.Rd b/man/MCMCquantreg.Rd
index e9bdddb..d718315 100644
--- a/man/MCMCquantreg.Rd
+++ b/man/MCMCquantreg.Rd
@@ -7,27 +7,25 @@
\description{
This function fits quantile regression models under Bayesian inference.
The function samples from the posterior distribution using Gibbs sampling with data augmentation.
- A multivariate normal prior is assumed for \eqn{\beta}{beta} and an inverse Gamma prior
- is assumed for \eqn{\sigma}{sigma}. The user supplies the prior parameters.
+ A multivariate normal prior is assumed for \eqn{\beta}{beta}. The user supplies the prior parameters.
A sample of the posterior distribution is returned as an mcmc object,
which can then be analysed by functions in the coda package.
}
\usage{
-MCMCquantreg(formula, tau=0.5, data = NULL, burnin = 1000,
- mcmc = 10000, thin = 1, verbose = 0, seed = NA,
- beta.start = NA, b0 = 0, B0 = 0, c0 = 0.001, d0 = 0.001, ...)
+MCMCquantreg(formula, data = NULL, tau=0.5, burnin = 1000,
+ mcmc = 10000, thin = 1, verbose = 0, seed = sample(1:1000000,1),
+ beta.start = NA, b0 = 0, B0 = 0, ...)
}
\arguments{
\item{formula}{ Model formula. }
- \item{tau}{The quantile of interest. Must be between 0 and 1. The default value of 0.5 corresponds to median regression.}
-
-
\item{data}{ Data frame. }
+ \item{tau}{The quantile of interest. Must be between 0 and 1. The default value of 0.5 corresponds to median regression.}
+
\item{burnin}{ The number of burn-in iterations for the sampler. }
\item{mcmc}{ The number of MCMC iterations after burnin. }
@@ -41,9 +39,9 @@ MCMCquantreg(formula, tau=0.5, data = NULL, burnin = 1000,
\eqn{\beta}{beta} and \eqn{\sigma}{sigma} are printed to
the screen every \code{verbose}th iteration. }
- \item{seed}{ The seed for the random number generator. If NA, the Mersenne
+ \item{seed}{ The seed for the random number generator. If NA, the Mersenne
Twister generator is used with default seed 12345; if an integer is
- passed it is used to seed the Mersenne twister. The user can also
+ passed it is used to seed the Mersenne twister. The default value for this argument is a random integer between 1 and 1,000,000. This default value ensures that if the function is used again with a different value of \eqn{\tau}{tau}, it is extremely unlikely that the seed will be identical. The user can also
pass a list of length two to use the L'Ecuyer random number generator,
which is suitable for parallel computation. The first element of the
list is the L'Ecuyer seed, which is a vector of length six or NA (if NA
@@ -78,51 +76,37 @@ MCMCquantreg(formula, tau=0.5, data = NULL, burnin = 1000,
as the prior precision of \eqn{\beta}{beta}. Default value of 0 is equivalent to
an improper uniform prior for \eqn{\beta}{beta}. }
- \item{c0}{ \eqn{c_0/2}{c0/2} is the shape parameter for the inverse
- Gamma prior on \eqn{\sigma}{sigma}. }
-
- \item{d0}{ \eqn{d_0/2}{d0/2} is the scale parameter for the
- inverse Gamma prior on \eqn{\sigma}{sigma}.}
+ \item{\dots}{ further arguments to be passed }
+}
- \item{\dots}{ further arguments to be passed }
+\value{
+ An mcmc object that contains the posterior sample. This
+ object can be summarised by functions provided by the coda package.
}
\details{
\code{MCMCquantreg} simulates from the posterior distribution using
- a partially collapsed Gibbs sampler with data augmentation (see \url{http://people.brunel.ac.uk/~mastkky/}).
- This involves a multivariate Normal draw for \eqn{\beta}{beta}, and an
- inverse Gamma draw for \eqn{\sigma}{sigma}.
- The augmented data are drawn conditionally from the inverse Gaussian distribution. The simulation
+ Gibbs sampling with data augmentation (see \url{http://people.brunel.ac.uk/~mastkky/}).
+ \eqn{\beta}{beta} are drawn from a multivariate normal distribution. The augmented data are drawn conditionally from the inverse Gaussian distribution. The simulation
is carried out in compiled C++ code to maximise efficiency. Please consult
the coda documentation for a comprehensive list of functions that can be
used to analyse the posterior sample.
- The model takes the following form:
- \deqn{y_i = x_i ' \beta + \varepsilon_i}{y_i = x_i'beta + epsilon_i}
- The errors are assumed to have an Asymmetric Laplace distribution with \eqn{\tau}{tau}th quantile equal to zero
- and shape parameter equal to \eqn{\sigma}{sigma}:
- \deqn{\varepsilon_i(\tau) \sim \mathcal{AL}(0, \sigma, \tau)}{epsilon_i(tau) ~ AL(0,
- sigma, tau),}
- where \eqn{\beta}{beta} and \eqn{\sigma}{sigma} depend on \eqn{\tau}{tau}.
- We assume standard, semi-conjugate priors:
- \deqn{\beta \sim \mathcal{N}(b_0,B_0^{-1})}{beta ~ N(b0,B0^(-1)),}
- and
- \deqn{\sigma^{-1} \sim \mathcal{G}amma(c_0/2, d_0/2)}{sigma^(-1) ~
- Gamma(c0/2, d0/2),}
- where \eqn{\beta}{beta} and \eqn{\sigma^{-1}}{sigma^(-1)} are assumed
- \emph{a priori} independent of each other and the parameters associated with the other quantiles. Only starting values for
+ We assume the model
+\deqn{Q_{\tau}(y_i|x_i) = x_i'\beta}{Q_tau(y_i|x_i) = x_i'beta}, where \eqn{Q_{\tau}(y_i|x_i)}{Q_tau(y_i|x_i)} denotes the conditional \eqn{\tau}{tau}th quantile of \eqn{y_i}{y_i} given \eqn{x_i}{x_i}, and \eqn{\beta=\beta(\tau)}{beta=beta(tau)} are the regression parameters possibly dependent on \eqn{\tau}{tau}.
+The likelihood is formed based on assuming independent Asymmetric Laplace distributions on the \eqn{y_i}{y_i} with skewness parameter \eqn{\tau}{tau} and location parameters \eqn{x_i'\beta}{x_i'beta}. This assumption ensures that the likelihood function is maximised by the \eqn{\tau}{tau}th conditional quantile of the response variable.
+ We assume standard, semi-conjugate priors on \eqn{\beta}{beta}:
+ \deqn{\beta \sim \mathcal{N}(b_0,B_0^{-1})}{beta ~ N(b0,B0^(-1))}.
+ Only starting values for
\eqn{\beta}{beta} are allowed for this sampler.
}
-\value{
- An mcmc object that contains the posterior sample. This
- object can be summarised by functions provided by the coda package.
-}
+\author{ Craig Reed}
\references{ Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007.
\emph{Scythe Statistical Library 1.2.} \url{http://scythe.wustl.edu}.
- Craig Reed and Keming Yu. 2009. "A Partially Collapsed Gibbs Sampler for Bayesian Quantile Regression." Technical Report.
+ Craig Reed and Keming Yu. 2009. "An Efficient Gibbs Sampler for Bayesian Quantile Regression." Technical Report.
Keming Yu and Jin Zhang. 2005. "A Three Parameter Asymmetric Laplace Distribution and it's extensions."
\emph{Communications in Statistics - Theory and Methods}, 34, 1867-1879.
@@ -130,11 +114,6 @@ MCMCquantreg(formula, tau=0.5, data = NULL, burnin = 1000,
Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2002.
\emph{Output Analysis and Diagnostics for MCMC (CODA)}.
\url{http://www-fis.iarc.fr/coda/}.}
-
-\author{ Craig Reed}
-
-\seealso{ \code{\link[MCMCpack]{MCMCregress}}, \code{\link[coda]{plot.mcmc}},
- \code{\link[coda]{summary.mcmc}}, \code{\link[stats]{lm}}, \code{\link[quantreg]{rq}} }
\examples{
\dontrun{
@@ -154,3 +133,10 @@ summary(posterior_95)
}
\keyword{models}
+
+\seealso{
+\code{\link[MCMCpack]{MCMCregress}},
+\code{\link[coda]{plot.mcmc}},
+\code{\link[coda]{summary.mcmc}},
+\code{\link[stats]{lm}},
+\code{\link[quantreg]{rq}}}
diff --git a/man/QRSSVSplot.Rd b/man/QRSSVSplot.Rd
new file mode 100644
index 0000000..991bd16
--- /dev/null
+++ b/man/QRSSVSplot.Rd
@@ -0,0 +1,50 @@
+\name{plot.qrssvs}
+\alias{plot.qrssvs}
+\title{Plot output from quantile regression stochastic search variable selection (QR-SSVS).}
+
+\description{This function produces a Trellis plot of the predictors on the y-axis versus the marginal posterior probability of inclusion on the x-axis.}
+
+\usage{
+\method{plot}{qrssvs}(x, \dots)
+}
+
+\arguments{
+ \item{x}{An object of class \code{qrssvs}. Typically this will be the \code{gamma} component of the list returned by \code{SSVSquantreg}.}
+
+ \item{\dots}{Further arguments}
+}
+
+\value{An object with class \code{"trellis"}. The associated \code{\link[lattice:update.trellis]{update}} and \code{\link[lattice:print.trellis]{print}} methods are documented in the "Lattice" package.}
+
+\author{Craig Reed}
+
+\references{
+
+ Deepayan Sarkar. 2008. \emph{lattice: Lattice Graphics.} R package version 0.17-17
+}
+
+\examples{
+\dontrun{
+set.seed(1)
+epsilon<-rnorm(100)
+set.seed(2)
+x<-matrix(rnorm(1000),100,10)
+y<-x[,1]+x[,10]+epsilon
+qrssvs<-SSVSquantreg(y~x)
+plot(qrssvs$gamma)
+## Modify the graph by increasing the fontsize on the axes
+qrssvsplot<-plot(qrssvs$gamma)
+update(qrssvsplot, scales=list(cex=3))
+}
+}
+
+\keyword{models}
+
+\seealso{
+\code{\link[MCMCpack]{SSVSquantreg}},
+\code{\link[MCMCpack]{mptable}},
+\code{\link[lattice:Lattice]{Lattice}} for a brief introduction to
+ lattice displays and links to further documentation.
+}
+
+
diff --git a/man/QRSSVSsummary.Rd b/man/QRSSVSsummary.Rd
new file mode 100644
index 0000000..f3da97d
--- /dev/null
+++ b/man/QRSSVSsummary.Rd
@@ -0,0 +1,42 @@
+\name{summary.qrssvs}
+\alias{summary.qrssvs}
+\alias{print.summary.qrssvs}
+\title{Summarising the results of quantile regression stochastic search variable selection (QR-SSVS).}
+
+\description{This function produces a table of predictors and their associated marginal posterior probability of inclusion. It also returns the median probability model (see the details section).}
+
+\usage{\method{summary}{qrssvs}(object, \dots)}
+
+\arguments{
+ \item{object}{An object of class \code{qrssvs}. Typically this will be the \code{gamma} component of the list returned by \code{SSVSquantreg}.}
+
+ \item{\dots}{Further arguments.}
+}
+
+\details{The median probability model is defined to be the model that contains any predictor with marginal posterior probability greater than or equal to 0.5. If the goal is to select a single model e.g. for prediction, Barbieri and Berger (2004) recommend the median probability model. In some cases, this will coincide with the maximum probability model.}
+
+\author{Craig Reed}
+
+\references{
+
+ Maria M. Barbieri, and James O. Berger (2004). "Optimal predictive model selection". \emph{Annals of Statistics}, 32, 870-897.
+}
+
+\examples{
+\dontrun{
+set.seed(1)
+epsilon<-rnorm(100)
+set.seed(2)
+x<-matrix(rnorm(1000),100,10)
+y<-x[,1]+x[,10]+epsilon
+qrssvs<-SSVSquantreg(y~x)
+summary(qrssvs$gamma)
+}
+}
+
+\keyword{models}
+
+\seealso{
+\code{\link[MCMCpack]{SSVSquantreg}},
+\code{\link[MCMCpack]{mptable}},
+\code{\link[MCMCpack]{topmodels}}}
diff --git a/man/SSVSquantreg.Rd b/man/SSVSquantreg.Rd
new file mode 100644
index 0000000..a24ef01
--- /dev/null
+++ b/man/SSVSquantreg.Rd
@@ -0,0 +1,136 @@
+\name{SSVSquantreg}
+
+\alias{SSVSquantreg}
+
+\title{ Stochastic search variable selection for quantile regression }
+
+\description{
+ This function uses stochastic search to select promising regression models at a fixed quantile \eqn{\tau}{tau}.
+ Indicator variables \eqn{\gamma}{gamma} are used to represent whether a predictor is included in the model or not.
+The user supplies the data and the prior distribution on the model size.
+ A list is returned containing the posterior sample of \eqn{\gamma}{gamma} and the associated regression parameters \eqn{\beta}{beta}.}
+
+\usage{
+SSVSquantreg(formula, data = NULL, tau = 0.5, include=NULL, burnin = 1000,
+ mcmc = 10000, thin = 1, verbose = 0, seed = sample(1:1000000,1),
+ pi0a0 = 1, pi0b0 = 1, ...)
+}
+
+\arguments{
+
+ \item{formula}{ Model formula. }
+
+ \item{data}{ Data frame. }
+
+ \item{tau}{The quantile of interest. Must be between 0 and 1. The default value of 0.5 corresponds to median regression model selection.}
+
+ \item{include}{The predictor(s) that should definitely appear in the model. Can be specified by name, or their position in the formula (taking into account the intercept).}
+
+ \item{burnin}{ The number of burn-in iterations for the sampler. }
+
+ \item{mcmc}{ The number of MCMC iterations after burnin. }
+
+ \item{thin}{The thinning interval used in the simulation. The number of
+ MCMC iterations must be divisible by this value.}
+
+ \item{verbose}{ A switch which determines whether or not the progress of
+ the sampler is printed to the screen. If \code{verbose} is greater
+ than 0 the iteration number, the most recently sampled values of \eqn{\gamma}{gamma} and the associated values of \eqn{\beta}{beta} are printed to
+ the screen every \code{verbose}th iteration. }
+
+ \item{seed}{ The seed for the random number generator. If NA, the Mersenne
+ Twister generator is used with default seed 12345; if an integer is
+ passed it is used to seed the Mersenne twister. The default value for this argument is a random integer between 1 and 1,000,000.
+ This default value ensures that if the function is used again with a different value of \eqn{\tau}{tau}, it is extremely unlikely that the seed will be identical. The user can also
+ pass a list of length two to use the L'Ecuyer random number generator,
+ which is suitable for parallel computation. The first element of the
+ list is the L'Ecuyer seed, which is a vector of length six or NA (if NA
+ a default seed of \code{rep(12345,6)} is used). The second element of
+ list is a positive substream number. See the MCMCpack
+ specification for more details. }
+
+ \item{pi0a0, pi0b0}{Hyperparameters of the beta prior on \eqn{\pi_0}{pi_0}, the prior probability of including a predictor. Default values of (1,1) are equivalent to a uniform distribution.}
+
+ \item{\dots}{ Further arguments }
+}
+
+\value{
+ A list containing:
+
+\item{gamma}{The posterior sample of \eqn{\gamma}{gamma}. This has associated summary and plot methods.}
+
+\item{beta}{The posterior sample of the associated regression parameters \eqn{\beta}{beta}. This can be analysed with functions from the coda package.}
+
+ }
+
+\details{
+
+ \code{SSVSquantreg} implements stochastic search variable selection over a set of potential predictors to obtain promising models.
+ The models considered take the following form:
+ \deqn{Q_{\tau}(y_i|x_{i\gamma}) = x_{i\gamma} ' \beta_{\gamma},}{Q_tau(y_i|x_igamma) = x_igamma'beta_gamma,}
+ where \eqn{Q_{\tau}(y_i|x_{i\gamma})}{Q_tau(y_i|x_igamma)} denotes the conditional \eqn{\tau}{tau}th quantile of \eqn{y_i}{y_i} given \eqn{x_{i\gamma}}{x_igamma},
+ \eqn{x_{i\gamma} }{x_igamma} denotes \eqn{x_i}{x_i} with those predictors \eqn{x_{ij}}{x_ij} for which \eqn{\gamma_j=0}{gamma_j=0} removed
+and \eqn{\beta_{\gamma}}{beta_gamma} denotes the model specific regression parameters.
+
+ The likelihood is formed based on the assumption of independent asymmetric Laplace distributions
+ on the \eqn{y_i}{y_i} with skewness parameter \eqn{\tau}{tau} and location parameters \eqn{ x_{i\gamma} ' \beta_{\gamma}}{x_igamma'beta_gamma}. This assumption ensures that the likelihood function is maximised by the \eqn{\tau}{tau}th conditional quantile of the response variable.
+
+ The prior on each \eqn{\beta_j}{beta_j} is
+\deqn{(1-\gamma_j)\delta_0+\gamma_j\mbox{Cauchy}(0,1),}{(1-gamma_j)delta_0+gamma_jCauchy(0,1),}
+where \eqn{\delta_0}{delta_0} denotes a degenerate distribution with all mass at 0.
+A standard Cauchy distribution is chosen conditional on \eqn{\gamma_j=1}{gamma_j=1}.
+ This allows for a wider range of nonzero values of \eqn{\beta_j}{beta_j} than a standard Normal distribution, improving the robustness of the method.
+ Each of the indicator variables \eqn{\gamma_j}{gamma_j} is independently assigned a Bernoulli prior, with prior probability of inclusion \eqn{\pi_0}{pi_0}.
+ This in turn is assigned a beta distribution, resulting in a beta-binomial prior on the model size.
+ The user can supply the hyperparameters for the beta distribution.
+Starting values are randomly generated from the prior distribution.
+
+It is recommended to standardise any non-binary predictors in order to compare these predictors on the same scale.
+This can be achieved using the \code{scale} function.
+
+If it is certain that a predictor should be included, all predictors specified are brought to the first positions for computational convenience. The regression parameters associated with these predictors are given independent improper priors. Users may notice a small speed advantage if they specify the predictors that they feel certain should appear in the model, particularly for large models with a large number of observations.
+}
+
+\author{ Craig Reed}
+
+\references{
+ Craig Reed, David B. Dunson and Keming Yu. 2010. "Bayesian Variable Selection for Quantile Regression" Technical Report.
+
+ Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007.
+ \emph{Scythe Statistical Library 1.2.} \url{http://scythe.wustl.edu}.
+
+ Keming Yu and Jin Zhang. 2005. "A Three Parameter Asymmetric Laplace Distribution and it's extensions."
+\emph{Communications in Statistics - Theory and Methods}, 34, 1867-1879.
+
+ Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2002.
+ \emph{Output Analysis and Diagnostics for MCMC (CODA)}.
+ \url{http://www-fis.iarc.fr/coda/}.}
+
+\examples{
+\dontrun{
+
+set.seed(1)
+epsilon<-rnorm(100)
+set.seed(2)
+x<-matrix(rnorm(1000),100,10)
+y<-x[,1]+x[,10]+epsilon
+qrssvs<-SSVSquantreg(y~x)
+model.50pc<-SSVSquantreg(y~x)
+model.90pc<-SSVSquantreg(y~x,tau=0.9)
+summary(model.50pc) ## Intercept not in median probability model
+summary(model.90pc) ## Intercept appears in median probability model
+}
+}
+
+\keyword{models}
+
+\seealso{
+\code{\link[MCMCpack]{MCMCquantreg}},
+\code{\link[MCMCpack]{summary.qrssvs}},
+\code{\link[MCMCpack]{plot.qrssvs}},
+\code{\link[MCMCpack]{mptable}},
+\code{\link[MCMCpack]{topmodels}},
+\code{\link[base]{scale}},
+\code{\link[quantreg]{rq}}}
+
+
diff --git a/man/mptable.Rd b/man/mptable.Rd
new file mode 100644
index 0000000..f3a4ecd
--- /dev/null
+++ b/man/mptable.Rd
@@ -0,0 +1,33 @@
+\name{mptable}
+\alias{mptable}
+
+\title{Calculate the marginal posterior probabilities of predictors being included in a quantile regression model.}
+
+\description{This function extracts the marginal probability table produced by \code{summary.qrssvs}.}
+
+\usage{mptable(qrssvs)}
+
+\arguments{
+ \item{qrssvs}{An object of class \code{qrssvs}. Typically this will be the \code{gamma} component of the list returned by \code{SSVSquantreg}.}
+}
+
+\value{A table with the predictors listed together with their posterior marginal posterior probability of inclusion.}
+
+\author{Craig Reed}
+
+\examples{
+\dontrun{
+set.seed(1)
+epsilon<-rnorm(100)
+set.seed(2)
+x<-matrix(rnorm(1000),100,10)
+y<-x[,1]+x[,10]+epsilon
+qrssvs<-SSVSquantreg(y~x)
+mptable(qrssvs$gamma)
+}
+}
+
+\keyword{models}
+
+\seealso{
+\code{\link[MCMCpack]{SSVSquantreg}}}
diff --git a/man/topmodels.Rd b/man/topmodels.Rd
new file mode 100644
index 0000000..b6b1ea8
--- /dev/null
+++ b/man/topmodels.Rd
@@ -0,0 +1,42 @@
+\name{topmodels}
+\alias{topmodels}
+
+\title{Shows an ordered list of the most frequently visited models sampled during quantile regression stochastic search variable selection (QR-SSVS).}
+
+\description{Given output from quantile regression stochastic search variable selection, this function returns a table of the 'best' models together with their associated empirical posterior probability.}
+
+\usage{
+topmodels(qrssvs, nmodels=5, abbreviate=FALSE, minlength=3)
+}
+
+\arguments{
+
+ \item{qrssvs}{An object of class \code{qrssvs}. Typically this will be the \code{gamma} component of the list returned by \code{SSVSquantreg}.}
+
+ \item{nmodels}{The number of models to tabulate.}
+
+ \item{abbreviate}{Logical: should the names of the predictors be abbreviated?}
+
+ \item{minlength}{If \code{abbreviate} is set to \code{TRUE}, the minimum length of the abbreviations. }
+}
+
+\value{A table with the models and their associated posterior probability. The models are arranged in descending order of probability.}
+
+\author{Craig Reed}
+
+\examples{
+\dontrun{
+set.seed(1)
+epsilon<-rnorm(100)
+set.seed(2)
+x<-matrix(rnorm(1000),100,10)
+y<-x[,1]+x[,10]+epsilon
+qrssvs<-SSVSquantreg(y~x)
+topmodels(qrssvs$gamma)
+}
+}
+
+\keyword{models}
+
+\seealso{
+\code{\link[MCMCpack]{SSVSquantreg}}}
diff --git a/src/MCMCdynamicIRT1d-b.cc b/src/MCMCdynamicIRT1d-b.cc
index 48422e6..6f2be70 100644
--- a/src/MCMCdynamicIRT1d-b.cc
+++ b/src/MCMCdynamicIRT1d-b.cc
@@ -92,8 +92,10 @@ void MCMCdynamicIRT1d_b_impl(rng<RNGTYPE>& stream,
const int tot_iter = *burnin + *mcmc;
const int nsamp = *mcmc / *thin;
//sparse matrix of latent outcome variables
- double Z[*nsubj][*nitems];
+ double** Z;
+ Z = new double*[*nsubj];
for (int s=0; s<*nsubj; ++s){
+ Z[s] = new double[*nitems];
for (int i=0; i<*nitems; ++i){
Z[s][i] = -999;
}
@@ -624,8 +626,10 @@ void MCMCdynamicIRT1d_b_impl(rng<RNGTYPE>& stream,
} // end iter loop
-
-
+ for (int s=0; s<*nsubj; ++s)
+ delete[] Z[s];
+ delete[] Z;
+
} // end MCMCdynamicIRT1d_impl function
diff --git a/src/MCMCdynamicIRT1d.cc b/src/MCMCdynamicIRT1d.cc
index 8768205..bd000ef 100644
--- a/src/MCMCdynamicIRT1d.cc
+++ b/src/MCMCdynamicIRT1d.cc
@@ -92,8 +92,10 @@ void MCMCdynamicIRT1d_impl(rng<RNGTYPE>& stream,
const int tot_iter = *burnin + *mcmc;
const int nsamp = *mcmc / *thin;
//sparse matrix of latent outcome variables
- double Z[*nsubj][*nitems];
+ double** Z;
+ Z = new double*[*nsubj];
for (int s=0; s<*nsubj; ++s){
+ Z[s] = new double[*nitems];
for (int i=0; i<*nitems; ++i){
Z[s][i] = -999;
}
@@ -612,8 +614,10 @@ void MCMCdynamicIRT1d_impl(rng<RNGTYPE>& stream,
} // end iter loop
-
-
+ for (int s=0; s<*nsubj; ++s)
+ delete[] Z[s];
+ delete[] Z;
+
} // end MCMCdynamicIRT1d_impl function
diff --git a/src/MCMCfcds.h b/src/MCMCfcds.h
index bcfe162..6cae159 100644
--- a/src/MCMCfcds.h
+++ b/src/MCMCfcds.h
@@ -87,87 +87,24 @@ NormIGregress_sigma2_draw (const Matrix <> &X, const Matrix <> &Y,
return stream.rigamma (c_post, d_post);
}
-// linear regression with Laplace errors beta draw
+// Bayesian quantile (including median) regression beta draw
// (multivariate Normal prior)
-// regression model is y = X * beta + epsilon, epsilon ~ Laplace(0,sigma)
-// b0 is the prior mean of beta
-// B0 is the prior precision (the inverse variance) of beta
-template <typename RNGTYPE>
-Matrix<double>
-LaplaceNormregress_beta_draw (const Matrix<>& X, const Matrix<>& Y, const Matrix<>& weights,
- const Matrix<>& b0, const Matrix<>& B0, double sigma,
- rng<RNGTYPE>& stream)
-{
-
- const unsigned int k = X.cols();
- const unsigned int n_obs = X.rows();
- const double one_over_two_sigma = 1.0/(2.0*sigma);
- Matrix<> XtwX(k,k,false);
- Matrix<> XtwY(k,1,false);
- double temp_x = 0.0;
- double temp_y = 0.0;
-
-//Calculate XtwY, where w denotes a diagonal matrix with the augmented data (weights) on the diagonal
- for (unsigned int i=0; i<k; ++i){
- for (unsigned int m=0; m<n_obs; ++m){
- temp_y = temp_y + weights(m)*X(m,i)*Y(m);
- }
- XtwY(i) = temp_y;
- temp_y = 0.0;
- }
-//Calculate XtwX
- for (unsigned int i=0; i<k; ++i){
- for (unsigned int j=0; j<(i+1); ++j){
- for (unsigned int m=0; m<n_obs; ++m){
- temp_x = temp_x + weights(m)*X(m,i)*X(m,j);
- }
- XtwX(i,j) = temp_x;
- XtwX(j,i) = temp_x;
- temp_x = 0.0;
- }
- }
-
- const Matrix<> var_matrix_beta = invpd(B0+one_over_two_sigma*XtwX);
- const Matrix<> C = cholesky(var_matrix_beta);
- const Matrix<> betahat = var_matrix_beta*gaxpy(B0,b0,one_over_two_sigma*XtwY);
-
- return( gaxpy(C, stream.rnorm(k,1, 0, 1), betahat) );
-}
-
-// linear regression with Laplace errors sigma draw
-// (inverse-Gamma prior)
-// regression model is y = X * beta + epsilon, epsilon ~ Laplace(0,sigma)
-// c0/2 is the prior shape parameter for sigma
-// d0/2 is the prior scale parameter for sigma
-template <typename RNGTYPE>
-double
-LaplaceIGammaregress_sigma_draw (const Matrix <> &abse, double c0, double d0,
- rng<RNGTYPE>& stream)
-
-{
- const double c_post = 0.5*c0 + abse.rows();
- const double d_post = 0.5*(d0 + sum(abse));
-
- return stream.rigamma (c_post, d_post);
-
-}
-
-// linear regression with Asymmetric Laplace errors beta draw
-// (multivariate Normal prior)
-// regression model is y = X * beta + epsilon, epsilon ~ ALaplace(0,sigma,tau)
-// b0 is the prior mean of beta
-// B0 is the prior precision (the inverse variance) of beta
+//
+// b0 is the prior mean of beta(tau)
+// B0 is the prior precision (the inverse variance) of beta(tau)
template <typename RNGTYPE>
Matrix<double>
ALaplaceNormregress_beta_draw (double tau, const Matrix<>& X, const Matrix<>& Y, const Matrix<>& weights,
- const Matrix<>& b0, const Matrix<>& B0, double sigma,
+ const Matrix<>& b0, const Matrix<>& B0,
rng<RNGTYPE>& stream)
{
const unsigned int k = X.cols();
const unsigned int n_obs = X.rows();
- const double one_over_two_sigma = 1.0/(2.0*sigma);
- const Matrix<> U = Y - (1.0-2.0*tau)*pow(weights,-1.0);
+ Matrix<> U(Y);
+if (tau!=0.5){
+ U -= (1.0-2.0*tau)*weights;
+}
Matrix<> XtwX(k,k,false);
Matrix<> XtwU(k,1,false);
double temp_x = 0.0;
@@ -176,16 +113,17 @@ ALaplaceNormregress_beta_draw (double tau, const Matrix<>& X, const Matrix<>& Y,
//Calculate XtwU where w denotes a diagonal matrix with the augmented data (weights) on the diagonal
for (unsigned int i=0; i<k; ++i){
for (unsigned int m=0; m<n_obs; ++m){
- temp_u = temp_u + weights(m)*X(m,i)*U(m);
+ temp_u += X(m,i)*U(m)/weights(m);
}
XtwU(i) = temp_u;
temp_u = 0.0;
}
+
//Calculate XtwX
for (unsigned int i=0; i<k; ++i){
for (unsigned int j=0; j<(i+1); ++j){
for (unsigned int m=0; m<n_obs; ++m){
- temp_x = temp_x + weights(m)*X(m,i)*X(m,j);
+ temp_x += X(m,i)*X(m,j)/weights(m);
}
XtwX(i,j) = temp_x;
XtwX(j,i) = temp_x;
@@ -193,50 +131,31 @@ ALaplaceNormregress_beta_draw (double tau, const Matrix<>& X, const Matrix<>& Y,
}
}
- const Matrix<> var_matrix_beta = invpd(B0+one_over_two_sigma*XtwX);
+ const Matrix<> var_matrix_beta = invpd(B0+0.5*XtwX);
const Matrix<> C = cholesky(var_matrix_beta);
- const Matrix<> betahat = var_matrix_beta*gaxpy(B0,b0,one_over_two_sigma*XtwU);
+ const Matrix<> betahat = var_matrix_beta*gaxpy(B0,b0,0.5*XtwU);
return( gaxpy(C, stream.rnorm(k,1, 0, 1), betahat) );
}
-// linear regression with Asymmetric Laplace errors sigma draw
-// (inverse-Gamma prior)
-// regression model is y = X * beta + epsilon, epsilon ~ ALaplace(0,sigma,tau)
-// c0/2 is the prior shape parameter for sigma
-// d0/2 is the prior scale parameter for sigma
-template <typename RNGTYPE>
-double
-ALaplaceIGammaregress_sigma_draw (double tau, const Matrix<> &e, const Matrix <> &abse, double c0, double d0,
- rng<RNGTYPE>& stream)
-
-{
- const double c_post = 0.5*c0 + abse.rows();
- const double d_post = 0.5*(d0 + sum(abse)+(2.0*tau-1.0)*sum(e));
-
- return stream.rigamma (c_post, d_post);
-
-}
-
// This function draws from the full conditional distribution of the latent random variables (weights) under quantile regression (including median regression) and returns a column vector of those weights.
template <typename RNGTYPE>
Matrix<double>
-ALaplaceIGaussregress_weights_draw (const Matrix <> &abse, double sigma,
+ALaplaceIGaussregress_weights_draw (const Matrix <> &abse,
rng<RNGTYPE>& stream)
{
- const double lambda = 1.0/(2.0*sigma);
- const Matrix<double> nu_params = pow(abse,-1.0);
+ const Matrix<double> nu_params = pow(abse,-1.0);
Matrix<> w(abse);
- unsigned int n_obs = abse.rows();
+ const unsigned int n_obs = abse.rows();
// The inverse Gaussian distribution
for (unsigned int i=0; i<n_obs; ++i){
double chisq = stream.rchisq(1);
double nu = nu_params(i);
- double smallroot = nu/(2.0*lambda)*(nu*chisq+2.0*lambda-std::sqrt(nu*nu*chisq*chisq+4.0*nu*lambda*chisq));
+ double smallroot = nu*(nu*chisq+1.0-std::sqrt(nu*nu*chisq*chisq+2.0*nu*chisq));
unsigned int q = stream.rbern(nu/(nu+smallroot));
if (q == 1){
w(i) = smallroot;
@@ -245,9 +164,267 @@ ALaplaceIGaussregress_weights_draw (const Matrix <> &abse, double sigma,
w(i) = nu*nu/smallroot;
}
}
- return w;
+ return(pow(w,-1.0));
}
+/////////////////////////////////////////////////////
+
+// Functions for the quantile regression stochastic search variable selection
+
+struct COV_TRIAL{
+Matrix<> Cnew;
+bool newtrial;
+double logdetminhalf;
+};
+
+// updating the indicator variable corresponding to whether a
+// covariate is included in the model, given that it was
+// previously absent
+template <typename RNGTYPE>
+COV_TRIAL
+QR_SSVS_covariate_trials_draw_absent(const Matrix<>& C, const Matrix<>& X_gamma, const Matrix<>& U,
+const Matrix<>& newXcol, unsigned int row_index, const Matrix<>& weights, double pi0, double newlambda, double logolddetminhalf, rng<RNGTYPE>& stream){
+ const unsigned int n_obs = U.rows();
+ const unsigned int k = C.rows();
+
+//Calculate new row required to update the Cholesky decomposition
+ Matrix<> XUXnewtXnew(k+1,1,false);
+ double temp_xux1 = 0.0;
+ double temp_xux2 = 0.0;
+ double temp_xux3 = 0.0;
+
+//Calculate XUXnewtXnew
+ for (unsigned int i=0; i<k-1; ++i){
+ for (unsigned int m=0; m<n_obs; ++m){
+ temp_xux1 += X_gamma(m,i)*newXcol(m)/weights(m);
+ }
+ XUXnewtXnew(i) = 0.5*temp_xux1;
+ temp_xux1 = 0.0;
+ }
+ for (unsigned int m=0; m<n_obs; ++m){
+temp_xux2 += U(m)*newXcol(m)/weights(m);
+temp_xux3 += newXcol(m)*newXcol(m)/weights(m);
+}
+XUXnewtXnew(k-1) = 0.5*temp_xux2;
+XUXnewtXnew(k) = 0.5*temp_xux3+newlambda;
+
+//Obtain the Cholesky Decomposition of the new matrix formed by adding newXcol onto the right hand side
+
+Matrix<> z(k,1,false);
+for (unsigned int i = 0; i < k; ++i) {
+ double sum = 0;
+ for (unsigned int j = 0; j < i; ++j) {
+ sum += C(i,j) * z(j);
+ }
+ z(i) = (XUXnewtXnew(i) - sum) / C(i, i);
+ }
+double rho = std::sqrt(XUXnewtXnew(k)-crossprod(z)(0));
+
+Matrix<> Cnew(k+1, k+1, true, 0.0);
+Cnew(0,0,k-1,k-1) = C;
+Cnew(k,0,k,k-1) = z;
+Cnew(k,k) = rho;
+
+// Permuting the Cholesky decomposition so that it corresponds to what would be obtained if the X matrix included the covariate
+
+ Matrix<> temp(Cnew);
+if (row_index != 0){
+ temp(0,0,row_index-1,k) = Cnew(0,0,row_index-1,k);
+}
+ temp(row_index,_) = Cnew(k,_);
+ temp(row_index+1,0,k,k) = Cnew(row_index,0,k-1,k);
+
+// Givens rotations
+
+ Matrix<> Q(2,2,false);
+ for (unsigned int i=k; i>row_index; --i)
+ {
+ double two_norm = std::sqrt(temp(row_index,i)*temp(row_index,i)
+ +temp(row_index,i-1)*temp(row_index,i-1));
+ Q(0,0) = temp(row_index,i-1)/two_norm;
+ Q(1,0) = temp(row_index,i)/two_norm;
+ Q(1,1) = Q(0,0);
+ Q(0,1) = -1.0*Q(1,0);
+ if (i!=k){
+ temp(i+1,i-1,k,i) = temp(i+1,i-1,k,i) * Q;
+ }
+ double temp2 = temp(i,i-1);
+ temp(i,i-1) = Q(0,0)*temp2;
+ temp(i,i) = Q(0,1)*temp2;
+ if (temp(i,i) < 0){
+ temp(i,i,k,i) = -1.0*temp(i,i,k,i);
+ }
+ temp(row_index,i-1) = two_norm;
+ temp(row_index,i) = 0.0;
+ }
+ Cnew=temp;
+
+//Work out -0.5*log(det(Cnew'Cnew))
+ double lognewdetminhalf = 0.0;
+ for (unsigned int i=0; i<k; ++i){
+ lognewdetminhalf -= std::log(Cnew(i,i));
+}
+
+
+
+double log_g0 = logolddetminhalf-0.5*C(k-1,k-1)*C(k-1,k-1);
+double log_g1 = 0.5*std::log(newlambda)+lognewdetminhalf-0.5*Cnew(k,k)*Cnew(k,k);
+
+ double log_ratio = log_g0+std::log(1.0-pi0)-log_g1-std::log(pi0);
+ double success_prob = 1.0/(1.0+std::exp(log_ratio));
+ bool new_covariate_trial = stream.rbern(success_prob);
+COV_TRIAL result;
+result.newtrial = new_covariate_trial;
+ if (new_covariate_trial == false){
+ result.Cnew = C;
+ result.logdetminhalf = logolddetminhalf;
+ }
+ else {
+ result.Cnew = Cnew;
+ result.logdetminhalf = lognewdetminhalf;
+ }
+return result;
+}
+
+// updating the indicator variable corresponding to whether a
+// covariate is included in the model, given that it was
+// previously present in the model
+
+template <typename RNGTYPE>
+COV_TRIAL
+QR_SSVS_covariate_trials_draw_present(const Matrix<>& C, unsigned int row_index, unsigned int n_obs, double pi0, double oldlambda, double logolddetminhalf, rng<RNGTYPE>& stream){
+ unsigned int k = C.rows();
+
+// Permuting the Cholesky decomposition so that it corresponds to what would be obtained if the X matrix had the covariate in the final column
+
+ Matrix<> temp(C);
+if (row_index != 0){
+ temp(0,0,row_index-1,k-1) = C(0,0,row_index-1,k-1);
+}
+ temp(k-1,_) = C(row_index,_);
+ temp(row_index,0,k-2,k-1) = C(row_index+1,0,k-1,k-1);
+
+// Givens rotations
+
+ Matrix<> Q(2,2,false);
+ for (unsigned int i=row_index; i<k-1; ++i)
+ {
+ double two_norm = std::sqrt(temp(i,i)*temp(i,i)
+ +temp(i,i+1)*temp(i,i+1));
+ Q(0,0) = temp(i,i)/two_norm;
+ Q(1,0) = temp(i,i+1)/two_norm;
+ Q(1,1) = Q(0,0);
+ Q(0,1) = -1.0*Q(1,0);
+ if (i!=k-2){
+ temp(i+1,i,k-2,i+1) = temp(i+1,i,k-2,i+1) * Q;
+ }
+ double temp2 = temp(k-1,i);
+ temp(k-1,i) = Q(0,0)*temp2;
+ temp(k-1,i+1) = Q(0,1)*temp2;
+ temp(i,i) = two_norm;
+ temp(i,i+1) = 0.0;
+ }
+ if (temp(k-1,k-1) < 0){
+ temp(k-1,k-1) = -1.0*temp(k-1,k-1);
+}
+
+Matrix<> Cnew = temp(0,0,k-2,k-2);
+
+// Work out -1/2*log(det(Cnew'Cnew))
+ double lognewdetminhalf = 0.0;
+ for (unsigned int i=0; i<k-2; ++i){
+ lognewdetminhalf -= std::log(Cnew(i,i));
+}
+
+
+ double log_g0 = lognewdetminhalf-0.5*Cnew(k-2,k-2)*Cnew(k-2,k-2);
+ double log_g1 = 0.5*std::log(oldlambda)+logolddetminhalf-0.5*C(k-1,k-1)*C(k-1,k-1);
+
+
+double log_ratio = log_g0+std::log(1.0-pi0)-log_g1-std::log(pi0);
+ double success_prob = 1.0/(1.0+std::exp(log_ratio));
+ bool new_covariate_trial = stream.rbern(success_prob);
+COV_TRIAL result;
+result.newtrial = new_covariate_trial;
+ if (new_covariate_trial == false){
+ result.Cnew = Cnew;
+ result.logdetminhalf = lognewdetminhalf;
+ }
+ else {
+ result.Cnew = C;
+ result.logdetminhalf = logolddetminhalf;
+ }
+return result;
+}
+
+// update betas using Cholesky decomposition
+template <typename RNGTYPE>
+Matrix<>
+QR_SSVS_beta_draw(const Matrix<>& C, rng<RNGTYPE>& stream){
+ unsigned int k = C.rows();
+ Matrix<> standnorm = stream.rnorm(k-1,1,0,1);
+ Matrix<> z(k-1,1,false);
+ z = t(C(k-1,0,k-1,k-2));
+ Matrix<> Q = z+standnorm*std::sqrt(2.0);
+ Matrix<> result(k-1,1,false);
+ for (int i = k-2; i >= 0; --i) {
+ double sum = 0;
+ for (unsigned int j = i+1; j < k-1; ++j) {
+ sum += C(j,i) * result(j);
+ }
+ result(i) = (Q(i) - sum) / C(i, i);
+ }
+ return result;
+}
+
+//hyperparameter pi0 updating
+
+template <typename RNGTYPE>
+double
+QR_SSVS_pi0_draw(unsigned int n_uncert_cov, unsigned int tot_n_uncert_cov,
+double pi0a0, double pi0b0, rng<RNGTYPE>& stream){
+ double pi0a1 = pi0a0 + n_uncert_cov;
+ double pi0b1 = pi0b0 + tot_n_uncert_cov - n_uncert_cov;
+ return(stream.rbeta(pi0a1,pi0b1));
+}
+
+//update latent lambdas
+
+template<typename RNGTYPE>
+Matrix<double>
+QR_SSVS_lambda_draw(const Matrix<>& beta_red, const Matrix<>& gamma, unsigned int tot_n_cov, unsigned int n_cert_cov, rng<RNGTYPE>& stream)
+ {
+unsigned int n_uncert_cov = tot_n_cov - n_cert_cov;
+Matrix<> newlambda(n_uncert_cov,1,false);
+
+for (unsigned int i=n_cert_cov; i<tot_n_cov; ++i){
+
+unsigned int j = i-n_cert_cov;
+
+if (gamma(i) == true){
+
+ unsigned int col_index = n_cert_cov;
+ //obtain column index of betas
+ for (unsigned int m=n_cert_cov; m<i; ++m){
+ if (gamma(m) == true){
+ ++col_index;
+ }
+ }
+
+newlambda(j) = stream.rexp(0.5*(1.0+beta_red(col_index)*beta_red(col_index)));
+}
+
+else {
+newlambda(j) = stream.rexp(0.5);
+}
+
+}
+ return newlambda;
+ }
+
+
+/////////////////////////////////////////////////////
+
// update latent data for standard item response models
// only works for 1 dimensional case
template <typename RNGTYPE>
diff --git a/src/MCMCirtHier1d.cc b/src/MCMCirtHier1d.cc
index bca391f..7370b37 100644
--- a/src/MCMCirtHier1d.cc
+++ b/src/MCMCirtHier1d.cc
@@ -86,7 +86,7 @@ double irt_W_update(Matrix<>& Z, const Matrix<>& X, const Matrix<>& theta,
//Rprintf("\nRSS: %5f alpha0: %5f alpha1: %5f\n",RSS,alpha,alpha1);
return(std::sqrt(alpha1/alpha));
}
-
+/*
template <typename RNGTYPE>
void hirt_level0_metrop (Matrix<>& theta, Matrix<>& thetahat,
const Matrix<>& Z,
@@ -101,7 +101,8 @@ void hirt_level0_metrop (Matrix<>& theta, Matrix<>& thetahat,
//exp( prop - cur)
//runif & accept if > ratio
//note the acceptance
-}
+ }
+*/
/* MCMCirt1d implementation. */
@@ -267,9 +268,11 @@ void MCMCirtHier1d_impl (rng<RNGTYPE>& stream,
double sigma2fcdmean = sigma2fcdsum / static_cast<double>(tot_iter);
const double sig2_inv = 1.0 / sigma2star;
const Matrix<> sig_beta = invpd (B0 + XpX * sig2_inv);
- const Matrix<> betahat = sig_beta * gaxpy(B0, b0, XpY*sig2_inv);
- double logbetafcd = lndmvn(betastar, betahat, sig_beta);
- if(L==1) { logbetafcd = lndnorm(betastar[0], betahat[0],sig_beta[0]);}
+ Matrix<> betahat = sig_beta * gaxpy(B0, b0, XpY*sig2_inv);
+ double logbetafcd = 0.0;
+ for (unsigned int i=0; i<L; ++i) {
+ logbetafcd += lndnorm(betastar[i], betahat[i], sig_beta[i]);
+ }
Rprintf("logPost: beta %10.5f sigma2 %10.5f\n",logbetafcd,sigma2fcdmean);
// calculate loglikelihood at (betastar, sigma2star)
double sigmastar = sqrt(sigma2star);
diff --git a/src/MCMCmedreg.cc b/src/MCMCmedreg.cc
deleted file mode 100644
index 81c52fb..0000000
--- a/src/MCMCmedreg.cc
+++ /dev/null
@@ -1,138 +0,0 @@
-// MCMCmedreg.cc is a program that simulates draws from the
-// posterior distribution of a linear regression model with
-// errors that come from a Laplace (double exponential) distribution.
-//
-// The initial version of this file was generated by the
-// auto.Scythe.call() function in the MCMCpack R package
-// written by:
-//
-// Andrew D. Martin
-// Dept. of Political Science
-// Washington University in St. Louis
-// admartin at wustl.edu
-//
-// Kevin M. Quinn
-// Dept. of Government
-// Harvard University
-// kevin_quinn at harvard.edu
-//
-// This software is distributed under the terms of the GNU GENERAL
-// PUBLIC LICENSE Version 2, June 1991. See the package LICENSE
-// file for more information.
-//
-// Copyright (C) 2009 Andrew D. Martin and Kevin M. Quinn
-//
-// This file was initially generated on Wed Apr 1 11:39:12 2009
-// The function was rewritten by:
-//
-// Craig Reed
-// Department of Mathematical Sciences
-// Brunel University
-// craig.reed at brunel.ac.uk
-//
-// CR 7/4/09 [Rewritten using templates]
-
-#ifndef MCMCMEDREG_CC
-#define MCMCMEDREG_CC
-
-#include "MCMCrng.h"
-#include "MCMCfcds.h"
-#include "matrix.h"
-#include "distributions.h"
-#include "stat.h"
-#include "la.h"
-#include "ide.h"
-#include "smath.h"
-#include "rng.h"
-
-#include <R.h> // needed to use Rprintf()
-#include <R_ext/Utils.h> // needed to allow user interrupts
-
-using namespace std;
-using namespace scythe;
-
-/* MCMCmedreg implementation. Takes Matrix<> reference which it
- * fills with the posterior.
- */
-template <typename RNGTYPE>
-void MCMCmedreg_impl (rng<RNGTYPE>& stream, const Matrix<>& Y,
- const Matrix<>& X, Matrix<>& beta, const Matrix<>& b0,
- const Matrix<>& B0, double c0, double d0,
- unsigned int burnin, unsigned int mcmc, unsigned int thin,
- unsigned int verbose,
- Matrix<>& result)
-{
- // define constants
- const unsigned int tot_iter = burnin + mcmc; //total iterations
- const unsigned int nstore = mcmc / thin; // number of draws to store
- const unsigned int k = X.cols ();
- const unsigned int n_obs = X.rows();
-
- // storage matrices
- Matrix<> betamatrix (k, nstore);
- Matrix<> sigmamatrix (1, nstore);
-
- // Gibbs sampler
- unsigned int count = 0;
- for (unsigned int iter = 0; iter < tot_iter; ++iter) {
- Matrix<> abse = fabs(gaxpy(X, (-1*beta), Y));
- double sigma = LaplaceIGammaregress_sigma_draw (abse, c0, d0, stream);
- Matrix<> weights = ALaplaceIGaussregress_weights_draw (abse, sigma, stream);
- beta = LaplaceNormregress_beta_draw (X, Y, weights, b0, B0, sigma, stream);
-
- // store draws in storage matrix (or matrices)
- if (iter >= burnin && (iter % thin == 0)) {
- sigmamatrix (0, count) = sigma;
- betamatrix(_, count) = beta;
- ++count;
- }
-
- // print output to stdout
- if(verbose > 0 && iter % verbose == 0) {
- Rprintf("\n\nMCMCmedreg iteration %i of %i \n",
- (iter+1), tot_iter);
- Rprintf("beta = \n");
- for (unsigned int j=0; j<k; ++j)
- Rprintf("%10.5f\n", beta(j));
- Rprintf("sigma = %10.5f\n", sigma);
- }
-
- R_CheckUserInterrupt(); // allow user interrupts
-
- } // end MCMC loop
-
- result = cbind(t(betamatrix), t(sigmamatrix));
- } // end MCMCmedreg_impl
-
-extern "C" {
- void MCMCmedreg(double *sampledata, const int *samplerow,
- const int *samplecol, const double *Ydata, const int *Yrow,
- const int *Ycol, const double *Xdata, const int *Xrow,
- const int *Xcol, const int *burnin, const int *mcmc,
- const int *thin, const int *uselecuyer, const int *seedarray,
- const int *lecuyerstream, const int *verbose,
- const double *betastartdata, const int *betastartrow,
- const int *betastartcol, const double *b0data,
- const int *b0row, const int *b0col,
- const double *B0data, const int *B0row,
- const int *B0col, const double *c0, const double *d0)
- {
- // pull together Matrix objects
- Matrix<> Y(*Yrow, *Ycol, Ydata);
- Matrix<> X(*Xrow, *Xcol, Xdata);
- Matrix<> betastart(*betastartrow, *betastartcol, betastartdata);
- Matrix<> b0(*b0row, *b0col, b0data);
- Matrix<> B0(*B0row, *B0col, B0data);
-
- Matrix<> storagematrix;
- MCMCPACK_PASSRNG2MODEL(MCMCmedreg_impl, Y, X, betastart, b0, B0,
- *c0, *d0, *burnin, *mcmc, *thin, *verbose,
- storagematrix);
-
- const unsigned int size = *samplerow * *samplecol;
- for (unsigned int i = 0; i < size; ++i)
- sampledata[i] = storagematrix(i);
- }
-}
-
-#endif
diff --git a/src/MCMCquantreg.cc b/src/MCMCquantreg.cc
index 35ecab3..e1297be 100644
--- a/src/MCMCquantreg.cc
+++ b/src/MCMCquantreg.cc
@@ -1,6 +1,5 @@
// MCMCquantreg.cc is a function that draws from the posterior
-// distribution of a linear regression model with errors that
-// come from an Asymmetric Laplace distribution having pth quantile equal to zero.
+// distribution of a linear regression model with errors having tauth quantile equal to zero.
//
// The initial version of this file was generated by the
// auto.Scythe.call() function in the MCMCpack R package
@@ -32,6 +31,8 @@
// craig.reed at brunel.ac.uk
//
// CR 9/4/09 [Rewritten function using templates]
+//
+// CR 7/12/09 [Placed MCMCmedreg within MCMCquantreg]
#ifndef MCMCQUANTREG_CC
#define MCMCQUANTREG_CC
@@ -56,9 +57,9 @@ using namespace scythe;
* fills with the posterior.
*/
template <typename RNGTYPE>
-void MCMCquantreg_impl (rng<RNGTYPE>& stream, double tau, const Matrix<>& Y,
+void MCMCquantreg_impl (rng<RNGTYPE>& stream, double tau, Matrix<>& Y,
const Matrix<>& X, Matrix<>& beta, const Matrix<>& b0,
- const Matrix<>& B0, double c0, double d0,
+ const Matrix<>& B0,
unsigned int burnin, unsigned int mcmc, unsigned int thin,
unsigned int verbose,
Matrix<>& result)
@@ -71,20 +72,17 @@ void MCMCquantreg_impl (rng<RNGTYPE>& stream, double tau, const Matrix<>& Y,
// storage matrices
Matrix<> betamatrix (k, nstore);
- Matrix<> sigmamatrix (1, nstore);
-
+
// Gibbs sampler
unsigned int count = 0;
for (unsigned int iter = 0; iter < tot_iter; ++iter) {
Matrix<> e = gaxpy(X, (-1*beta), Y);
Matrix<> abse = fabs(e);
- double sigma = ALaplaceIGammaregress_sigma_draw (tau, e, abse, c0, d0, stream);
- Matrix<> weights = ALaplaceIGaussregress_weights_draw (abse, sigma, stream);
- beta = ALaplaceNormregress_beta_draw (tau, X, Y, weights, b0, B0, sigma, stream);
+ Matrix<> weights = ALaplaceIGaussregress_weights_draw (abse, stream);
+ beta = ALaplaceNormregress_beta_draw (tau, X, Y, weights, b0, B0, stream);
- // store draws in storage matrix (or matrices)
+ // store draws in storage matrix
if (iter >= burnin && (iter % thin == 0)) {
- sigmamatrix (0, count) = sigma;
betamatrix(_, count) = beta;
++count;
}
@@ -96,14 +94,13 @@ void MCMCquantreg_impl (rng<RNGTYPE>& stream, double tau, const Matrix<>& Y,
Rprintf("beta = \n");
for (unsigned int j=0; j<k; ++j)
Rprintf("%10.5f\n", beta(j));
- Rprintf("sigma = %10.5f\n", sigma);
}
R_CheckUserInterrupt(); // allow user interrupts
} // end MCMC loop
- result = cbind(t(betamatrix), t(sigmamatrix));
+ result = t(betamatrix);
} // end MCMCquantreg_impl
extern "C" {
@@ -117,7 +114,7 @@ extern "C" {
const int *betastartcol, const double *b0data,
const int *b0row, const int *b0col,
const double *B0data, const int *B0row,
- const int *B0col, const double *c0, const double *d0)
+ const int *B0col)
{
// pull together Matrix objects
Matrix<> Y(*Yrow, *Ycol, Ydata);
@@ -128,7 +125,7 @@ extern "C" {
Matrix<> storagematrix;
MCMCPACK_PASSRNG2MODEL(MCMCquantreg_impl, *tau, Y, X, betastart, b0, B0,
- *c0, *d0, *burnin, *mcmc, *thin, *verbose,
+ *burnin, *mcmc, *thin, *verbose,
storagematrix);
const unsigned int size = *samplerow * *samplecol;
diff --git a/src/SSVSquantreg.cc b/src/SSVSquantreg.cc
new file mode 100644
index 0000000..28907c6
--- /dev/null
+++ b/src/SSVSquantreg.cc
@@ -0,0 +1,304 @@
+// SSVSquantreg.cc is a function that uses stochastic search variable selection
+// to select promising models at a pre-specified quantile.
+//
+// The initial version of this file was generated by the
+// auto.Scythe.call() function in the MCMCpack R package
+// written by:
+//
+// Andrew D. Martin
+// Dept. of Political Science
+// Washington University in St. Louis
+// admartin at wustl.edu
+//
+// Kevin M. Quinn
+// Dept. of Government
+// Harvard University
+// kevin_quinn at harvard.edu
+//
+// This software is distributed under the terms of the GNU GENERAL
+// PUBLIC LICENSE Version 2, June 1991. See the package LICENSE
+// file for more information.
+//
+// Copyright (C) 2009 Andrew D. Martin and Kevin M. Quinn
+//
+// This file was initially generated on June 1 2009
+//
+// The function was rewritten by:
+//
+// Craig Reed
+// Department of Mathematical Sciences
+// Brunel University
+// craig.reed at brunel.ac.uk
+
+#ifndef SSVSQUANTREG_CC
+#define SSVSQUANTREG_CC
+
+#include "MCMCrng.h"
+#include "MCMCfcds.h"
+#include "matrix.h"
+#include "distributions.h"
+#include "stat.h"
+#include "la.h"
+#include "ide.h"
+#include "smath.h"
+#include "rng.h"
+
+#include <R.h> // needed to use Rprintf()
+#include <R_ext/Utils.h> // needed to allow user interrupts
+
+using namespace std;
+using namespace scythe;
+
+struct COV_TRIAL_PREP{
+Matrix<> C;
+Matrix<> U;
+double logdetminhalf;
+};
+
+// Helper function to prepare for the sequence of draws from the covariate indicators
+
+static inline COV_TRIAL_PREP
+QR_SSVS_covariate_trials_prep(const Matrix<>& X_gamma, const Matrix<>& Y, const Matrix<bool>& gamma, const Matrix<>& weights, const Matrix<>& lambda, double tau, unsigned int n_cov, unsigned int q){
+
+const unsigned int n_obs=Y.rows();
+Matrix<> U(Y);
+
+// obtain cholesky decomposition of (Xtilde, utilde)
+ if (tau!=0.5){
+ U -= (1.0-2.0*tau)*weights;
+ }
+ Matrix<> XU(U.rows(),n_cov+1,false);
+ if (n_cov == 0){
+ XU = U;
+ }
+ else{
+ XU = cbind(X_gamma,U);
+ }
+
+ Matrix<> XUtwXU(n_cov+1,n_cov+1,false);
+
+ double temp_xu = 0.0;
+
+//Calculate XUtwXU
+ for (unsigned int i=0; i<n_cov+1; ++i){
+ for (unsigned int j=0; j<(i+1); ++j){
+ for (unsigned int m=0; m<n_obs; ++m){
+ temp_xu += XU(m,i) * XU(m,j) / weights(m);
+ }
+ XUtwXU(i,j) = temp_xu;
+ XUtwXU(j,i) = temp_xu;
+ temp_xu = 0.0;
+ }
+ }
+
+XUtwXU *= 0.5;
+const Matrix<> lambda_gamma = selif(lambda,gamma(q,0,gamma.rows()-1,0));
+
+for (unsigned int j=q; j<n_cov; ++j){
+
+XUtwXU(j,j) += lambda_gamma(j-q);
+}
+
+Matrix<> C = cholesky(XUtwXU);
+
+double logdetminhalf = 0.0;
+// Work out -1/2 log(det(C'C))
+for (unsigned int r=0; r<n_cov; ++r){
+ logdetminhalf -= std::log(C(r,r));
+}
+COV_TRIAL_PREP result;
+result.logdetminhalf = logdetminhalf;
+result.U = U;
+result.C = C;
+return result;
+}
+
+/* SSVSquantreg implementation. Takes Matrix<> reference which it
+ * fills with the posterior of the indicator variables.
+ */
+template <typename RNGTYPE>
+void SSVSquantreg_impl (rng<RNGTYPE>& stream, double tau, Matrix<>& Y,
+ const Matrix<>& X, unsigned int q,
+ double pi0a0, double pi0b0,
+ unsigned int burnin, unsigned int mcmc,
+ unsigned int thin, unsigned int verbose,
+ Matrix<>& result)
+
+{
+ // define constants
+ unsigned int n_obs = X.rows();
+ unsigned int k = X.cols();
+ double tau_onemintau = tau*(1.0-tau);
+ const unsigned int tot_iter = burnin + mcmc; //total iterations
+ const unsigned int nstore = mcmc / thin; // number of draws to store
+
+ // storage matrices
+ Matrix<bool> gamma(k,1,true,true);
+ Matrix<> gamma_matrix (k, nstore);
+ Matrix<> beta(k,1,true,0.0);
+ Matrix<> beta_matrix(k, nstore);
+
+ COV_TRIAL trial;
+ COV_TRIAL_PREP trial_prep;
+
+ unsigned int n_cov = k;
+
+ // Matrices that will change dimension
+ Matrix<> *X_gamma = new Matrix<>(X);
+ Matrix<> *C = new Matrix<>(n_cov+1,n_cov+1,false);
+ Matrix<> *beta_gamma = new Matrix<>(n_cov,1,false);
+
+ // Matrix with same dimensions
+ Matrix<> U(n_obs,1,false);
+
+ // Gibbs sampler
+ unsigned int count = 0;
+ unsigned int col_index = q;
+
+ // Initial value of pi0
+ double pi0 = stream.rbeta(pi0a0,pi0b0);
+ // Initial values for lambda
+ Matrix<> lambda = stream.rgamma(n_cov-q,1,0.5,0.5);
+ // Initial values for weights
+ Matrix<> weights = stream.rexp(n_obs,1,tau_onemintau);
+ for (unsigned int iter = 0; iter < tot_iter; ++iter) {
+
+trial_prep = QR_SSVS_covariate_trials_prep(*X_gamma, Y, gamma, weights, lambda, tau, n_cov, q);
+C = new Matrix<>(trial_prep.C);
+double logdetminhalf = trial_prep.logdetminhalf;
+U = trial_prep.U;
+
+// updating the indicator variables corresponding to each covariate
+
+ for (unsigned int j=q; j<k; ++j){
+
+ //obtain column index
+ for (unsigned int m=q; m<j; ++m){
+ if (gamma(m) == true){
+ ++col_index;
+ }
+ }
+
+ if (gamma(j) == false){
+ trial = QR_SSVS_covariate_trials_draw_absent(*C,*X_gamma,trial_prep.U,X(_,j),col_index,weights,pi0,lambda(j-q),logdetminhalf,stream);
+ gamma(j) = trial.newtrial;
+ logdetminhalf = trial.logdetminhalf;
+ if (trial.newtrial == true){
+ delete C;
+ C = new Matrix<> (trial.Cnew);
+ delete X_gamma;
+ X_gamma = new Matrix<> (t(selif(t(X),gamma)));
+ ++n_cov;
+ }
+ }
+
+ else{
+ trial = QR_SSVS_covariate_trials_draw_present(*C,col_index,n_obs,pi0,lambda(j-q),logdetminhalf,stream);
+ gamma(j) = trial.newtrial;
+ logdetminhalf = trial.logdetminhalf;
+ if (trial.newtrial == false){
+ delete C;
+ C = new Matrix<> (trial.Cnew);
+ delete X_gamma;
+ X_gamma = new Matrix<> (t(selif(t(X),gamma)));
+ --n_cov;
+ }
+ }
+
+ col_index = q;
+
+ } // end covariate trials loop
+
+if (n_cov == 0){ // Support for null model
+weights = ALaplaceIGaussregress_weights_draw (fabs(Y), stream);
+lambda = stream.rexp(k-q,1,0.5);
+}
+
+else{
+beta_gamma = new Matrix<> (QR_SSVS_beta_draw (*C, stream));
+
+if (q != 0){
+beta(0,0,q-1,0) = (*beta_gamma)(0,0,q-1,0);
+}
+
+for (unsigned int j=q; j<k; ++j){
+
+if (gamma(j) == true){
+beta(j) = (*beta_gamma)(col_index);
+++col_index;
+}
+}
+
+col_index = q;
+
+Matrix<> e = gaxpy(*X_gamma, (-1*(*beta_gamma)), Y);
+Matrix<> abse = fabs(e);
+
+weights = ALaplaceIGaussregress_weights_draw (abse, stream);
+
+lambda = QR_SSVS_lambda_draw(*beta_gamma, gamma, k, q, stream);
+
+delete beta_gamma;
+
+}
+
+delete C;
+
+pi0 = QR_SSVS_pi0_draw(n_cov-q, k-q, pi0a0, pi0b0, stream);
+
+ // store draws in storage matrices
+ if (iter >= burnin && (iter % thin == 0)) {
+ gamma_matrix(_, count) = gamma;
+ beta_matrix(_, count) = beta;
+ ++count;
+ }
+
+ // print output to stdout
+ if(verbose > 0 && iter % verbose == 0) {
+ Rprintf("\n\nSSVSquantreg iteration %i of %i \n",
+ (iter+1), tot_iter);
+ Rprintf("gamma = \n");
+ for (unsigned int r=0; r<k; ++r)
+ Rprintf("%10.5f\n", gamma(r));
+ Rprintf("beta = \n");
+ for (unsigned int r=0; r<k; ++r)
+ Rprintf("%10.5f\n", beta(r));
+ }
+
+ beta = Matrix<>(k,1,true,0.0);
+
+ R_CheckUserInterrupt(); // allow user interrupts
+
+ } // end MCMC loop
+
+ delete X_gamma;
+
+ result = cbind(t(gamma_matrix),t(beta_matrix));
+ } // end SSVSquantreg_impl
+
+extern "C" {
+ void SSVSquantreg(double *sampledata, const int *samplerow,
+ const int *samplecol, const double *tau, const double *Ydata, const int *Yrow,
+ const int *Ycol, const double *Xdata, const int *Xrow,
+ const int *Xcol, const int *q, const int *burnin, const int *mcmc, const int *thin,
+ const int *uselecuyer, const int *seedarray,
+ const int *lecuyerstream, const int *verbose,
+ const double *pi0a0, const double *pi0b0)
+ {
+ // pull together Matrix objects
+ Matrix<> Y(*Yrow, *Ycol, Ydata);
+ Matrix<> X(*Xrow, *Xcol, Xdata);
+ Matrix<> storagematrix;
+ MCMCPACK_PASSRNG2MODEL(SSVSquantreg_impl, *tau, Y, X, *q,
+ *pi0a0, *pi0b0, *burnin, *mcmc, *thin, *verbose,
+ storagematrix);
+
+ const unsigned int size = *samplerow * *samplecol;
+ for (unsigned int h = 0; h < size; ++h){
+ sampledata[h] = storagematrix(h);
+ }
+}
+}
+
+#endif
+
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/r-cran-mcmcpack.git
More information about the debian-science-commits
mailing list