[r-cran-mcmcpack] 43/90: Imported Upstream version 1.0-5
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 c5faba9f7a98f777b37ae7c8e6f08393ac6b65bd
Author: Andreas Tille <tille at debian.org>
Date: Fri Dec 16 08:07:16 2016 +0100
Imported Upstream version 1.0-5
---
DESCRIPTION | 10 +-
HISTORY | 4 +
NAMESPACE | 1 +
R/MCMCirtHier1d.R | 4 +-
R/MCMCirtKdHet.R | 242 ++++++++++++++++++++++
man/MCMCirtKdHet.Rd | 190 ++++++++++++++++++
src/MCMCirtKdHet.cc | 399 ++++++++++++++++++++++++++++++++++++
src/MCMCoprobit.cc | 568 ++++++++++++++++++++++++++--------------------------
8 files changed, 1127 insertions(+), 291 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index 63c905f..cbe1645 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,11 +1,11 @@
Package: MCMCpack
-Version: 1.0-4
-Date: 2009-09-07
+Version: 1.0-5
+Date: 2009-11-28
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>
Maintainer: Andrew D. Martin <admartin at wustl.edu>
-Depends: R (>= 2.8.0), coda (>= 0.11-3), MASS, stats
+Depends: R (>= 2.10.0), coda (>= 0.11-3), MASS, stats
Description: This package contains functions to perform Bayesian
inference using posterior simulation for a number of
statistical models. Most simulation is done in compiled C++
@@ -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-09-08 00:08:21 UTC; adm
+Packaged: 2009-11-28 16:57:19 UTC; adm
Repository: CRAN
-Date/Publication: 2009-09-08 05:04:49
+Date/Publication: 2009-11-29 17:22:58
diff --git a/HISTORY b/HISTORY
index ceedbee..b875b3b 100644
--- a/HISTORY
+++ b/HISTORY
@@ -2,6 +2,10 @@
// Changes and Bug Fixes
//
+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]
+
1.0-3 to 1.0-4
* little change to parameterization of BQR [contributed by Craig Reed]
diff --git a/NAMESPACE b/NAMESPACE
index fb9cbc9..6c4cddf 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -24,6 +24,7 @@ export(
MCMCirt1d,
MCMCirtHier1d,
MCMCirtKd,
+ MCMCirtKdHet,
MCMCirtKdRob,
MCMClogit,
MCMCmetrop1R,
diff --git a/R/MCMCirtHier1d.R b/R/MCMCirtHier1d.R
index 05dfd83..713afa2 100644
--- a/R/MCMCirtHier1d.R
+++ b/R/MCMCirtHier1d.R
@@ -121,7 +121,7 @@
call.=FALSE)
}
- cat("doing betastart ")
+ cat("Generating starting values (glm) for hierarchical parameters:\n")
## starting values are regression of theta.start on Xj
## or passed vector, or same value all beta
if (max(is.na(beta.start))==1) { # beta.start NA
@@ -238,7 +238,7 @@
beta.names <- paste("beta.",beta.names,sep="")
theta.names <- paste("theta.", subject.names, sep = "")
- alpha.beta.names <- paste(rep(c("alpha.","beta."), K),
+ alpha.beta.names <- paste(rep(c("a.","b."), K),
rep(item.names, each = 2),
sep = "")
diff --git a/R/MCMCirtKdHet.R b/R/MCMCirtKdHet.R
new file mode 100644
index 0000000..1c6445a
--- /dev/null
+++ b/R/MCMCirtKdHet.R
@@ -0,0 +1,242 @@
+MCMCirtKdHet <-
+function(datamatrix, dimensions, item.constraints=list(),
+ burnin = 1000, mcmc = 1000, thin=1, verbose = 0, seed = NA,
+ alphabeta.start = NA, b0 = 0, B0=0.04, c0=0, d0=0, store.item = FALSE,
+ store.ability=TRUE,store.sigma=TRUE, drop.constant.items=TRUE) {
+
+ echo.name <- "MCMCirtKdHet"
+
+ # translate variable names to MCMCordfactanal naming convention
+ x <- as.matrix(datamatrix)
+ factors <- dimensions
+ lambda.constraints <- item.constraints
+ lambda.start <- alphabeta.start
+ l0 <- b0
+ L0 <- B0
+ sigma.c0 <- c0
+ sigma.d0 <- d0
+ store.lambda <- store.item
+ store.scores <- store.ability
+ drop.constantvars <- drop.constant.items
+
+ # extract X and variable names from the model formula and frame
+
+ if (is.matrix(x)){
+ if (drop.constantvars==TRUE){
+ x.col.var <- apply(x, 2, var, na.rm=TRUE)
+ keep.inds <- x.col.var>0
+ keep.inds[is.na(keep.inds)] <- FALSE
+ x <- x[,keep.inds]
+ }
+ X <- as.data.frame(x)
+ xvars <- dimnames(X)[[2]]
+ xobs <- dimnames(X)[[1]]
+ N <- nrow(X) # number of observations
+ K <- ncol(X) # number of manifest variables
+ for (i in 1:K){
+ X[,i] <- as.integer(X[,i])
+ if (sum(X[,i] != 0 & X[,i] != 1,na.rm=TRUE) != 0){
+ stop("Data must be 0, 1, and NA only.\n")
+ }
+ X[is.na(X[,i]), i] <- -999
+ }
+ X <- as.matrix(X)
+ }
+ else {
+ stop("Please provide data as a matrix.\n")
+ }
+
+ ## take care of the case where X has no row names
+ if (is.null(xobs)){
+ xobs <- 1:N
+ }
+
+ #check.offset(list(...))
+ check.mcmc.parameters(burnin, mcmc, thin)
+
+ ## setup constraints on Lambda
+ holder <- build.factor.constraints(lambda.constraints, X, K, factors+1)
+ Lambda.eq.constraints <- holder[[1]]
+ Lambda.ineq.constraints <- holder[[2]]
+ X.names <- holder[[3]]
+
+ ## setup prior on Lambda
+ holder <- form.factload.norm.prior(l0, L0, K, factors+1, X.names)
+ Lambda.prior.mean <- holder[[1]]
+ Lambda.prior.prec <- holder[[2]]
+ Lambda.prior.mean[,1] <- Lambda.prior.mean[,1] * -1
+
+
+
+ # seeds
+ seeds <- form.seeds(seed)
+ lecuyer <- seeds[[1]]
+ seed.array <- seeds[[2]]
+ lecuyer.stream <- seeds[[3]]
+
+
+ ## Starting values for Lambda
+ Lambda <- matrix(0, K, factors+1)
+ if (is.na(lambda.start)){# sets Lambda to equality constraints & 0s
+ for (i in 1:K){
+ for (j in 1:(factors+1)){
+ if (Lambda.eq.constraints[i,j]==-999){
+ if(Lambda.ineq.constraints[i,j]==0){
+ if (j==1){
+ probit.out <- glm(as.factor(X[X[,i]!=-999,i])~1,
+ family=binomial(link=probit))
+ probit.beta <- coef(probit.out)
+ Lambda[i,j] <- probit.beta[1]
+ }
+ }
+ if(Lambda.ineq.constraints[i,j]>0){
+ Lambda[i,j] <- 1.0
+ }
+ if(Lambda.ineq.constraints[i,j]<0){
+ Lambda[i,j] <- -1.0
+ }
+ }
+ else Lambda[i,j] <- Lambda.eq.constraints[i,j]
+ }
+ }
+ }
+ else if (is.matrix(lambda.start)){
+ if (nrow(lambda.start)==K && ncol(lambda.start)==(factors+1))
+ Lambda <- lambda.start
+ else {
+ cat("Starting values not of correct size for model specification.\n")
+ stop("Please respecify and call ", echo.name, "() again\n")
+ }
+ }
+ else if (length(lambda.start)==1 && is.numeric(lambda.start)){
+ Lambda <- matrix(lambda.start, K, factors+1)
+ for (i in 1:K){
+ for (j in 1:(factors+1)){
+ if (Lambda.eq.constraints[i,j] != -999)
+ Lambda[i,j] <- Lambda.eq.constraints[i,j]
+ }
+ }
+ }
+ else {
+ cat("Starting values neither NA, matrix, nor scalar.\n")
+ stop("Please respecify and call ", echo.name, "() again\n")
+ }
+
+
+ ## define holder for posterior sample
+ if (store.scores == FALSE && store.lambda == FALSE && store.sigma == FALSE){
+ stop("Please specify parameters to be stored.\n")
+ }
+ else if (store.scores == TRUE && store.lambda == FALSE && store.sigma == FALSE){
+ sample <- matrix(data=0, mcmc/thin, (factors+1)*N)
+ }
+ else if(store.scores == FALSE && store.lambda == TRUE && store.sigma == FALSE) {
+ sample <- matrix(data=0, mcmc/thin, K*(factors+1))
+ }
+ else if(store.scores==TRUE && store.lambda==TRUE && store.sigma == FALSE) {
+ sample <- matrix(data=0, mcmc/thin, K*(factors+1)+(factors+1)*N)
+ }
+ else if (store.scores == FALSE && store.lambda == FALSE && store.sigma == TRUE){
+ sample <- matrix(data=0, mcmc/thin, N)
+ }
+ else if (store.scores == TRUE && store.lambda == FALSE && store.sigma == TRUE){
+ sample <- matrix(data=0, mcmc/thin, (factors+1)*N + N)
+ }
+ else if(store.scores == FALSE && store.lambda == TRUE && store.sigma == TRUE) {
+ sample <- matrix(data=0, mcmc/thin, K*(factors+1)+N)
+ }
+ else if(store.scores==TRUE && store.lambda==TRUE && store.sigma == TRUE) {
+ sample <- matrix(data=0, mcmc/thin, K*(factors+1)+(factors+1)*N + N)
+ }
+
+
+ ## create templates
+ posterior <- .C("irtKdHetpost",
+ samdata = as.double(sample),
+ samrow = as.integer(nrow(sample)),
+ samcol = as.integer(ncol(sample)),
+ Xdata = as.integer(X),
+ Xrow = as.integer(nrow(X)),
+ Xcol = as.integer(ncol(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),
+ Lamstartdata = as.double(Lambda),
+ Lamstartrow = as.integer(nrow(Lambda)),
+ Lamstartcol = as.integer(ncol(Lambda)),
+ Lameqdata = as.double(Lambda.eq.constraints),
+ Lameqrow = as.integer(nrow(Lambda.eq.constraints)),
+ Lameqcol = as.integer(ncol(Lambda.ineq.constraints)),
+ Lamineqdata = as.double(Lambda.ineq.constraints),
+ Lamineqrow = as.integer(nrow(Lambda.ineq.constraints)),
+ Lamineqcol = as.integer(ncol(Lambda.ineq.constraints)),
+ Lampmeandata = as.double(Lambda.prior.mean),
+ Lampmeanrow = as.integer(nrow(Lambda.prior.mean)),
+ Lampmeancol = as.integer(ncol(Lambda.prior.prec)),
+ Lampprecdata = as.double(Lambda.prior.prec),
+ Lampprecrow = as.integer(nrow(Lambda.prior.prec)),
+ Lamppreccol = as.integer(ncol(Lambda.prior.prec)),
+ storelambda = as.integer(store.lambda),
+ storescores = as.integer(store.scores),
+ storesigma = as.integer(store.sigma),
+ sigmapriorc = as.double(sigma.c0),
+ sigmapriord = as.double(sigma.d0),
+ package="MCMCpack"
+ )
+
+ output <- mcmc(data=matrix(posterior$samdata, posterior$samrow, posterior$samcol,byrow=FALSE),start=1, end=mcmc, thin=thin)
+
+ par.names <- NULL
+ if (store.lambda==TRUE){
+
+ alpha.hold <- paste("alpha", X.names, sep=".")
+ beta.hold <- paste("beta", X.names, sep = ".")
+ beta.hold <- rep(beta.hold, factors, each=factors)
+ beta.hold <- paste(beta.hold, 1:factors, sep=".")
+
+ Lambda.names <- t(cbind(matrix(alpha.hold, K, 1),
+ matrix(beta.hold,K,factors,byrow=TRUE)))
+ dim(Lambda.names) <- NULL
+
+ par.names <- c(par.names, Lambda.names)
+ }
+
+ if (store.scores==TRUE){
+ phi.names <- paste(paste("theta",
+ rep(xobs, each=(factors+1)), sep="."),
+ rep(0:factors,(factors+1)), sep=".")
+ par.names <- c(par.names, phi.names)
+
+ }
+
+ if (store.sigma==TRUE){
+ sigma.names <- paste("sigma", rep(xobs), sep=".")
+ par.names <- c(par.names, sigma.names)
+ }
+
+ varnames(output) <- par.names
+
+ # get rid of columns for constrained parameters
+ output.df <- as.data.frame(as.matrix(output))
+
+ output.var <- sd(output.df)
+ output.df <- output.df[,output.var != 0]
+ output <- mcmc(as.matrix(output.df), start=burnin+1, end=burnin+mcmc,
+ thin=thin)
+
+ # add constraint info so this isn't lost
+ attr(output, "constraints") <- lambda.constraints
+ attr(output, "n.manifest") <- K
+ attr(output, "n.factors") <- factors
+ attr(output, "title") <-
+ "MCMCpack Heteroskedastic K-Dimensional Item Response Theory Model Posterior Sample"
+
+ return(output)
+
+
+ }
+
diff --git a/man/MCMCirtKdHet.Rd b/man/MCMCirtKdHet.Rd
new file mode 100644
index 0000000..8acf1d6
--- /dev/null
+++ b/man/MCMCirtKdHet.Rd
@@ -0,0 +1,190 @@
+\name{MCMCirtKdHet}
+\Rdversion{1.1}
+\alias{MCMCirtKdHet}
+\title{
+Markov Chain Monte Carlo for Heteroskedastic K-Dimensional Item Response Theory Model}
+\description{
+This function generates a sample from the posterior distribution of a
+ heteroskedastic K-dimensional item response theory (IRT) model, with standard
+ normal priors on the subject abilities (ideal points),
+ normal priors on the item parameters, and inverse-gamma priors on
+ subject error variances. To maintain identification and comparability
+ with results of the homoskedastic estimator, the mean root subject
+ error precision is constrained to one. The user supplies data and
+ priors, and a sample from the posterior distribution is returned as an mcmc
+ object, which can be subsequently analyzed with functions provided in the coda
+ package.}
+\usage{
+MCMCirtKdHet(datamatrix, dimensions, item.constraints = list(), burnin = 1000,
+mcmc = 1000, thin = 1, verbose = 0, seed = NA, alphabeta.start = NA, b0 = 0,
+B0 = 0.04, c0 = 0, d0 = 0, store.item = FALSE, store.ability = TRUE,
+store.sigma = TRUE, drop.constant.items = TRUE)
+}
+\arguments{
+ \item{datamatrix}{The matrix of data. Must be 0, 1, or NA.
+ It is of dimensionality subjects by items.}
+
+ \item{dimensions}{The number of dimensions in the latent space.}
+
+ \item{item.constraints}{List of lists specifying possible equality
+ or simple inequality constraints on the item parameters. A typical
+ entry in the list has one of three forms: \code{rowname=list(d,c)}
+ which will constrain the dth item parameter for the item named
+ rowname to be equal to c, \code{rowname=list(d,"+")} which will
+ constrain the dth item parameter for the item named rowname to be
+ positive, and\code{rowname=list(d, "-")} which will constrain the dth
+ item parameter for the item named rowname to be negative. If x is a
+ matrix without row names defaults names of ``V1", ``V2", ... , etc
+ will be used. In a K dimensional model, the first item parameter for
+ item \eqn{i}{i} is the difficulty parameter (\eqn{\alpha_i}{alpha_i}),
+ the second item parameter is the discrimation parameter on dimension
+ 1 (\eqn{\beta_{i,1}}{beta_{i,1}}), the third item parameter is the
+ discrimation parameter on dimension 2
+ (\eqn{\beta_{i,2}}{beta_{i,2}}), ..., and the (K+1)th item parameter
+ is the discrimation parameter on dimension K
+ (\eqn{\beta_{i,1}}{beta_{i,1}}).
+ The item difficulty parameters (\eqn{\alpha}{alpha}) should
+ generally not be constrained.
+ }
+
+ \item{burnin}{The number of burn-in iterations for the sampler.}
+
+ \item{mcmc}{The number of iterations for the sampler.}
+
+ \item{thin}{The thinning interval used in the simulation. The number of
+ 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 then every \code{verbose}th iteration will be printed to the
+ screen.}
+
+ \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
+ 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{alphabeta.start}{The starting values for the
+ \eqn{\alpha}{alpha} and \eqn{\beta}{beta} difficulty and
+ discrimination parameters. If \code{alphabeta.start} is set to a
+ scalar the starting value for all unconstrained item parameters will
+ be set to that scalar. If \code{alphabeta.start} is a matrix of
+ dimension \eqn{(K+1) \times items}{(K+1) x items} then the
+ \code{alphabeta.start} matrix is used as the starting values (except
+ for equality-constrained elements). If \code{alphabeta.start} is set
+ to \code{NA} (the default) then starting values for unconstrained
+ elements are set to values generated from a series of proportional
+ odds logistic regression fits, and starting values for inequality
+ constrained elements are set to either 1.0 or -1.0 depending on the
+ nature of the constraints. }
+
+ \item{b0}{The prior means of the
+ \eqn{\alpha}{alpha} and \eqn{\beta}{beta} difficulty and
+ discrimination parameters, stacked for all items.
+ If a scalar is passed, it
+ is used as the prior mean for all items.}
+
+ \item{B0}{The prior precisions (inverse variances) of the
+ independent normal prior on the item parameters.
+ Can be either a scalar or a matrix of dimension
+ \eqn{(K+1) \times items}{(K+1) x items}.}
+
+ \item{c0}{The first parameter of the inverse gamma prior on the
+ subject-specific variance parameters. This can be thought of as
+ the number of bills that the prior information is equivalent to.
+ This scalar value is common across all subjects (legislators) and
+ defaults to an uninformative prior. NOTE: regardless of the value
+ provided, identification is provided by a constraint on the mean
+ root subject specific variance.}
+
+ \item{d0}{The second parameter of the inverse gamma prior on the
+ subject-specific variance parameters. This can be thought of as
+ the sum of square error that the prior information is equivalent to.
+ This scalar value is common across all subjects (legislators) and
+ defaults to an uninformative prior. NOTE: regardless of the value
+ provided, identification is provided by a constraint on the mean
+ root subject specific variance.}
+
+ \item{store.item}{A switch that determines whether or not to
+ store the item parameters for posterior analysis.
+ \emph{NOTE: In applications with many items
+ this takes an enormous amount of memory. If you have many items
+ and want to want to store the item parameters you may want to thin
+ the chain heavily}. By default, the
+ item parameters are not stored.}
+
+ \item{store.ability}{A switch that determines whether or not to store
+ the subject abilities for posterior analysis. \emph{NOTE: In
+ applications with many subjects this takes an enormous amount of
+ memory. If you have many subjects and want to want to store the ability
+ parameters you may want to thin the chain heavily}. By default, the
+ ability parameters are all stored.}
+
+ \item{store.sigma}{A switch that determines whether or not to store
+ the subject-specific variances for posterior analysis. \emph{NOTE: In
+ applications with many subjects this takes an enormous amount of
+ memory. If you have many subjects and want to want to store the ability
+ parameters you may want to thin the chain heavily}. By default, the
+ subject-specific variance parameters are all stored.}
+
+ \item{drop.constant.items}{A switch that determines whether or not
+ items that have no variation
+ should be deleted before fitting the model. Default = TRUE.}
+}
+\value{
+ An mcmc object that contains the posterior sample. This
+ object can be summarized by functions provided by the coda package.
+}
+\references{
+ Benjamin E. Lauderdale. ``Unpredictable Voters in Ideal Point Estimation''
+\emph{Political Analysis.} forthcoming.
+}
+\author{
+Benjamin E. Lauderdale, \email{blauderd at princeton.edu},
+\url{http://www.princeton.edu/~blauderd/}.
+
+Modified from \code{\link[MCMCpack]{MCMCirtKd}} and \code{\link[MCMCpack]{MCMCordfactanal}}. Suggestions for additional options are welcome.
+}
+
+\seealso{\code{\link[coda]{plot.mcmc}},\code{\link[coda]{summary.mcmc}},
+\code{\link[MCMCpack]{MCMCirtKd}}
+}
+
+\examples{
+
+data(Senate)
+Y <- as.matrix(Senate[,6:677])
+
+Hompost <- MCMCirtKd(Y,1,b0=0,B0=0.04,burn=1000,mcmc=1000,thin=1,verbose=250)
+Hetpost <- MCMCirtKdHet(Y,1,b0=0,B0=0.04,burn=1000,mcmc=1000,thin=1,verbose=250)
+
+SenatorNames <- Senate[,5]
+HomoskedasticIdealPointEstimates <- colMeans(Hompost)[1:102]
+HeteroskedasticIdealPointEstimates <- colMeans(Hetpost)[1:102]
+HeteroskedasticSigmaEstimates <- colMeans(Hetpost)[103:204]
+
+plot(HomoskedasticIdealPointEstimates, HeteroskedasticIdealPointEstimates,
+cex= HeteroskedasticSigmaEstimates,xlab="Ideal Points (Homoskedastic)",
+ylab="Ideal Points (Heteroskedastic)",
+main="Comparison of Ideal Point Estimates for the 106th Senate",
+xlim=c(-2.5,2.5),ylim=c(-2.5,2.5))
+for (i in 1:102){
+ if (rank(-HeteroskedasticSigmaEstimates)[i] <= 10){
+ text(HomoskedasticIdealPointEstimates[i],
+ HeteroskedasticIdealPointEstimates[i],SenatorNames[i],
+ pos=3-sign(HomoskedasticIdealPointEstimates[i]),cex=0.75)
+ }
+}
+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/src/MCMCirtKdHet.cc b/src/MCMCirtKdHet.cc
new file mode 100644
index 0000000..dc3dca6
--- /dev/null
+++ b/src/MCMCirtKdHet.cc
@@ -0,0 +1,399 @@
+//////////////////////////////////////////////////////////////////////////
+// MCMCordfactanal.cc is C++ code to estimate an ordinal data
+// factor analysis model
+//
+// 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.
+//
+// revised version of older MCMCordfactanal
+// 7/16/2004 KQ
+// updated to new version of Scythe ADM 7/24/2004
+// fixed a bug pointed out by Alexander Raach 1/16/2005 KQ
+//
+// Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+// Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+// and Jong Hee Park
+//////////////////////////////////////////////////////////////////////////
+
+
+#ifndef MCMCIRTKDHET_CC
+#define MCMCIRTKDHET_CC
+
+#include <iostream>
+
+#include "matrix.h"
+#include "algorithm.h"
+#include "distributions.h"
+#include "stat.h"
+#include "la.h"
+#include "ide.h"
+#include "smath.h"
+#include "MCMCrng.h"
+#include "MCMCfcds.h"
+
+#include <R.h> // needed to use Rprintf()
+#include <R_ext/Utils.h> // needed to allow user interrupts
+
+typedef Matrix<double,Row,View> rmview;
+
+using namespace std;
+using namespace scythe;
+
+template <typename RNGTYPE>
+void MCMCirtKdHet_impl(rng<RNGTYPE>& stream,
+ const Matrix<int>& X, Matrix<>& Lambda,
+ const Matrix<>& Lambda_eq,
+ const Matrix<>& Lambda_ineq,
+ const Matrix<>& Lambda_prior_mean,
+ const Matrix<>& Lambda_prior_prec,
+ const double sigmapriorc,
+ const double sigmapriord,
+ bool storelambda, bool storescores,
+ bool storesigma, unsigned int burnin,
+ unsigned int mcmc, unsigned int thin,
+ unsigned int verbose,
+ Matrix<>& output) {
+
+ // constants
+ const unsigned int K = X.cols(); // number of manifest variables
+ const unsigned int N = X.rows(); // number of observations
+ const unsigned int D = Lambda.cols(); // # of factors (incl constant)
+ const unsigned int tot_iter = burnin + mcmc;
+ const unsigned int nsamp = mcmc / thin;
+ const Matrix<> I = eye<double>(D-1);
+ const Matrix<bool> Lambda_free_indic(K, D);
+ for (unsigned int i = 0; i < (K * D); ++i)
+ if (Lambda_eq(i) == -999)
+ Lambda_free_indic(i) = true;
+
+ const Matrix<> Psi = eye<double>(K);
+ const Matrix<> Psi_inv = eye<double>(K);
+ const double c0 = sigmapriorc;
+ const double d0 = sigmapriord;
+
+ //Rprintf("Switches are %i %i %i\n", storelambda, storescores, storesigma);
+
+ // starting values for phi, Xstar, and sigma
+ Matrix<> phi(N, D-1);
+ phi = cbind(ones<double>(N,1), phi);
+ Matrix<> sigma2 = ones<double>(N,1);
+ Matrix<> sigma = ones<double>(N,1);
+ Matrix<> Xstar(N, K);
+ double sigma_norm = 1;
+
+ // storage matrices (row major order)
+ Matrix<> Lambda_store;
+ if (storelambda){
+ Lambda_store = Matrix<>(nsamp, K*D);
+ }
+ Matrix<> phi_store;
+ if (storescores){
+ phi_store = Matrix<>(nsamp, N*D);
+ }
+ Matrix<> sigma_store;
+ if (storesigma){
+ sigma_store = Matrix<>(nsamp, N*1);
+ }
+
+ ///////////////////
+ // Gibbs Sampler //
+ ///////////////////
+ int count = 0;
+ for (unsigned int iter = 0; iter < tot_iter; ++iter) {
+
+ // sample Xstar
+ for (unsigned int i = 0; i < N; ++i) {
+ Matrix<> X_mean = Lambda * t(phi(i,_));
+ for (unsigned int j = 0; j < K; ++j) {
+ if (X(i,j) == -999) // if missing
+ Xstar(i,j) = stream.rnorm(X_mean[j], sigma(i,0));
+ if (X(i,j) == 0) // if not missing
+ Xstar(i,j) = stream.rtnorm_combo(X_mean[j], sigma2(i,0), -300, 0);
+ if (X(i,j) == 1) // if not missing
+ Xstar(i,j) = stream.rtnorm_combo(X_mean[j], sigma2(i,0), 0, 300);
+ }
+ }
+
+ // sample phi
+ Matrix<> Lambda_const = Lambda(_,0);
+ Matrix<> Lambda_rest = Lambda(0, 1, K-1, D-1);
+
+ for (unsigned int i=0; i<N; ++i){
+ Matrix<> phi_post_var = invpd(I + crossprod(Lambda_rest) / sigma2(i,0) );
+ Matrix<> phi_post_C = cholesky(phi_post_var);
+ Matrix<> phi_post_mean = phi_post_var *
+ (t(Lambda_rest) * (t(Xstar(i,_))-Lambda_const)) / sigma2(i,0);
+ Matrix<> phi_samp = gaxpy(phi_post_C, stream.rnorm(D-1, 1, 0, 1),
+ phi_post_mean);
+ for (unsigned int j=0; j<(D-1); ++j)
+ phi(i,j+1) = phi_samp[j];
+ }
+
+
+ // sample Lambda
+
+ for (unsigned int i=0; i<K; ++i){
+ const Matrix<bool> free_indic = t(Lambda_free_indic(i,_));
+ const Matrix<bool> not_free_indic = 1 - free_indic;
+ if (sumc(free_indic)[0] > 0 &&
+ sumc(not_free_indic)[0] > 0){ // both constrnd & unconstrnd
+ const Matrix<> phifree_i = t(selif(t(phi), free_indic));
+ const Matrix<> mulamfree_i = selif(t(Lambda_prior_mean(i,_)),
+ free_indic); // prior mean
+ const Matrix<> hold = selif(t(Lambda_prior_prec(i,_)),
+ free_indic);
+ Matrix<> sig2lamfree_inv_i =
+ eye<double>(hold.rows()); // prior prec
+ for (unsigned int j=0; j<(hold.rows()); ++j)
+ sig2lamfree_inv_i(j,j) = hold[j];
+ const Matrix<> Lambdacon_i =
+ selif(t(Lambda(i,_)), not_free_indic);
+ const Matrix<> phicon_i = t(selif(t(phi), not_free_indic));
+ const Matrix<> X_i = Xstar(_,i);
+ const Matrix<> newX_i = gaxpy((-1.0*phicon_i), Lambdacon_i, X_i);
+
+ for (unsigned int R=0; R<N; ++R){ // multiply through by inverse sigma
+ for (int Q=0; Q<D; ++ Q){
+ phifree_i(R,Q) = phifree_i(R,Q) / sigma(R,0);
+ }
+ newX_i(R,0) = newX_i(R,0) /sigma(R,0);
+ }
+
+ const Matrix<> Lam_post_var = invpd(sig2lamfree_inv_i +
+ Psi_inv(i,i) * crossprod(phifree_i));
+ const Matrix<> Lam_post_C = cholesky(Lam_post_var);
+ const Matrix<> Lam_post_mean = Lam_post_var *
+ (sig2lamfree_inv_i * mulamfree_i + Psi_inv(i,i) *
+ t(phifree_i) * newX_i);
+
+ Matrix<double> Lambdafree_i =
+ gaxpy(Lam_post_C, stream.rnorm(hold.rows(), 1, 0, 1), Lam_post_mean);
+
+ // check to see if inequality constraints hold
+ const Matrix<> Lambda_ineq_vec = Lambda_ineq(i,_);
+ double ineq_holds = 0;
+ int Lam_count = 0;
+ for (unsigned int j=0; j<D; ++j){
+ if (free_indic[j]==1)
+ ineq_holds = std::min(ineq_holds,
+ Lambda_ineq_vec[j] *
+ Lambdafree_i[Lam_count]);
+ ++Lam_count;
+ }
+ while (ineq_holds < 0){
+ Lambdafree_i =
+ gaxpy(Lam_post_C, stream.rnorm(hold.rows(), 1, 0, 1), Lam_post_mean);
+ Lam_count = 0;
+ double test = 0;
+ for (int j=0; j<D; ++j){
+ if (free_indic[j]==1){
+ Matrix<> prodcheck =
+ Lambda_ineq_vec[j]*Lambdafree_i[Lam_count];
+ test = std::min(test, prodcheck[0]);
+ ++Lam_count;
+ }
+ }
+ ineq_holds = test;
+ }
+
+ // put draw into Lambda
+ Lam_count = 0;
+ for (unsigned int j=0; j<D; ++j){
+ if (free_indic[j] == 1){
+ Lambda(i,j) = Lambdafree_i[Lam_count];
+ ++Lam_count;
+ }
+ }
+ }
+ else if (sumc(free_indic)[0] > 0){ // just unconstrained
+ const Matrix<> phifree_i = t(selif(t(phi), free_indic));
+ const Matrix<> mulamfree_i = selif(t(Lambda_prior_mean(i,_)),
+ free_indic); // prior mean
+ const Matrix<> hold = selif(t(Lambda_prior_prec(i,_)),
+ free_indic);
+ Matrix<> sig2lamfree_inv_i =
+ eye<double>(hold.rows()); // prior prec
+ for (unsigned int j=0; j<hold.rows(); ++j)
+ sig2lamfree_inv_i(j,j) = hold[j];
+
+ Matrix<> sX_i = ones<double>(N,1);
+ for (unsigned int R=0; R<N; ++R){ // multiply through by inverse sigma
+ for (int Q=0; Q<D; ++ Q){
+ phifree_i(R,Q) = phifree_i(R,Q) / sigma(R,0);
+ }
+ sX_i(R,0) = Xstar(R,i) / sigma(R,0);
+ }
+
+ const Matrix<> Lam_post_var = invpd(sig2lamfree_inv_i +
+ Psi_inv(i,i) * crossprod(phifree_i)); //
+
+ const Matrix<> Lam_post_C = cholesky(Lam_post_var);
+ const Matrix<> Lam_post_mean = Lam_post_var *
+ (sig2lamfree_inv_i * mulamfree_i + Psi_inv(i,i) *
+ t(phifree_i) * sX_i); //
+ Matrix<> Lambdafree_i =
+ gaxpy(Lam_post_C, stream.rnorm(hold.rows(), 1, 0, 1), Lam_post_mean);
+
+ // check to see if inequality constraints hold
+ Matrix<> Lambda_ineq_vec = Lambda_ineq(i,_);
+ double ineq_holds = 0;
+ for (unsigned int j=0; j<D; ++j){
+ ineq_holds =
+ std::min(ineq_holds, Lambda_ineq_vec[j]*Lambdafree_i[j]);
+ }
+ while (ineq_holds < 0){
+ Lambdafree_i =
+ gaxpy(Lam_post_C, stream.rnorm(hold.rows(), 1, 0, 1), Lam_post_mean);
+ double test = 0;
+ for (unsigned int j=0; j<D; ++j){
+ //if (free_indic[j]==1)
+ double prodcheck = Lambda_ineq_vec[j]*Lambdafree_i[j];
+ test = std::min(test, prodcheck);
+ }
+ ineq_holds = test;
+ }
+
+ // put draw into Lambda
+ for (unsigned int j=0; j<D; ++j){
+ Lambda(i,j) = Lambdafree_i[j];
+ }
+ }
+ }
+
+
+ // sample sigma
+ for (unsigned int i=0; i<N; ++i){
+ Matrix<> phi_i = phi(i,_);
+ Matrix<> Xstar_i = Xstar(i,_);
+ const Matrix<> e = gaxpy(phi_i, (-1*t(Lambda)), Xstar_i);
+ const Matrix<> SSE = crossprod (t(e));
+ const double c_post = (c0 + K) * 0.5;
+ const double d_post = (d0 + SSE[0]) * 0.5;
+ sigma2(i,0) = stream.rigamma(c_post, d_post);
+ }
+
+ Matrix<> sigma = sqrt(sigma2); // store as sigma
+ sigma_norm = N / sum(pow(sigma,-1)); // renormalize inverse sigma to sum to N
+ for (unsigned int i=0; i<N; ++i){
+ sigma(i,0) = sigma(i,0) / sigma_norm;
+ sigma2(i,0) = pow(sigma(i,0),2);
+ }
+
+
+ // print results to screen
+ if (verbose > 0 && iter % verbose == 0){
+
+ //Rprintf("xstar 22, xstar 23, xstar 24, phi, sigma = \n");
+ //for (int i=0; i<N; ++i){
+ //Rprintf("%10.3f",Xstar(i,21));
+ //Rprintf("%10.3f",Xstar(i,22));
+ //Rprintf("%10.3f",Xstar(i,23));
+ //Rprintf("%10.3f",phi(i,1));
+ //Rprintf("%10.3f",sigma(i,0));
+ //Rprintf("\n");
+ //}
+
+ //Rprintf("Lambda = \n");
+ //for (unsigned int i=0; i<K; ++i){
+ //for (unsigned int j=0; j<D; ++j){
+ //Rprintf("%10.5f", Lambda(i,j));
+ //}
+ //Rprintf("\n");
+ //}
+
+ //Rprintf("\nsigma renormalization factor");
+ //Rprintf("%10.3f \n",sigma_norm);
+
+ Rprintf("\n\nMCMCirtKdHet iteration %i of %i \n", (iter+1), tot_iter);
+
+ }
+
+ // store results
+ if ((iter >= burnin) && ((iter % thin==0))) {
+
+ // store Lambda
+ if (storelambda==1)
+ rmview(Lambda_store(count, _)) = Lambda;
+
+ // store phi
+ if (storescores==1)
+ rmview(phi_store(count, _)) = phi;
+
+ // store sigma
+ if (storesigma==1)
+ rmview(sigma_store(count, _)) = sigma;
+
+ count++;
+
+ }
+
+ // allow user interrupts
+ R_CheckUserInterrupt();
+
+ } // end MCMC loop
+
+ if (storelambda == 1){
+ output = Lambda_store;
+ if (storescores == 1)
+ output = cbind(output, phi_store);
+ if (storesigma == 1)
+ output = cbind(output, sigma_store);
+ }
+ if (storelambda == 0) {
+ if (storescores == 1){
+ output = phi_store;
+ if (storesigma == 1)
+ output = cbind(output, sigma_store);
+ }
+ if (storescores == 0) {
+ output = sigma_store;
+ }
+ }
+
+}
+
+
+extern "C"{
+
+ // function called by R to fit model
+ void irtKdHetpost(double *samdata, const int *samrow, const int *samcol, const int *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 *Lamstartdata, const int *Lamstartrow, const int *Lamstartcol, const double *Lameqdata, const int *Lameqrow, const int *Lameqcol, const double *Lamineqdata, const int *Lamineqrow, const int *Lamineqcol, const [...]
+
+ // put together matrices
+ const Matrix<int> X(*Xrow, *Xcol, Xdata);
+ Matrix<> Lambda(*Lamstartrow, *Lamstartcol, Lamstartdata);
+ const Matrix<> Lambda_eq(*Lameqrow, *Lameqcol, Lameqdata);
+ const Matrix<> Lambda_ineq(*Lamineqrow, *Lamineqcol, Lamineqdata);
+ const Matrix<> Lambda_prior_mean(*Lampmeanrow, *Lampmeancol,
+ Lampmeandata);
+ const Matrix<> Lambda_prior_prec(*Lampprecrow, *Lamppreccol,
+ Lampprecdata);
+
+
+
+ // return output
+ Matrix<double> output;
+ MCMCPACK_PASSRNG2MODEL(MCMCirtKdHet_impl, X, Lambda, Lambda_eq, Lambda_ineq, Lambda_prior_mean,
+ Lambda_prior_prec, *sigmapriorc, *sigmapriord, *storelambda,
+ *storescores, *storesigma,
+ *burnin, *mcmc, *thin, *verbose, output);
+
+
+ for (unsigned int i = 0; i < output.size(); ++i)
+ samdata[i] = output(i);
+
+ }
+
+}
+
+#endif
diff --git a/src/MCMCoprobit.cc b/src/MCMCoprobit.cc
index 7cf4192..aaef110 100644
--- a/src/MCMCoprobit.cc
+++ b/src/MCMCoprobit.cc
@@ -50,146 +50,145 @@ using namespace std;
using namespace scythe;
static inline double lnmulttdens(const Matrix<>& theta,
- const Matrix<>& mu,
- const Matrix<>& C,
- const double df){
-
- const int d = theta.size();
- //const Matrix<> z = t(theta - mu) * C;
- // C is now C' if VC mat is C C'
- const Matrix<> z = C * (theta - mu);
- double zsumsq = 0;
- for (int i=0; i<d; ++i){
- zsumsq += std::pow(z[i], 2);
- }
-
- return ( (-(df + d)/2) * std::log(1 + (1/df) * zsumsq) );
+ const Matrix<>& mu,
+ const Matrix<>& C,
+ const double df){
+
+ const int d = theta.size();
+ //const Matrix<> z = t(theta - mu) * C;
+ // C is now C' if VC mat is C C'
+ const Matrix<> z = C * (theta - mu);
+ double zsumsq = 0;
+ for (int i=0; i<d; ++i){
+ zsumsq += std::pow(z[i], 2);
+ }
+
+ return ( (-(df + d)/2) * std::log(1 + (1/df) * zsumsq) );
}
// function that transforms alpha to gamma
Matrix<> gamma2alpha(const Matrix<>& gamma){
- const int m = gamma.rows() - 2;
- Matrix<> alpha(m, 1);
- alpha[0] = std::log(gamma[1]);
- for (int j=1; j< m ; ++j){
- alpha[j] = std::log(gamma[j+1] - gamma[j]);
- }
- return alpha;
+ const int m = gamma.rows() - 2;
+ Matrix<> alpha(m, 1);
+ alpha[0] = std::log(gamma[1]);
+ for (int j=1; j< m ; ++j){
+ alpha[j] = std::log(gamma[j+1] - gamma[j]);
+ }
+ return alpha;
}
// function that transforms alpha to gamma
Matrix<> alpha2gamma(const Matrix<>& alpha){
- const int m = alpha.rows();
- Matrix<> gamma(m+2, 1);
- gamma[0] = -300;
- gamma[m+1] = 300;
- for (int j=1; j<m+1 ; ++j){
- double gamma_sum = 0.0;
- for(int i=0;i<j; ++i){
- gamma_sum +=exp(alpha[i]);
- }
- gamma[j] = gamma_sum;
- }
- return gamma;
+ const int m = alpha.rows();
+ Matrix<> gamma(m+2, 1);
+ gamma[0] = -300;
+ gamma[m+1] = 300;
+ for (int j=1; j<m+1 ; ++j){
+ double gamma_sum = 0.0;
+ for(int i=0;i<j; ++i){
+ gamma_sum +=exp(alpha[i]);
+ }
+ gamma[j] = gamma_sum;
+ }
+ return gamma;
}
double lndmvn_jhp(const Matrix<double> &x, const Matrix<double> &mu,
- const Matrix<double> &Sigma)
+ const Matrix<double> &Sigma)
{
- int k = Sigma.cols();
- double first = (-k/2.0) * ::log(2 * M_PI) -0.5 * ::log(det(Sigma));
- Matrix< > second = ::t(x-mu)*invpd(Sigma)*(x-mu);
- return (first - 0.5*second[0]);
+ int k = Sigma.cols();
+ double first = (-k/2.0) * ::log(2 * M_PI) -0.5 * ::log(det(Sigma));
+ Matrix< > second = ::t(x-mu)*invpd(Sigma)*(x-mu);
+ return (first - 0.5*second[0]);
}
// orpobit_logpost
static inline double oprobit_logpost(const Matrix<double>& nY, const Matrix<double>& X,
- const Matrix<double>& alpha,
- const Matrix<double>& alpha_prior_mean,
- const Matrix<double>& alpha_prior_var,
- const Matrix<double>& beta){
-
- // likelihood
- double loglike = 0.0;
- const int n = nY.rows();
- const int ncat = alpha.rows() + 1;
-
- // the linear predictor
- Matrix<> mu = X*beta;
- Matrix<> gamma = alpha2gamma(alpha);
-
- // compute prob
- Matrix<> cat_prob(n, ncat-1);
- //cat_prob: CATegory specific PROBability
- //the first col of cat_prob = pnorm(gamma_1 - mu)
- //thus, the col number is ncat - 1
- Matrix<> prob(n, ncat);
- for (int j=0; j< ncat-1; ++j){
- for (int i=0; i<n ; ++i){
- cat_prob(i, j) = pnorm(gamma[j+1] - mu[i], 0.0, 1.0);
- }
- }
- prob(_, ncat-1) = 1 - cat_prob(_, ncat-2); // top category
- prob(_, 0) = cat_prob(_, 0); // bottom category
- for (int j=1; j<(ncat-1); ++j){ // middle categories are actually differences of cumulatives
- prob(_, j) = cat_prob(_,j) - cat_prob(_, j-1);
- }
-
- // the loglikelihood
-
- for (int i=0; i<n; ++i){
- int ind = (int) nY[i];
- loglike += std::log(prob(i,ind-1));
- }
-
- // prior
- double logprior = 0.0;
- logprior = lndmvn_jhp(alpha, alpha_prior_mean, alpha_prior_var);
-
- return (loglike + logprior);
-
+ const Matrix<double>& alpha,
+ const Matrix<double>& alpha_prior_mean,
+ const Matrix<double>& alpha_prior_var,
+ const Matrix<double>& beta){
+
+ // likelihood
+ double loglike = 0.0;
+ const int n = nY.rows();
+ const int ncat = alpha.rows() + 1;
+
+ // the linear predictor
+ Matrix<> mu = X*beta;
+ Matrix<> gamma = alpha2gamma(alpha);
+
+ // compute prob
+ Matrix<> cat_prob(n, ncat-1);
+ //cat_prob: CATegory specific PROBability
+ //the first col of cat_prob = pnorm(gamma_1 - mu)
+ //thus, the col number is ncat - 1
+ Matrix<> prob(n, ncat);
+ for (int j=0; j< ncat-1; ++j){
+ for (int i=0; i<n ; ++i){
+ cat_prob(i, j) = pnorm(gamma[j+1] - mu[i], 0.0, 1.0);
+ }
+ }
+ prob(_, ncat-1) = 1 - cat_prob(_, ncat-2); // top category
+ prob(_, 0) = cat_prob(_, 0); // bottom category
+ for (int j=1; j<(ncat-1); ++j){ // middle categories are actually differences of cumulatives
+ prob(_, j) = cat_prob(_,j) - cat_prob(_, j-1);
+ }
+
+ // the loglikelihood
+ for (int i=0; i<n; ++i){
+ int ind = (int) nY[i];
+ loglike += std::log(prob(i,ind-1));
+ }
+
+ // prior
+ double logprior = 0.0;
+ logprior = lndmvn_jhp(alpha, alpha_prior_mean, alpha_prior_var);
+
+ return (loglike + logprior);
+
}// end of statis oprobit_logpost
// The oprobitModel class stores additional data to be passed to BFGS
class oprobitModel {
public:
- double operator() (const Matrix<double> alpha){
- const int n = y_.rows();
- const int ncat = alpha.rows() + 1;
-
- // the linear predictor
- Matrix<> mu = X_ * beta_;
- Matrix<> gamma = alpha2gamma(alpha);
-
- Matrix<> cat_prob(n, ncat-1);
- //cat_prob: CATegory specific PROBability
- //the first col of cat_prob = pnorm(gamma_1 - mu)
- //thus, the col number is ncat - 1
- Matrix<> prob(n, ncat);
- for (int j=0; j< ncat-1; ++j){
- for (int i=0; i<n ; ++i){
- cat_prob(i, j) = pnorm(gamma[j+1] - mu[i], 0.0, 1.0);
- }
- }
- prob(_, ncat-1) = 1 - cat_prob(_, ncat-2); // top category
- prob(_, 0) = cat_prob(_, 0); // bottom category
- for (int j=1; j<(ncat-1); ++j){ // middle categories are actually differences of cumulatives
- prob(_, j) = cat_prob(_,j) - cat_prob(_, j-1);
- }
-
- // the loglikelihood
- double loglike = 0.0;
- for (int i=0; i<n; ++i){
- int ind = y_[i];
- loglike += std::log(prob(i,ind-1));
- }
- return -1 * loglike;
- }
- Matrix<double> y_;
- Matrix<double> X_;
- Matrix<double> beta_;
+ double operator() (const Matrix<double> alpha){
+ const int n = y_.rows();
+ const int ncat = alpha.rows() + 1;
+
+ // the linear predictor
+ Matrix<> mu = X_ * beta_;
+ Matrix<> gamma = alpha2gamma(alpha);
+
+ Matrix<> cat_prob(n, ncat-1);
+ //cat_prob: CATegory specific PROBability
+ //the first col of cat_prob = pnorm(gamma_1 - mu)
+ //thus, the col number is ncat - 1
+ Matrix<> prob(n, ncat);
+ for (int j=0; j< ncat-1; ++j){
+ for (int i=0; i<n ; ++i){
+ cat_prob(i, j) = pnorm(gamma[j+1] - mu[i], 0.0, 1.0);
+ }
+ }
+ prob(_, ncat-1) = 1 - cat_prob(_, ncat-2); // top category
+ prob(_, 0) = cat_prob(_, 0); // bottom category
+ for (int j=1; j<(ncat-1); ++j){ // middle categories are actually differences of cumulatives
+ prob(_, j) = cat_prob(_,j) - cat_prob(_, j-1);
+ }
+
+ // the loglikelihood
+ double loglike = 0.0;
+ for (int i=0; i<n; ++i){
+ int ind = y_[i];
+ loglike += std::log(prob(i,ind-1));
+ }
+ return -1 * loglike;
+ }
+ Matrix<double> y_;
+ Matrix<double> X_;
+ Matrix<double> beta_;
};
/* MCMCoprobit implementation. Takes Matrix<> reference which it
@@ -197,167 +196,168 @@ public:
*/
template <typename RNGTYPE>
void MCMCoprobit_impl (rng<RNGTYPE>& stream, const int * Y,
- const Matrix<>& nY, const Matrix<>& X, Matrix<>& beta, Matrix<>& gamma,
- const Matrix<>& b0, const Matrix<>& B0,
- const Matrix<>& alpha_prior_mean, const Matrix<>& alpha_prior_var,
- const unsigned int burnin, const unsigned int mcmc, const unsigned int thin,
- const unsigned int verbose, const Matrix<>& tune, const double tdf,
- const unsigned int cowles, Matrix<>& result) {
-
- // define constants and from cross-product matrices
- const unsigned int tot_iter = burnin + mcmc; // total number of mcmc iterations
- const unsigned int nstore = mcmc / thin; // number of draws to store
- const unsigned int k = X.cols();
- const unsigned int N = X.rows();
- const int ncat = gamma.rows() - 1;
- const Matrix<> XpX = crossprod(X);
-
- // storage matrix or matrices
- Matrix<> storemat(nstore, k+ncat+1);
-
- // initialize Z
- Matrix<> Z(N,1);
- Matrix<> Xbeta = X * beta;
-
- // Gibbs loop
- int count = 0;
- int accepts = 0;
-
- Matrix<> gamma_p = gamma;
- Matrix<> gamma_new = gamma + 1;
- Matrix<> alpha = gamma2alpha(gamma_new);
- Matrix<> alpha_hat = alpha;
-
- // initialize current value stuff
- Matrix<> propV = tune * alpha_prior_var * tune;
- Matrix<> propCinvT = t(cholesky(invpd(propV)));
- double logpost_cur = oprobit_logpost(nY, X, alpha, alpha_prior_mean, alpha_prior_var, beta);
- double logjump_cur = lnmulttdens(alpha_prior_mean, alpha_hat, propCinvT, tdf);
-
- double tune_double = tune[0];
- for (unsigned int iter = 0; iter < tot_iter; ++iter) {
-
- //////////////////
- if (cowles == 1){
- //////////////////
- // Cowles method [gamma | Z, beta]
- for (int i=2; i<(ncat); ++i){
- if (i==(ncat-1)){
- gamma_p[i] = stream.rtbnorm_combo(gamma[i], ::pow(tune_double, 2.0),
- gamma_p[i-1]);
- }
- else {
- gamma_p[i] = stream.rtnorm_combo(gamma[i], ::pow(tune_double, 2.0),
- gamma_p[i-1],
- gamma[i+1]);
- }
- }
- double loglikerat = 0.0;
- double loggendenrat = 0.0;
-
- // loop over observations and construct the acceptance ratio
- for (unsigned int i=0; i<N; ++i){
- if (Y[i] == ncat){
- loglikerat = loglikerat
- + log(1.0 - pnorm(gamma_p[Y[i]-1] - Xbeta[i], 0.0, 1.0) )
- - log(1.0 - pnorm(gamma[Y[i]-1] - Xbeta[i], 0.0, 1.0) );
- }
- else if (Y[i] == 1){
- loglikerat = loglikerat + log(pnorm(gamma_p[Y[i]] - Xbeta[i], 0.0, 1.0) )
- - log(pnorm(gamma[Y[i]] - Xbeta[i], 0.0, 1.0) );
- }
- else{
- loglikerat = loglikerat
- + log(pnorm(gamma_p[Y[i]] - Xbeta[i], 0.0, 1.0) -
- pnorm(gamma_p[Y[i]-1] - Xbeta[i], 0.0, 1.0) )
- - log(pnorm(gamma[Y[i]] - Xbeta[i], 0.0, 1.0) -
- pnorm(gamma[Y[i]-1] - Xbeta[i], 0.0, 1.0) );
- }
- }
- double logacceptrat = loglikerat + loggendenrat;
- if (stream.runif() <= exp(logacceptrat)){
- gamma = gamma_p;
- ++accepts;
- }
-
- }// end of Cowles
-
- //////////////////
+ const Matrix<>& nY, const Matrix<>& X, Matrix<>& beta, Matrix<>& gamma,
+ const Matrix<>& b0, const Matrix<>& B0,
+ const Matrix<>& alpha_prior_mean, const Matrix<>& alpha_prior_var,
+ const unsigned int burnin, const unsigned int mcmc, const unsigned int thin,
+ const unsigned int verbose, const Matrix<>& tune, const double tdf,
+ const unsigned int cowles, Matrix<>& result) {
+
+ // define constants and from cross-product matrices
+ const unsigned int tot_iter = burnin + mcmc; // total number of mcmc iterations
+ const unsigned int nstore = mcmc / thin; // number of draws to store
+ const unsigned int k = X.cols();
+ const unsigned int N = X.rows();
+ const int ncat = gamma.rows() - 1;
+ const Matrix<> XpX = crossprod(X);
+
+ // storage matrix or matrices
+ Matrix<> storemat(nstore, k+ncat+1);
+
+ // initialize Z
+ Matrix<> Z(N,1);
+ Matrix<> Xbeta = X * beta;
+
+ // Gibbs loop
+ int count = 0;
+ int accepts = 0;
+
+ Matrix<> gamma_p = gamma;
+ Matrix<> gamma_new = gamma + 1;
+ Matrix<> alpha = gamma2alpha(gamma_new);
+ Matrix<> alpha_hat = alpha;
+
+ // initialize current value stuff
+ Matrix<> propV = tune * alpha_prior_var * tune;
+ Matrix<> propCinvT = t(cholesky(invpd(propV)));
+ double logpost_cur = oprobit_logpost(nY, X, alpha, alpha_prior_mean, alpha_prior_var, beta);
+ double logjump_cur = lnmulttdens(alpha_prior_mean, alpha_hat, propCinvT, tdf);
+
+ double tune_double = tune[0];
+ for (unsigned int iter = 0; iter < tot_iter; ++iter) {
+
+ //////////////////
+ if (cowles == 1){
+ //////////////////
+ // Cowles method [gamma | Z, beta]
+ for (int i=2; i<(ncat); ++i){
+ if (i==(ncat-1)){
+ gamma_p[i] = stream.rtbnorm_combo(gamma[i], ::pow(tune_double, 2.0),
+ gamma_p[i-1]);
+ }
else {
- //////////////////
- // Albert and Chib (2001) method
- // Step 1: [gamma | Z, beta]
- oprobitModel oprobit_model;
- oprobit_model.y_ = nY;
- oprobit_model.X_ = X;
- oprobit_model.beta_ = beta;
-
- // tailored proposal MH
- Matrix<> alpha_hat = BFGS(oprobit_model, alpha, stream, 100, 1e-5, false);
- Matrix<> alpha_V = invpd(hesscdif(oprobit_model, alpha_hat)); //note that oprobit_model contains the multiplication by -1
- Matrix<> propV = tune * alpha_V * tune;
- Matrix<> propCinvT = ::t(cholesky(invpd(propV)));
-
- // Draw alpha_can from multivariate t
- Matrix<> alpha_can = alpha_hat + stream.rmvt(propV, tdf);
-
- // compute components of transition kernel
- double logpost_can = oprobit_logpost(nY, X, alpha_can, alpha_prior_mean, alpha_prior_var, beta);
- double logjump_can = lnmulttdens(alpha_can, alpha_hat, propCinvT, tdf);
- double logpost_cur = oprobit_logpost(nY, X, alpha, alpha_prior_mean, alpha_prior_var, beta);
- double logjump_cur = lnmulttdens(alpha, alpha_hat, propCinvT, tdf);
- double ratio = exp(logpost_can - logjump_can - logpost_cur + logjump_cur);
- const double u = stream();
- if (u < ratio) {
- alpha = alpha_can;
- gamma = alpha2gamma(alpha);
- logpost_cur = logpost_can;
- logjump_cur = logjump_can;
- ++accepts;
- }
- }// end of AC method
-
- // Step 2: [Z| gamma, beta, y]
- Xbeta = X * beta;
- // Matrix<> Z_noconst(N,1);
- for (unsigned int i=0; i<N; ++i){
- Z[i] = stream.rtnorm_combo(Xbeta[i], 1.0, gamma[Y[i]-1], gamma[Y[i]]);
- }
-
- // Step 3: [beta|Z, gamma]
- const Matrix<> XpZ = t(X) * Z;
- beta = NormNormregress_beta_draw(XpX, XpZ, b0, B0, 1.0, stream);
-
- // store values in matrices
- if (iter >= burnin && ((iter % thin)==0)){
- for (unsigned int j=0; j<k; ++j)
- storemat(count, j) = beta[j];
- for (int j=0; j<(ncat+1); ++j)
- storemat(count, j+k) = gamma[j];
- ++count;
- }
-
- // print output to stdout
- if(verbose > 0 && iter % verbose == 0){
- Rprintf("\n\nMCMCoprobit iteration %i of %i \n", (iter+1), tot_iter);
- Rprintf("beta = \n");
- Rprintf("%10.5f\n", beta[0]-gamma[1]);
- for (unsigned int j=1; j<k; ++j)
- Rprintf("%10.5f\n", beta[j]);
- Rprintf("Metropolis acceptance rate for gamma = %3.5f\n\n",
- static_cast<double>(accepts)/static_cast<double>(iter+1));
- }
-
- R_CheckUserInterrupt(); // allow user interrupts
- }// end of MCMC
-
- Rprintf("\n\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n");
- Rprintf("The Metropolis acceptance rate for beta was %3.5f",
- static_cast<double>(accepts) / static_cast<double>(tot_iter));
- Rprintf("\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n");
-
- result = storemat;
-
+ gamma_p[i] = stream.rtnorm_combo(gamma[i], ::pow(tune_double, 2.0),
+ gamma_p[i-1],
+ gamma[i+1]);
+ }
+ }
+ double loglikerat = 0.0;
+ double loggendenrat = 0.0;
+
+ // loop over observations and construct the acceptance ratio
+ for (unsigned int i=0; i<N; ++i){
+ if (Y[i] == ncat){
+ loglikerat = loglikerat
+ + log(1.0 - pnorm(gamma_p[Y[i]-1] - Xbeta[i], 0.0, 1.0) )
+ - log(1.0 - pnorm(gamma[Y[i]-1] - Xbeta[i], 0.0, 1.0) );
+ }
+ else if (Y[i] == 1){
+ loglikerat = loglikerat + log(pnorm(gamma_p[Y[i]] - Xbeta[i], 0.0, 1.0) )
+ - log(pnorm(gamma[Y[i]] - Xbeta[i], 0.0, 1.0) );
+ }
+ else{
+ loglikerat = loglikerat
+ + log(pnorm(gamma_p[Y[i]] - Xbeta[i], 0.0, 1.0) -
+ pnorm(gamma_p[Y[i]-1] - Xbeta[i], 0.0, 1.0) )
+ - log(pnorm(gamma[Y[i]] - Xbeta[i], 0.0, 1.0) -
+ pnorm(gamma[Y[i]-1] - Xbeta[i], 0.0, 1.0) );
+ }
+ }
+ double logacceptrat = loglikerat + loggendenrat;
+ if (stream.runif() <= exp(logacceptrat)){
+ gamma = gamma_p;
+ ++accepts;
+ }
+
+ }// end of Cowles
+
+ //////////////////
+ else {
+ //////////////////
+ // Albert and Chib (2001) method
+ // Step 1: [gamma | Z, beta]
+ oprobitModel oprobit_model;
+ oprobit_model.y_ = nY;
+ oprobit_model.X_ = X;
+ oprobit_model.beta_ = beta;
+
+ // tailored proposal MH
+ Matrix<> alpha_hat = BFGS(oprobit_model, alpha, stream, 100, 1e-5, false);
+ Matrix<> alpha_V = invpd(hesscdif(oprobit_model, alpha_hat));
+ //note that oprobit_model contains the multiplication by -1
+ Matrix<> propV = tune * alpha_V * tune;
+ Matrix<> propCinvT = ::t(cholesky(invpd(propV)));
+
+ // Draw alpha_can from multivariate t
+ Matrix<> alpha_can = alpha_hat + stream.rmvt(propV, tdf);
+
+ // compute components of transition kernel
+ double logpost_can = oprobit_logpost(nY, X, alpha_can, alpha_prior_mean, alpha_prior_var, beta);
+ double logjump_can = lnmulttdens(alpha_can, alpha_hat, propCinvT, tdf);
+ double logpost_cur = oprobit_logpost(nY, X, alpha, alpha_prior_mean, alpha_prior_var, beta);
+ double logjump_cur = lnmulttdens(alpha, alpha_hat, propCinvT, tdf);
+ double ratio = exp(logpost_can - logjump_can - logpost_cur + logjump_cur);
+ const double u = stream();
+ if (u < ratio) {
+ alpha = alpha_can;
+ gamma = alpha2gamma(alpha);
+ logpost_cur = logpost_can;
+ logjump_cur = logjump_can;
+ ++accepts;
+ }
+ }// end of AC method
+
+ // Step 2: [Z| gamma, beta, y]
+ Xbeta = X * beta;
+ // Matrix<> Z_noconst(N,1);
+ for (unsigned int i=0; i<N; ++i){
+ Z[i] = stream.rtnorm_combo(Xbeta[i], 1.0, gamma[Y[i]-1], gamma[Y[i]]);
+ }
+
+ // Step 3: [beta|Z, gamma]
+ const Matrix<> XpZ = t(X) * Z;
+ beta = NormNormregress_beta_draw(XpX, XpZ, b0, B0, 1.0, stream);
+
+ // store values in matrices
+ if (iter >= burnin && ((iter % thin)==0)){
+ for (unsigned int j=0; j<k; ++j)
+ storemat(count, j) = beta[j];
+ for (int j=0; j<(ncat+1); ++j)
+ storemat(count, j+k) = gamma[j];
+ ++count;
+ }
+
+ // print output to stdout
+ if(verbose > 0 && iter % verbose == 0){
+ Rprintf("\n\nMCMCoprobit iteration %i of %i \n", (iter+1), tot_iter);
+ Rprintf("beta = \n");
+ Rprintf("%10.5f\n", beta[0]-gamma[1]);
+ for (unsigned int j=1; j<k; ++j)
+ Rprintf("%10.5f\n", beta[j]);
+ Rprintf("Metropolis acceptance rate for gamma = %3.5f\n\n",
+ static_cast<double>(accepts)/static_cast<double>(iter+1));
+ }
+
+ R_CheckUserInterrupt(); // allow user interrupts
+ }// end of MCMC
+
+ Rprintf("\n\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n");
+ Rprintf("The Metropolis acceptance rate for beta was %3.5f",
+ static_cast<double>(accepts) / static_cast<double>(tot_iter));
+ Rprintf("\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n");
+
+ result = storemat;
+
}
@@ -389,20 +389,20 @@ extern "C"{
Matrix <> gamma(*gammarow, *gammacol, gammadata);
const Matrix <> b0(*b0row, *b0col, b0data);
const Matrix <> B0(*B0row, *B0col, B0data);
- const Matrix <> alpha_prior_mean(*a0row, *a0col, a0data);
+ const Matrix <> alpha_prior_mean(*a0row, *a0col, a0data);
const Matrix <> alpha_prior_var(*A0row, *A0col, A0data);
const Matrix<> tune(*tunerow, *tunecol, tunedata);
-
+
Matrix<> storagematrix;
MCMCPACK_PASSRNG2MODEL(MCMCoprobit_impl, Y, nY, X, beta, gamma, b0, B0,
- alpha_prior_mean, alpha_prior_var,
- *burnin, *mcmc, *thin, *verbose, tune, *tdf,
- *cowles, storagematrix);
+ alpha_prior_mean, alpha_prior_var,
+ *burnin, *mcmc, *thin, *verbose, tune, *tdf,
+ *cowles, storagematrix);
const unsigned int size = *samplerow * *samplecol;
for (unsigned int i = 0; i < size; ++i)
- sampledata[i] = storagematrix(i);
- }
+ sampledata[i] = storagematrix(i);
+ }
}
#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