[r-cran-mcmcpack] 34/90: Imported Upstream version 1.0-1
Andreas Tille
tille at debian.org
Fri Dec 16 09:07:37 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 622e43a0ad4b580d57241714048e9fe7d6c0ee78
Author: Andreas Tille <tille at debian.org>
Date: Fri Dec 16 08:07:13 2016 +0100
Imported Upstream version 1.0-1
---
DESCRIPTION | 32 ++-
HISTORY | 1 +
NAMESPACE | 1 +
R/MCMCdynamicIRT1d-b.R | 417 +++++++++++++++++++++++++++
R/MCMCdynamicIRT1d.R | 415 +++++++++++++++++++++++++++
data/Rehnquist.rda | Bin 0 -> 1696 bytes
man/MCMCdynamicIRT1d.Rd | 252 ++++++++++++++++
man/MCMCirtKd.Rd | 2 +-
man/Rehnquist.Rd | 32 +++
src/MCMCdynamicIRT1d-b.cc | 714 ++++++++++++++++++++++++++++++++++++++++++++++
src/MCMCdynamicIRT1d.cc | 699 +++++++++++++++++++++++++++++++++++++++++++++
11 files changed, 2549 insertions(+), 16 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index 7953237..4190fa1 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,22 +1,24 @@
Package: MCMCpack
-Version: 0.9-6
-Date: 2009-2-17
+Version: 1.0-1
+Date: 2009-6-29
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>
+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
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++ written in the Scythe
- Statistical Library Version 1.0.2. All models return coda mcmc objects
- that can then be summarized using the coda package. MCMCpack
- also contains some useful utility functions, including some
- additional density functions and pseudo-random number generators
- for statistical distributions, a general purpose Metropolis
- sampling algorithm, and tools for visualization.
-License: GPL version 2
+ inference using posterior simulation for a number of
+ statistical models. Most simulation is done in compiled C++
+ written in the Scythe Statistical Library Version 1.0.2. All
+ models return coda mcmc objects that can then be summarized
+ using the coda package. MCMCpack also contains some useful
+ utility functions, including some additional density functions
+ and pseudo-random number generators for statistical
+ distributions, a general purpose Metropolis sampling algorithm,
+ and tools for visualization.
+License: GPL-2
SystemRequirements: gcc (>= 4.0)
URL: http://mcmcpack.wustl.edu
-Packaged: Tue Feb 17 10:39:42 2009; adm
+Packaged: 2009-06-30 09:23:24 UTC; adm
+Repository: CRAN
+Date/Publication: 2009-07-01 12:02:06
diff --git a/HISTORY b/HISTORY
index d88800d..94a4ebf 100644
--- a/HISTORY
+++ b/HISTORY
@@ -4,6 +4,7 @@
0.9-6 to 1.0-1
* added one-dimensional dynamic IRT model
+ * added Rehnquist Court data to illustrate dynamic IRT model
0.9-5 to 0.9-6
* fixed bug in MCMCmetrop1R.R [thanks to many for spotting memory issue]
diff --git a/NAMESPACE b/NAMESPACE
index 57c881e..eade73c 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -17,6 +17,7 @@ export(
MCnormalnormal,
MCmultinomdirichlet,
MCMCdynamicEI,
+ MCMCdynamicIRT1d,
MCMCfactanal,
MCMChierEI,
MCMCirt1d,
diff --git a/R/MCMCdynamicIRT1d-b.R b/R/MCMCdynamicIRT1d-b.R
new file mode 100644
index 0000000..78a25b1
--- /dev/null
+++ b/R/MCMCdynamicIRT1d-b.R
@@ -0,0 +1,417 @@
+##########################################################################
+## samples from the posterior distribution of a dynamic 1d IRT model
+## a la Martin and Quinn (2002)
+##
+## 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.
+##
+## Assumes a local level model for the evolution of theta
+##
+## y_{jkt}^* = -alpha_k + beta_k * theta_{jt} + epsilon_{jkt}
+## theta_{jt} ~ N(theta_{j(t-1)}, tau^2)
+##
+## Kevin Quinn
+## 1/28/2008
+##
+## 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
+##########################################################################
+
+## prior for alpha and beta that coincides with a uniform prior on the
+## cutpoints
+
+"MCMCdynamicIRT1d_b" <- function(datamatrix, item.time.map,
+ theta.constraints=list(),
+ burnin=1000, mcmc=20000, thin=1,
+ verbose=0, seed=NA, theta.start=NA,
+ alpha.start=NA, beta.start=NA,
+ tau2.start=1,
+ a0=0, A0=.1, b0=0, B0=.1,
+ c0=-1, d0=-1, e0=0, E0=1,
+ store.ability=TRUE,
+ store.item=TRUE, ...
+ ){
+
+
+ datamatrix <- as.matrix(datamatrix)
+
+ nitems <- ncol(datamatrix)
+ nsubj <- nrow(datamatrix)
+ ntime <- max(item.time.map)
+
+ ## checks
+ check.offset(list(...))
+ check.mcmc.parameters(burnin, mcmc, thin)
+
+
+ if (nitems != length(item.time.map)){
+ cat("Number of rows of datamatrix not equal to length of item.time.map\n")
+ stop("Please check data and try MCMCdynamicIRT1d() again.\n",
+ call.=FALSE)
+ }
+ if (min(item.time.map) != 1){
+ cat("Minimum value in item.time.map not equal to 1\n")
+ stop("Please check data and try MCMCdynamicIRT1d() again.\n",
+ call.=FALSE)
+ }
+ if(sum(datamatrix==1 | datamatrix==0 | is.na(datamatrix)) !=
+ (nitems * nsubj)) {
+ cat("Error: Data matrix contains elements other than 0, 1 or NA.\n")
+ stop("Please check data and try MCMCdynamicIRT1d() again.\n",
+ call.=FALSE)
+ }
+
+
+ if (A0 < 0){
+ cat("Error: A0 (prior precision for alpha) is less than 0).\n")
+ stop("Please respecify and try MCMCdynamicIRT1d() again.\n")
+ }
+ if (B0 < 0){
+ cat("Error: B0 (prior precision for beta) is less than 0).\n")
+ stop("Please respecify and try MCMCdynamicIRT1d() again.\n")
+ }
+
+
+
+
+
+
+ ## setup constraints on theta
+ if(length(theta.constraints) != 0) {
+ for (i in 1:length(theta.constraints)){
+ theta.constraints[[i]] <-
+ list(as.integer(1), theta.constraints[[i]][1])
+ }
+ }
+ holder <- build.factor.constraints(theta.constraints, t(datamatrix), nsubj, 1)
+ theta.eq.constraints <- holder[[1]]
+ theta.ineq.constraints <- holder[[2]]
+ ##subject.names <- holder[[3]]
+ ## names
+ item.names <- colnames(datamatrix)
+ if (is.null(item.names)){
+ item.names <- paste("item", 1:nitems, sep="")
+ }
+
+
+
+
+
+
+ ## starting values for theta error checking
+ theta.start <- factor.score.start.check(theta.start, datamatrix,
+ matrix(0,1,1), matrix(1,1,1),
+ theta.eq.constraints,
+ theta.ineq.constraints, 1)
+
+ ## starting values for (alpha, beta)
+ ab.starts <- matrix(NA, nitems, 2)
+ for (i in 1:nitems){
+ local.y <- datamatrix[,i]
+ ##local.y[local.y==9] <- NA
+ if (length(na.omit(local.y)) <= 1){
+ ab.starts[i,] <- c(0, 10)
+ }
+ else if (var(na.omit(local.y))==0){
+ ab.starts[i,] <- c(0,10)
+
+ }
+ else {
+ ab.starts[i,] <- coef(suppressWarnings(glm(local.y~theta.start,
+ family=binomial(probit),
+ control=glm.control(
+ maxit=8, epsilon=1e-3)
+ )))
+ }
+ }
+ ab.starts[,1] <- -1 * ab.starts[,1] # make this into a difficulty param
+
+ ## starting values for alpha and beta error checking
+ if (is.na(alpha.start)) {
+ alpha.start <- ab.starts[,1]
+ }
+ else if(is.null(dim(alpha.start))) {
+ alpha.start <- alpha.start * matrix(1,nitems,1)
+ }
+ else if((dim(alpha.start)[1] != nitems) || (dim(alpha.start)[2] != 1)) {
+ cat("Error: Starting value for alpha not conformable.\n")
+ stop("Please respecify and call MCMCdynamicIRT1d() again.\n",
+ call.=FALSE)
+ }
+ if (is.na(beta.start)) {
+ beta.start <- ab.starts[,2]
+ }
+ else if(is.null(dim(beta.start))) {
+ beta.start <- beta.start * matrix(1,nitems,1)
+ }
+ else if((dim(beta.start)[1] != nitems) || (dim(beta.start)[2] != 1)) {
+ cat("Error: Starting value for beta not conformable.\n")
+ stop("Please respecify and call MCMCdynamicIRT1d() again.\n",
+ call.=FALSE)
+ }
+
+
+
+
+ ## generate time-specific theta information and create extended theta.start
+ subject.names <- NULL
+ theta.start.new <- NULL
+ ## theta.info has:
+ ## col1: subj ID, col2: #time periods, col3: offset (first term C indexing)
+ theta.info <- matrix(NA, nsubj, 3)
+ for (s in 1:nsubj){
+ namestub <- rownames(datamatrix)[s]
+ theta.info[s,1] <- s
+ count <- 0
+ holder <- NULL
+ for (i in 1:nitems){
+ if (!is.na(datamatrix[s,i])){
+ holder <- c(holder, item.time.map[i])
+ }
+ }
+ minholder <- min(holder)
+ maxholder <- max(holder)
+ theta.info[s,2] <- maxholder - minholder + 1
+ theta.info[s,3] <- minholder - 1
+ theta.start.new <- c(theta.start.new, rep(theta.start[s], theta.info[s,2]))
+ subject.names <- c(subject.names,
+ paste(namestub, ".t", minholder:maxholder, sep=""))
+ }
+ nthetas <- length(subject.names)
+ theta.start <- theta.start.new
+
+
+
+ if (length(c0) < nsubj){
+ c0 <- array(c0, nsubj)
+ }
+ if (length(d0) < nsubj){
+ d0 <- array(d0, nsubj)
+ }
+ if (length(tau2.start) < nsubj){
+ tau2.start <- array(tau2.start, nsubj)
+ }
+ if (length(e0) < nsubj){
+ e0 <- array(e0, nsubj)
+ }
+ if (length(E0) < nsubj){
+ E0 <- array(E0, nsubj)
+ }
+ E0inv <- 1/E0
+
+
+ subject.names.short <- rownames(datamatrix)
+
+
+ ## convert datamatrix into a sparse format where the missing values are
+ ## not explicitly represented
+ ##
+ ## 1st column: subject index (C indexing)
+ ## 2nd column: item index (C indexing)
+ ## 3rd column: vote
+ data.sparse.si <- NULL
+ for (i in 1:nsubj){
+ for (j in 1:nitems){
+ if (!is.na(datamatrix[i,j])){
+ data.sparse.si <- rbind(data.sparse.si, c(i-1, j-1, datamatrix[i,j]))
+ }
+ }
+ }
+ ## 1st column: item index (C indexing)
+ ## 2nd column: subject index (C indexing)
+ ## 3rd column: vote
+## data.sparse.is <- NULL
+## for (i in 1:nitems){
+## for (j in 1:nsubj){
+## if (!is.na(datamatrix[j,i])){
+## data.sparse.is <- rbind(data.sparse.is, c(i-1, j-1, datamatrix[j,i]))
+## }
+## }
+## }
+
+ rm(datamatrix)
+
+
+
+
+
+
+ ## define storage matrix for posterior theta draws
+ thetadraws <- matrix(0, mcmc/thin, nthetas)
+ if (store.ability != TRUE){
+ thetadraws <- matrix(0, 1, 1)
+ }
+ alphadraws <- matrix(1, mcmc/thin, nitems)
+ betadraws <- matrix(2, mcmc/thin, nitems)
+ if (store.item != TRUE){
+ alphadraws <- matrix(1, 1, 1)
+ betadraws <- matrix(2, 1, 1)
+ }
+ tau2draws <- matrix(0, mcmc/thin, nsubj)
+
+
+ ## seeds
+ seeds <- form.seeds(seed)
+ lecuyer <- seeds[[1]]
+ seed.array <- seeds[[2]]
+ lecuyer.stream <- seeds[[3]]
+
+ ## print(theta.eq.constraints)
+ ## print(theta.ineq.constraints)
+
+
+# return(list(theta=theta.start, alpha=alpha.start, beta=beta.start))
+
+
+ ## call C++ code to draw sample
+ posterior <- .C("MCMCdynamicIRT1d_b",
+
+ thetadata = as.double(thetadraws),
+ thetarow = as.integer(nrow(thetadraws)),
+ thetacol = as.integer(ncol(thetadraws)),
+
+ alphadata = as.double(alphadraws),
+ alpharow = as.integer(nrow(alphadraws)),
+ alphacol = as.integer(ncol(alphadraws)),
+
+ betadata = as.double(betadraws),
+ betarow = as.integer(nrow(betadraws)),
+ betacol = as.integer(ncol(betadraws)),
+
+ tau2data = as.double(tau2draws),
+ tau2row = as.integer(nrow(tau2draws)),
+ tau2col = as.integer(ncol(tau2draws)),
+
+ nsubj = as.integer(nsubj),
+ nitems = as.integer(nitems),
+ ntime = as.integer(ntime),
+
+ Ydata = as.integer(data.sparse.si),
+ Yrow = as.integer(nrow(data.sparse.si)),
+ Ycol = as.integer(ncol(data.sparse.si)),
+
+ IT = as.integer(item.time.map-1),
+ ITlength = as.integer(length(item.time.map)),
+
+ 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),
+
+ thetastartdata = as.double(theta.start),
+ thetastartlength = as.integer(length(theta.start)),
+
+ thetainfo = as.integer(theta.info),
+ thetainforow = as.integer(nrow(theta.info)),
+ thetainfocol = as.integer(ncol(theta.info)),
+
+ astartdata = as.double(alpha.start),
+ astartlength = as.integer(length(alpha.start)),
+
+ bstartdata = as.double(beta.start),
+ bstartlength = as.integer(length(beta.start)),
+
+ tau2data = as.double(tau2.start),
+ tau2length = as.integer(length(tau2.start)),
+
+ c0data = as.double(c0),
+ c0length = as.integer(length(c0)),
+
+ d0data = as.double(d0),
+ d0length = as.integer(length(d0)),
+
+ a0data = as.double(a0),
+ A0data = as.double(A0),
+
+ b0data = as.double(b0),
+ B0data = as.double(B0),
+
+ e0data = as.double(e0),
+ E0invdata = as.double(E0inv),
+
+ thetaeqdata = as.double(theta.eq.constraints),
+ thetaeqrow = as.integer(nrow(theta.eq.constraints)),
+ thetaeqcol = as.integer(ncol(theta.eq.constraints)),
+
+ thetaineqdata = as.double(theta.ineq.constraints),
+ thetaineqrow = as.integer(nrow(theta.ineq.constraints)),
+ thetaineqcol = as.integer(ncol(theta.ineq.constraints)),
+
+ storei = as.integer(store.item),
+ storea = as.integer(store.ability),
+ PACKAGE="MCMCpack"
+ )
+
+
+
+
+
+ tau2samp <- matrix(posterior$tau2data,
+ posterior$tau2row,
+ posterior$tau2col)
+ colnames(tau2samp) <- paste("tau2.", subject.names.short, sep="")
+
+
+
+ if (store.item & store.ability){
+ thetasamp <- matrix(posterior$thetadata,
+ posterior$thetarow,
+ posterior$thetacol)
+ colnames(thetasamp) <- paste("theta.", subject.names, sep="")
+
+ alphasamp <- matrix(posterior$alphadata,
+ posterior$alpharow,
+ posterior$alphacol)
+ colnames(alphasamp) <- paste("alpha.", item.names, sep="")
+
+
+ betasamp <- matrix(posterior$betadata,
+ posterior$betarow,
+ posterior$betacol)
+ colnames(betasamp) <- paste("beta.", item.names, sep="")
+
+ outmat <- mcmc(cbind(thetasamp, alphasamp, betasamp, tau2samp),
+ start=1, end=mcmc, thin=thin)
+ }
+ else if (store.item){
+ alphasamp <- matrix(posterior$alphadata,
+ posterior$alpharow,
+ posterior$alphacol)
+ colnames(alphasamp) <- paste("alpha.", item.names, sep="")
+
+
+ betasamp <- matrix(posterior$betadata,
+ posterior$betarow,
+ posterior$betacol)
+ colnames(betasamp) <- paste("beta.", item.names, sep="")
+
+ outmat <- mcmc(cbind(alphasamp, betasamp, tau2samp),
+ start=1, end=mcmc, thin=thin)
+ }
+ else if (store.ability){
+ thetasamp <- matrix(posterior$thetadata,
+ posterior$thetarow,
+ posterior$thetacol)
+ colnames(thetasamp) <- paste("theta.", subject.names, sep="")
+
+ outmat <- mcmc(cbind(thetasamp, tau2samp),
+ start=1, end=mcmc, thin=thin)
+ }
+ else {
+ outmat <- mcmc(tau2samp,
+ start=1, end=mcmc, thin=thin)
+ }
+
+ return(outmat)
+
+
+}
+
+
diff --git a/R/MCMCdynamicIRT1d.R b/R/MCMCdynamicIRT1d.R
new file mode 100644
index 0000000..590aa48
--- /dev/null
+++ b/R/MCMCdynamicIRT1d.R
@@ -0,0 +1,415 @@
+##########################################################################
+## samples from the posterior distribution of a dynamic 1d IRT model
+## a la Martin and Quinn (2002)
+##
+## 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.
+##
+## Assumes a local level model for the evolution of theta
+##
+## y_{jkt}^* = -alpha_k + beta_k * theta_{jt} + epsilon_{jkt}
+## theta_{jt} ~ N(theta_{j(t-1)}, tau^2)
+##
+## Kevin Quinn
+## 1/28/2008
+##
+## 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
+##########################################################################
+
+
+"MCMCdynamicIRT1d" <- function(datamatrix, item.time.map,
+ theta.constraints=list(),
+ burnin=1000, mcmc=20000, thin=1,
+ verbose=0, seed=NA, theta.start=NA,
+ alpha.start=NA, beta.start=NA,
+ tau2.start=1,
+ a0=0, A0=.1, b0=0, B0=.1,
+ c0=-1, d0=-1, e0=0, E0=1,
+ store.ability=TRUE,
+ store.item=TRUE, ...
+ ){
+
+
+ datamatrix <- as.matrix(datamatrix)
+
+ nitems <- ncol(datamatrix)
+ nsubj <- nrow(datamatrix)
+ ntime <- max(item.time.map)
+
+ ## checks
+ check.offset(list(...))
+ check.mcmc.parameters(burnin, mcmc, thin)
+
+
+ if (nitems != length(item.time.map)){
+ cat("Number of rows of datamatrix not equal to length of item.time.map\n")
+ stop("Please check data and try MCMCdynamicIRT1d() again.\n",
+ call.=FALSE)
+ }
+ if (min(item.time.map) != 1){
+ cat("Minimum value in item.time.map not equal to 1\n")
+ stop("Please check data and try MCMCdynamicIRT1d() again.\n",
+ call.=FALSE)
+ }
+ if(sum(datamatrix==1 | datamatrix==0 | is.na(datamatrix)) !=
+ (nitems * nsubj)) {
+ cat("Error: Data matrix contains elements other than 0, 1 or NA.\n")
+ stop("Please check data and try MCMCdynamicIRT1d() again.\n",
+ call.=FALSE)
+ }
+
+
+ if (A0 < 0){
+ cat("Error: A0 (prior precision for alpha) is less than 0).\n")
+ stop("Please respecify and try MCMCdynamicIRT1d() again.\n")
+ }
+ if (B0 < 0){
+ cat("Error: B0 (prior precision for beta) is less than 0).\n")
+ stop("Please respecify and try MCMCdynamicIRT1d() again.\n")
+ }
+
+
+
+
+
+
+ ## setup constraints on theta
+ if(length(theta.constraints) != 0) {
+ for (i in 1:length(theta.constraints)){
+ theta.constraints[[i]] <-
+ list(as.integer(1), theta.constraints[[i]][1])
+ }
+ }
+ holder <- build.factor.constraints(theta.constraints, t(datamatrix), nsubj, 1)
+ theta.eq.constraints <- holder[[1]]
+ theta.ineq.constraints <- holder[[2]]
+ ##subject.names <- holder[[3]]
+ ## names
+ item.names <- colnames(datamatrix)
+ if (is.null(item.names)){
+ item.names <- paste("item", 1:nitems, sep="")
+ }
+
+
+
+
+
+
+ ## starting values for theta error checking
+ theta.start <- factor.score.start.check(theta.start, datamatrix,
+ matrix(0,1,1), matrix(1,1,1),
+ theta.eq.constraints,
+ theta.ineq.constraints, 1)
+
+ ## starting values for (alpha, beta)
+ ab.starts <- matrix(NA, nitems, 2)
+ for (i in 1:nitems){
+ local.y <- datamatrix[,i]
+ ##local.y[local.y==9] <- NA
+ if (length(na.omit(local.y)) <= 1){
+ ab.starts[i,] <- c(0, 10)
+ }
+ else if (var(na.omit(local.y))==0){
+ ab.starts[i,] <- c(0,10)
+
+ }
+ else {
+ ab.starts[i,] <- coef(suppressWarnings(glm(local.y~theta.start,
+ family=binomial(probit),
+ control=glm.control(
+ maxit=8, epsilon=1e-3)
+ )))
+ }
+ }
+ ab.starts[,1] <- -1 * ab.starts[,1] # make this into a difficulty param
+
+ ## starting values for alpha and beta error checking
+ if (is.na(alpha.start)) {
+ alpha.start <- ab.starts[,1]
+ }
+ else if(is.null(dim(alpha.start))) {
+ alpha.start <- alpha.start * matrix(1,nitems,1)
+ }
+ else if((dim(alpha.start)[1] != nitems) || (dim(alpha.start)[2] != 1)) {
+ cat("Error: Starting value for alpha not conformable.\n")
+ stop("Please respecify and call MCMCdynamicIRT1d() again.\n",
+ call.=FALSE)
+ }
+ if (is.na(beta.start)) {
+ beta.start <- ab.starts[,2]
+ }
+ else if(is.null(dim(beta.start))) {
+ beta.start <- beta.start * matrix(1,nitems,1)
+ }
+ else if((dim(beta.start)[1] != nitems) || (dim(beta.start)[2] != 1)) {
+ cat("Error: Starting value for beta not conformable.\n")
+ stop("Please respecify and call MCMCdynamicIRT1d() again.\n",
+ call.=FALSE)
+ }
+
+
+
+
+ ## generate time-specific theta information and create extended theta.start
+ subject.names <- NULL
+ theta.start.new <- NULL
+ ## theta.info has:
+ ## col1: subj ID, col2: #time periods, col3: offset (first term C indexing)
+ theta.info <- matrix(NA, nsubj, 3)
+ for (s in 1:nsubj){
+ namestub <- rownames(datamatrix)[s]
+ theta.info[s,1] <- s
+ count <- 0
+ holder <- NULL
+ for (i in 1:nitems){
+ if (!is.na(datamatrix[s,i])){
+ holder <- c(holder, item.time.map[i])
+ }
+ }
+ minholder <- min(holder)
+ maxholder <- max(holder)
+ theta.info[s,2] <- maxholder - minholder + 1
+ theta.info[s,3] <- minholder - 1
+ theta.start.new <- c(theta.start.new, rep(theta.start[s], theta.info[s,2]))
+ subject.names <- c(subject.names,
+ paste(namestub, ".t", minholder:maxholder, sep=""))
+ }
+ nthetas <- length(subject.names)
+ theta.start <- theta.start.new
+
+
+
+ if (length(c0) < nsubj){
+ c0 <- array(c0, nsubj)
+ }
+ if (length(d0) < nsubj){
+ d0 <- array(d0, nsubj)
+ }
+ if (length(tau2.start) < nsubj){
+ tau2.start <- array(tau2.start, nsubj)
+ }
+ if (length(e0) < nsubj){
+ e0 <- array(e0, nsubj)
+ }
+ if (length(E0) < nsubj){
+ E0 <- array(E0, nsubj)
+ }
+ E0inv <- 1/E0
+
+
+ subject.names.short <- rownames(datamatrix)
+
+
+ ## convert datamatrix into a sparse format where the missing values are
+ ## not explicitly represented
+ ##
+ ## 1st column: subject index (C indexing)
+ ## 2nd column: item index (C indexing)
+ ## 3rd column: vote
+ data.sparse.si <- NULL
+ for (i in 1:nsubj){
+ for (j in 1:nitems){
+ if (!is.na(datamatrix[i,j])){
+ data.sparse.si <- rbind(data.sparse.si, c(i-1, j-1, datamatrix[i,j]))
+ }
+ }
+ }
+ ## 1st column: item index (C indexing)
+ ## 2nd column: subject index (C indexing)
+ ## 3rd column: vote
+## data.sparse.is <- NULL
+## for (i in 1:nitems){
+## for (j in 1:nsubj){
+## if (!is.na(datamatrix[j,i])){
+## data.sparse.is <- rbind(data.sparse.is, c(i-1, j-1, datamatrix[j,i]))
+## }
+## }
+## }
+
+ rm(datamatrix)
+
+
+
+
+
+
+ ## define storage matrix for posterior theta draws
+ thetadraws <- matrix(0, mcmc/thin, nthetas)
+ if (store.ability != TRUE){
+ thetadraws <- matrix(0, 1, 1)
+ }
+ alphadraws <- matrix(1, mcmc/thin, nitems)
+ betadraws <- matrix(2, mcmc/thin, nitems)
+ if (store.item != TRUE){
+ alphadraws <- matrix(1, 1, 1)
+ betadraws <- matrix(2, 1, 1)
+ }
+ tau2draws <- matrix(0, mcmc/thin, nsubj)
+
+
+ ## seeds
+ seeds <- form.seeds(seed)
+ lecuyer <- seeds[[1]]
+ seed.array <- seeds[[2]]
+ lecuyer.stream <- seeds[[3]]
+
+ ## print(theta.eq.constraints)
+ ## print(theta.ineq.constraints)
+
+
+# return(list(theta=theta.start, alpha=alpha.start, beta=beta.start))
+
+
+ ## call C++ code to draw sample
+ posterior <- .C("MCMCdynamicIRT1d",
+
+ thetadata = as.double(thetadraws),
+ thetarow = as.integer(nrow(thetadraws)),
+ thetacol = as.integer(ncol(thetadraws)),
+
+ alphadata = as.double(alphadraws),
+ alpharow = as.integer(nrow(alphadraws)),
+ alphacol = as.integer(ncol(alphadraws)),
+
+ betadata = as.double(betadraws),
+ betarow = as.integer(nrow(betadraws)),
+ betacol = as.integer(ncol(betadraws)),
+
+ tau2data = as.double(tau2draws),
+ tau2row = as.integer(nrow(tau2draws)),
+ tau2col = as.integer(ncol(tau2draws)),
+
+ nsubj = as.integer(nsubj),
+ nitems = as.integer(nitems),
+ ntime = as.integer(ntime),
+
+ Ydata = as.integer(data.sparse.si),
+ Yrow = as.integer(nrow(data.sparse.si)),
+ Ycol = as.integer(ncol(data.sparse.si)),
+
+ IT = as.integer(item.time.map-1),
+ ITlength = as.integer(length(item.time.map)),
+
+ 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),
+
+ thetastartdata = as.double(theta.start),
+ thetastartlength = as.integer(length(theta.start)),
+
+ thetainfo = as.integer(theta.info),
+ thetainforow = as.integer(nrow(theta.info)),
+ thetainfocol = as.integer(ncol(theta.info)),
+
+ astartdata = as.double(alpha.start),
+ astartlength = as.integer(length(alpha.start)),
+
+ bstartdata = as.double(beta.start),
+ bstartlength = as.integer(length(beta.start)),
+
+ tau2data = as.double(tau2.start),
+ tau2length = as.integer(length(tau2.start)),
+
+ c0data = as.double(c0),
+ c0length = as.integer(length(c0)),
+
+ d0data = as.double(d0),
+ d0length = as.integer(length(d0)),
+
+ a0data = as.double(a0),
+ A0data = as.double(A0),
+
+ b0data = as.double(b0),
+ B0data = as.double(B0),
+
+ e0data = as.double(e0),
+ E0invdata = as.double(E0inv),
+
+ thetaeqdata = as.double(theta.eq.constraints),
+ thetaeqrow = as.integer(nrow(theta.eq.constraints)),
+ thetaeqcol = as.integer(ncol(theta.eq.constraints)),
+
+ thetaineqdata = as.double(theta.ineq.constraints),
+ thetaineqrow = as.integer(nrow(theta.ineq.constraints)),
+ thetaineqcol = as.integer(ncol(theta.ineq.constraints)),
+
+ storei = as.integer(store.item),
+ storea = as.integer(store.ability),
+ PACKAGE="MCMCpack"
+ )
+
+
+
+
+
+ tau2samp <- matrix(posterior$tau2data,
+ posterior$tau2row,
+ posterior$tau2col)
+ colnames(tau2samp) <- paste("tau2.", subject.names.short, sep="")
+
+
+
+ if (store.item & store.ability){
+ thetasamp <- matrix(posterior$thetadata,
+ posterior$thetarow,
+ posterior$thetacol)
+ colnames(thetasamp) <- paste("theta.", subject.names, sep="")
+
+ alphasamp <- matrix(posterior$alphadata,
+ posterior$alpharow,
+ posterior$alphacol)
+ colnames(alphasamp) <- paste("alpha.", item.names, sep="")
+
+
+ betasamp <- matrix(posterior$betadata,
+ posterior$betarow,
+ posterior$betacol)
+ colnames(betasamp) <- paste("beta.", item.names, sep="")
+
+ outmat <- mcmc(cbind(thetasamp, alphasamp, betasamp, tau2samp),
+ start=1, end=mcmc, thin=thin)
+ }
+ else if (store.item){
+ alphasamp <- matrix(posterior$alphadata,
+ posterior$alpharow,
+ posterior$alphacol)
+ colnames(alphasamp) <- paste("alpha.", item.names, sep="")
+
+
+ betasamp <- matrix(posterior$betadata,
+ posterior$betarow,
+ posterior$betacol)
+ colnames(betasamp) <- paste("beta.", item.names, sep="")
+
+ outmat <- mcmc(cbind(alphasamp, betasamp, tau2samp),
+ start=1, end=mcmc, thin=thin)
+ }
+ else if (store.ability){
+ thetasamp <- matrix(posterior$thetadata,
+ posterior$thetarow,
+ posterior$thetacol)
+ colnames(thetasamp) <- paste("theta.", subject.names, sep="")
+
+ outmat <- mcmc(cbind(thetasamp, tau2samp),
+ start=1, end=mcmc, thin=thin)
+ }
+ else {
+ outmat <- mcmc(tau2samp,
+ start=1, end=mcmc, thin=thin)
+ }
+
+ return(outmat)
+
+
+}
+
+
diff --git a/data/Rehnquist.rda b/data/Rehnquist.rda
new file mode 100644
index 0000000..073d6c0
Binary files /dev/null and b/data/Rehnquist.rda differ
diff --git a/man/MCMCdynamicIRT1d.Rd b/man/MCMCdynamicIRT1d.Rd
new file mode 100644
index 0000000..fdda299
--- /dev/null
+++ b/man/MCMCdynamicIRT1d.Rd
@@ -0,0 +1,252 @@
+\name{MCMCdynamicIRT1d}
+\alias{MCMCdynamicIRT1d}
+\alias{MCMCdynamicIRT1d_b}
+\title{Markov Chain Monte Carlo for Dynamic One Dimensional Item Response Theory Model}
+
+\description{ This function generates a sample from the posterior
+ distribution of a dynamic one dimensional item response theory (IRT)
+ model, with Normal random walk priors on the subject abilities (ideal
+ points), and multivariate Normal priors on the item parameters. 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{
+MCMCdynamicIRT1d(datamatrix, item.time.map, theta.constraints=list(),
+ burnin = 1000, mcmc = 20000, thin = 1, verbose = 0,
+ seed = NA, theta.start = NA, alpha.start = NA,
+ beta.start = NA, tau2.start = 1, a0 = 0, A0 = 0.1,
+ b0 = 0, B0 = 0.1, c0 = -1, d0 = -1, e0 = 0,
+ E0 = 1, store.ability = TRUE,
+ store.item = TRUE, ...)
+}
+
+\arguments{
+ \item{datamatrix}{The matrix of data. Must be 0, 1, or missing values.
+ The rows of \code{datamatrix} correspond to subjects and the
+ columns correspond to items.}
+
+ \item{item.time.map}{A vector that relates each item to a time period.
+ Each element of \code{item.time.map} gives the time period of the
+ corresponding column of \code{datamatrix}. It is assumed that the
+ minimum value of \code{item.time.map} is 1.}
+
+ \item{theta.constraints}{ A list specifying possible simple equality
+ or inequality constraints on the ability parameters. A typical
+ entry in the list has one of three forms: \code{varname=c} which
+ will constrain the ability parameter for the subject named
+ \code{varname} to be equal to c, \code{varname="+"} which will
+ constrain the ability parameter for the subject named \code{varname}
+ to be positive, and \code{varname="-"} which will constrain the
+ ability parameter for the subject named \code{varname} to be
+ negative. If x is a matrix without row names defaults names of
+ ``V1",``V2", ... , etc will be used. See Rivers (2003) for a
+ thorough discussion of identification of IRT models. }
+
+ \item{burnin}{ The number of burn-in iterations for the sampler. }
+
+ \item{mcmc}{The number of Gibbs iterations for the sampler. }
+
+ \item{thin}{The thinning interval used in the simulation. The number of
+ Gibbs 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{theta.start}{ The starting values for the subject
+ abilities (ideal points). This can either be a scalar or a
+ column vector with dimension equal to the number of voters.
+ If this takes a scalar value, then that value will serve as
+ the starting value for all of the thetas. The default value
+ of NA will choose the starting values based on an
+ eigenvalue-eigenvector decomposition of the aggreement score matrix
+ formed from the \code{datamatrix}. }
+
+ \item{alpha.start}{ The starting values for the
+ \eqn{\alpha}{alpha} difficulty parameters. This can either be
+ a scalar or a column vector with dimension equal to the
+ number of items. If this takes a scalar value, then that
+ value will serve as the starting value for all of the alphas.
+ The default value of NA will set the starting values based on
+ a series of probit regressions that condition on the starting
+ values of theta. }
+
+ \item{beta.start}{ The starting values for the
+ \eqn{\beta}{beta} discrimination parameters. This can either
+ be a scalar or a column vector with dimension equal to the
+ number of items. If this takes a scalar value, then that
+ value will serve as the starting value for all of the betas.
+ The default value of NA will set the starting values based on a
+ series of probit regressions that condition on the starting values
+ of theta. }
+
+ \item{tau2.start}{ The starting values for the evolution variances
+ (the variance of the random walk increments for the ability
+ parameters / ideal points. Order corresponds to the rows of
+ \code{datamatrix}.}
+
+ \item{a0}{ A vector containing the prior mean of each of the
+ difficulty parameters \eqn{\alpha}{alpha}. Should have as many
+ elements as items / roll calls. Order corresponds to the columns of
+ \code{datamatrix}. If a scalar is passed it is assumed that all
+ elements of \code{a0} are equal to the scalar. }
+
+ \item{A0}{ A vector containing the prior precision (inverse
+ variance) of each of the difficulty parameters \eqn{\alpha}{alpha}.
+ Should have as many elements as items / roll calls. Order
+ corresponds to the columns of \code{datamatrix}. If a scalar is
+ passed it is assumed that all elements of \code{A0} are equal to the
+ scalar. }
+
+ \item{b0}{ A vector containing the prior mean of each of the
+ discrimination parameters \eqn{\beta}{beta}. Should have as many
+ elements as items / roll calls. Order corresponds to the columns of
+ \code{datamatrix}. If a scalar is passed it is assumed that all
+ elements of \code{b0} are equal to the scalar. }
+
+ \item{B0}{ A vector containing the prior precision (inverse
+ variance) of each of the discrimination parameters
+ \eqn{\beta}{beta}. Should have as many elements as items / roll
+ calls. Order corresponds to the columns of \code{datamatrix}. If a
+ scalar is passed it is assumed that all elements of \code{B0} are
+ equal to the scalar. }
+
+ \item{c0}{\eqn{c_{0/2}}{c0/2} is the shape parameter for the inverse
+ Gamma prior on \eqn{\tau^2}{tau^2} (the variance of the random
+ walk increments). The amount of information in the inverse Gamma
+ prior is something like that from \eqn{c_0}{c0}
+ pseudo-observations. \code{c0} can be either a vector with an
+ element for each subject or a scalar. If \code{c0}
+ is negative then \eqn{\tau^2}{tau^2} is not estimated-- the values
+ in \code{tau2.start} are used throughout the sampling.}
+
+ \item{d0}{ \eqn{d_{0/2}}{d0/2} is the scale parameter for the inverse
+ Gamma prior on \eqn{\tau^2}{tau^2} (the variance of the random
+ walk increments). In constructing the inverse Gamma prior,
+ \eqn{d_0}{d0} acts like the sum of squared errors from the
+ \eqn{c_0}{c0} pseudo-observations. \code{d0} can be either a
+ vector with an element for each subject or a scalar. If \code{d0}
+ is negative then \eqn{\tau^2}{tau^2} is not estimated-- the values
+ in \code{tau2.start} are used throughout the sampling. }
+
+ \item{e0}{ A vector containing the prior mean of the initial ability
+ parameter / ideal point for each subject. Should have as many
+ elements as subjects. Order corresponds to the rows of \code{datamatrix}. If a
+ scalar is passed it is assumed that all elements of \code{e0} are
+ equal to the scalar.}
+
+ \item{E0}{ A vector containing the prior variance of the initial
+ ability parameter / ideal point for each subject. Should have as
+ many elements as subjects. Order corresponds to the rows of
+ \code{datamatrix}. If a scalar is passed it is assumed that all
+ elements of \code{E0} are equal to the scalar.}
+
+ \item{store.ability}{ A switch that determines whether or not to
+ store the ability parameters for posterior analysis.
+ \emph{NOTE}: In situations with many individuals storing the ability
+ parameters takes an enormous amount of memory, so
+ \code{store.ability} should only be \code{TRUE} if the chain is thinned
+ heavily, or for applications with a small number of individuals.
+ By default, the item parameters are stored. }
+
+ \item{store.item}{A switch that determines whether or not to
+ store the item parameters for posterior analysis.
+ \emph{NOTE}: In situations with many items storing the item
+ parameters takes an enormous amount of memory, so
+ \code{store.item} should only be \code{FALSE} if the chain is thinned
+ heavily, or for applications with a small number of items.
+ By default, the item parameters are not stored.}
+
+ \item{\dots}{further arguments to be passed }
+
+}
+
+\details{
+ \code{MCMCdynamicIRT1d} simulates from the posterior distribution using
+ the algorithm of Martin and Quinn (2002). The simulation proper is done in
+ compiled C++ code to maximize efficiency. Please consult the
+ coda documentation for a comprehensive list of functions that
+ can be used to analyze the posterior sample.
+
+ The model takes the following form. We assume that each subject has
+ an subject ability (ideal point) denoted
+ \eqn{\theta_{j,t}}{theta_jt} (where \eqn{j}{j} indexes subjects and
+ \eqn{t}{t} indexes time periods) and that each item has a difficulty
+ parameter \eqn{\alpha_i}{alpha_i} and discrimination parameter
+ \eqn{\beta_i}{beta_i}. The observed choice by subject \eqn{j}{j} on
+ item \eqn{i}{i} is the observed data matrix which is \eqn{(I \times
+ J)}{(I * J)}. We assume that the choice is dictated by an unobserved
+ utility: \deqn{z_{i,j,t} = -\alpha_i + \beta_i \theta_{j,t} +
+ \varepsilon_{i,j,t}}{z_ijt = -alpha_i + beta_i*theta_jt + epsilon_ijt}
+ Where the disturbances are assumed to be distributed standard Normal. The
+ parameters of interest are the subject abilities (ideal points) and
+ the item parameters.
+
+ We assume the following priors. For the subject abilities (ideal points):
+ \deqn{\theta_{j,t} \sim \mathcal{N}(\theta_{j,t-1}, \tau^2_j)}{theta_jt ~ N(theta_j(t-1), tau^2)} with \deqn{\theta_{j,0} \sim \mathcal{N}(e0, E0)}{theta_j0 ~ N(e0, E0).}
+
+ The evolution variance has the following prior: \deqn{\tau^2_j \sim \mathcal{IG}(c0/2, d0/2)}{tau^2_j ~ IG(c0/2, d0/2).}
+
+ For the item parameters in the standard model, the prior is:
+ \deqn{\alpha_i \sim \mathcal{N}(a0, A0^{-1})}{alpha_i ~ N(a0, A0^{-1})} and
+ \deqn{\beta_i \sim \mathcal{N}(b0, B0^{-1})}{beta_i ~ N(b0, B0^{-1}).}
+
+ The model is identified by the proper priors on the item parameters
+ and constraints placed on the ability parameters.
+
+ As is the case with all measurement models, make sure that you have plenty
+ of free memory, especially when storing the item parameters.
+}
+
+\value{ An mcmc object that contains the posterior sample. This object
+can be summarized by functions provided by the coda package. }
+
+\references{ Andrew D. Martin and Kevin M. Quinn. 2002. "Dynamic Ideal
+Point Estimation via Markov Chain Monte Carlo for the U.S. Supreme
+Court, 1953-1999." \emph{Political Analysis.} 10: 134-153.}
+
+\author{Kevin M. Quinn }
+
+\seealso{\code{\link[coda]{plot.mcmc}},\code{\link[coda]{summary.mcmc}},
+\code{\link[MCMCpack]{MCMCirt1d}} }
+
+
+\examples{
+ \dontrun{
+ data(Rehnquist)
+
+ ## assign starting values
+ theta.start <- rep(0, 9)
+ theta.start[2] <- -3 ## Stevens
+ theta.start[7] <- 2 ## Thomas
+
+ out <- MCMCdynamicIRT1d(t(Rehnquist[,1:9]),
+ item.time.map=Rehnquist$time,
+ theta.start=theta.start,
+ mcmc=50000, burnin=20000, thin=5,
+ verbose=500, tau2.start=rep(0.1, 9),
+ e0=0, E0=1,
+ a0=0, A0=1,
+ b0=0, B0=1, c0=-1, d0=-1,
+ store.item=FALSE,
+ theta.constraints=list(Stevens="-", Thomas="+"))
+
+ summary(out)
+ }
+}
+% Add one or more standard keywords, see file 'KEYWORDS' in the
+% R documentation directory.
+\keyword{ models }
+
diff --git a/man/MCMCirtKd.Rd b/man/MCMCirtKd.Rd
index 2325fce..149e388 100644
--- a/man/MCMCirtKd.Rd
+++ b/man/MCMCirtKd.Rd
@@ -168,7 +168,7 @@ MCMCirtKd(datamatrix, dimensions, item.constraints=list(),
for each item parameter.
The model is identified by the constraints on the item parameters
- (see Jackman 2001). The user cannot place constraints on the subect
+ (see Jackman 2001). The user cannot place constraints on the subject
abilities. This identification scheme differs from that in
\code{MCMCirt1d}, which uses constraints on
the subject abilities to identify the model.
diff --git a/man/Rehnquist.Rd b/man/Rehnquist.Rd
new file mode 100644
index 0000000..d147dd8
--- /dev/null
+++ b/man/Rehnquist.Rd
@@ -0,0 +1,32 @@
+\name{Rehnquist}
+\alias{Rehnquist}
+\title{
+ U.S. Supreme Court Vote Matrix, Rehnquist Court (1994-2004)
+}
+\description{
+ This dataframe contains a matrix of votes cast by U.S. Supreme
+ Court justices by all cases in the 1994-2004 terms.
+}
+\usage{
+data(SupremeCourt)
+}
+\format{
+ The dataframe has contains data for justices Rehnquist, Stevens,
+ O'Connor, Scalia, Kennedy, Souter, Thomas, Ginsburg, and Breyer
+ for the 1994-2004 terms of the U.S. Supreme Court. The dataframe
+ also contains the term of the case, and a time variable that counts
+ from term 1 to 11. The votes are coded liberal (1)
+ and conservative (0) using the protocol of Spaeth (2003).
+ The unit of analysis
+ is the case citation (ANALU=0). We are concerned with formally
+ decided cases issued with written opinions, after full oral
+ argument and cases decided by an equally divided vote
+ (DECTYPE=1,5,6,7).
+}
+
+\source{
+ Harold J. Spaeth. 2005. \emph{Original United States Supreme Court Database:
+ 1953-2004 Terms.} \url{http://www.as.uky.edu/polisci/ulmerproject/sctdata.htm}.
+}
+
+\keyword{datasets}
diff --git a/src/MCMCdynamicIRT1d-b.cc b/src/MCMCdynamicIRT1d-b.cc
new file mode 100644
index 0000000..48422e6
--- /dev/null
+++ b/src/MCMCdynamicIRT1d-b.cc
@@ -0,0 +1,714 @@
+//////////////////////////////////////////////////////////////////////////
+// MCMCdynamicIRT1d.cc is C++ code to estimate a dynamic 1d IRT model
+//
+// Kevin Quinn
+// 1/29/2008
+//
+// 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 MCMCDYNAMICIRT1DB_CC
+#define MCMCDYNAMICIRT1DB_CC
+
+#include<vector>
+#include<algorithm>
+
+#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 "mersenne.h"
+#include "lecuyer.h"
+
+#include <R.h> // needed to use Rprintf()
+#include <R_ext/Utils.h> // needed to allow user interrupts
+
+using namespace std;
+using namespace scythe;
+
+
+// used to access 1d arrays holding R matrices like a 2d array
+#define M(ROW,COL,NROWS) (COL*NROWS+ROW)
+
+
+// MCMCdynamicIRT1d implementation
+template <typename RNGTYPE>
+void MCMCdynamicIRT1d_b_impl(rng<RNGTYPE>& stream,
+ double* thetadraws, const int* nrowthetadraws,
+ const int* ncolthetadraws,
+ double* alphadraws, const int* nrowalphadraws,
+ const int* ncolalphadraws,
+ double* betadraws, const int* nrowbetadraws,
+ const int* ncolbetadraws,
+ double* tau2draws, const int* nrowtau2draws,
+ const int* ncoltau2draws,
+ const int* nsubj, const int* nitems,
+ const int* ntime,
+ const int* Ydata, const int* nrowYdata,
+ const int* ncolYdata,
+ const int* ITdata, const int* lengthITdata,
+ const int* burnin, const int* mcmc, const int* thin,
+ const int* verbose,
+ const double* thetadata,
+ const int* lengththeta,
+ const int* thetainfodata,
+ const int* nrowthetainfo,
+ const int* ncolthetainfo,
+ double* alphadata,
+ const int* lengthalpha,
+ double* betadata,
+ const int* lengthbeta,
+ double* tau2data,
+ const int* lengthtau2,
+ const double* c0,
+ const int* lengthc0,
+ const double* d0,
+ const int* lengthd0,
+ const double* a0,
+ const double* A0,
+ const double* b0,
+ const double* B0,
+ const double* e0,
+ const double* E0inv,
+ const double* thetaeqdata,
+ const int* nrowthetaeq,
+ const int* ncolthetaeq,
+ const double* thetaineqdata,
+ const int* nrowthetaineq,
+ const int* ncolthetaineq,
+ const int* storeitem,
+ const int* storeability){
+
+
+
+ const int tot_iter = *burnin + *mcmc;
+ const int nsamp = *mcmc / *thin;
+ //sparse matrix of latent outcome variables
+ double Z[*nsubj][*nitems];
+ for (int s=0; s<*nsubj; ++s){
+ for (int i=0; i<*nitems; ++i){
+ Z[s][i] = -999;
+ }
+ }
+ for (int j=0; j<*nrowYdata; ++j){
+ const int s = Ydata[M(j, 0, *nrowYdata)]; // subject id
+ const int i = Ydata[M(j, 1, *nrowYdata)]; // item id
+ const int y = Ydata[M(j, 2, *nrowYdata)]; // y value
+ Z[s][i] = 0.0;
+ }
+
+
+ // stuff to make working with theta easy
+ // theta[s][t] gives the tth theta for subject s
+ // the actual time period corresponding to theta[s][t] is theta_offset[s]+t
+ vector< vector<double> > theta;
+ vector<int> theta_offset;
+ int count = 0;
+ for (int s=0; s<*nsubj; ++s){
+ vector<double> holder;
+ int ntime_s = thetainfodata[M(s, 1, *nrowthetainfo)];
+ theta_offset.push_back(thetainfodata[M(s, 2, *nrowthetainfo)]);
+ for (int t=0; t<ntime_s; ++t){
+ holder.push_back(thetadata[count]);
+ ++count;
+ }
+ theta.push_back(holder);
+ }
+
+
+
+ // define constants used in the sampling
+ const double A0a0 = *A0 * *a0;
+ const double B0b0 = *B0 * *b0;
+
+
+ // set up mappings
+ //
+ // IS gives the mapping from items to subjects
+ // IS[i] provides a vector of integers corresponding to the subjects
+ // who voted on item i.
+ // IS[i][s] gets the subject index of the sth subject to vote on item i
+ vector< vector<int> > IS;
+ for (int i=0; i<*nitems; ++i){
+ vector<int> subjholder;
+ for (int j=0; j<*nrowYdata; ++j){
+ if (Ydata[M(j,1,*nrowYdata)] == i){
+ subjholder.push_back(Ydata[M(j,0,*nrowYdata)]);
+ }
+ }
+ sort(subjholder.begin(), subjholder.end());
+ IS.push_back(subjholder);
+ }
+
+ // SI gives the mapping from subjects to items
+ // SI[s] provides a vector of integers corresponding to the items
+ // voted on by subject s.
+ // SI[s][i] gets the item index of the ith item voted on by subject s
+ vector< vector<int> > SI;
+ for (int i=0; i<*nsubj; ++i){
+ vector<int> itemholder;
+ for (int j=0; j<*nrowYdata; ++j){
+ if (Ydata[M(j,0,*nrowYdata)] == i){
+ itemholder.push_back(Ydata[M(j,1,*nrowYdata)]);
+ }
+ }
+ sort(itemholder.begin(), itemholder.end());
+ SI.push_back(itemholder);
+ }
+
+ // TI gives the mapping from times to items
+ // TI[t] provides a vector of integers corresponding to the items
+ // voted on in time period t.
+ // TI[t][i] gets the item index of the ith item voted on in time period t
+ vector< vector<int> > TI;
+ for (int t=0; t<*ntime; ++t){
+ vector<int> itemholder;
+ for (int i=0; i<*lengthITdata; ++i){
+ if (ITdata[i] == t){
+ itemholder.push_back(i);
+ }
+ }
+ sort(itemholder.begin(), itemholder.end());
+ TI.push_back(itemholder);
+ }
+
+ // IT gives the mapping from items to times and is just fancy
+ // way of holding the stuff in *ITdata
+ vector<int> IT;
+ for (int i=0; i<*nitems; ++i){
+ IT.push_back(ITdata[i]);
+ }
+
+
+ // ST gives the mapping from subjects to times
+ // ST[s] provides a vector of integers corresponding to the times
+ // in which items were voted on by subject s.
+ // ST[s][t] gets the time index of the tth time period served by subject s
+ /*
+ vector <vector <int> > ST;
+ for (int s=0; s<*nsubj; ++s){
+ vector <int> timeholder;
+ for (int iind=0; iind<SI[s].size(); ++iind){
+ const int i = SI[s][iind];
+ const int t = IT[i];
+ if (timeholder.empty()){
+ timeholder.push_back(t);
+ }
+ vector<int>::iterator myiter;
+ myiter = find(timeholder.begin(), timeholder.end(), t);
+ if (myiter == timeholder.end()){ // t not currently in timeholder
+ timeholder.push_back(t);
+ }
+ }
+ sort(timeholder.begin(), timeholder.end());
+ ST.push_back(timeholder);
+ }
+ */
+
+ // STI gives the list of item indices of items voted on by by subject s
+ // in the tth time period served by s
+ vector< vector < vector <int> > > STI;
+ for (int s=0; s<*nsubj; ++s){
+ vector < vector <int> > timeitemholder;
+ for (int tt=0; tt<theta[s].size(); ++tt){
+ const int t = tt + theta_offset[s];
+ vector <int> itemholder;
+ for (int ii=0; ii<TI[t].size(); ++ii){
+ const int i = TI[t][ii];
+ if (Z[s][i] != -999){
+ itemholder.push_back(i);
+ }
+ }
+ timeitemholder.push_back(itemholder);
+ }
+ STI.push_back(timeitemholder);
+ }
+
+
+
+
+
+ // MCMC iterations start here
+ int sampcount = 0;
+ for (int iter=0; iter < tot_iter; ++iter){
+
+
+
+ // sample latent Z
+ for (int j=0; j<*nrowYdata; ++j){
+ const int s = Ydata[M(j, 0, *nrowYdata)]; // subject id
+ const int i = Ydata[M(j, 1, *nrowYdata)]; // item id
+ const int y = Ydata[M(j, 2, *nrowYdata)]; // y value
+ const int t = IT[i]; // time period
+ const double mu = -alphadata[i] +
+ betadata[i] * theta[s][t-theta_offset[s]];
+ if (y == 1.0){
+ Z[s][i] = stream.rtbnorm_combo(mu, 1.0, 0);
+ }
+ if (y == 0.0){
+ Z[s][i] = stream.rtanorm_combo(mu, 1.0, 0);
+ }
+ }
+
+
+
+
+
+ // sample item parameters (alpha, beta)
+ for (int i=0; i<*nitems; ++i){
+ const int nsubj_i = IS[i].size();
+ double sumthetasq = 0.0;
+ double sumtheta = 0.0;
+ double sumz = 0.0;
+ double sumthetaz = 0.0;
+ const int t = IT[i];
+
+ for (int ss=0; ss<nsubj_i; ++ss){
+ const int s = IS[i][ss];
+ const double theta_st = theta[s][t - theta_offset[s]];
+ sumthetasq += std::pow(theta_st, 2.0);
+ sumtheta += -1*theta_st;
+ sumz += -1*Z[s][i];
+ sumthetaz += Z[s][i] * theta_st;
+ } // end ss loop
+
+ const double sumthetaallsq = std::pow(sumtheta, 2.0);
+ const double invdet = 1.0 / ((nsubj_i + *A0) *
+ (sumthetasq + *B0) -
+ sumthetaallsq);
+ const double v11star = invdet * (sumthetasq + *B0);
+ const double v12star = invdet * (-1.0 * sumtheta);
+ const double v22star = invdet * (nsubj_i + *B0);
+ const double s1star = std::sqrt(v11star);
+ const double s2star = std::sqrt(v22star);
+ const double rho = v12star / (s1star * s2star);
+ const double holder1 = sumz + A0a0;
+ const double holder2 = sumthetaz + B0b0;
+ const double m1star = v11star * holder1 + v12star * holder2;
+ const double m2star = v12star * holder1 + v22star * holder2;
+ // (alpha[i], beta[i]) ~
+ // N((m1star, m2star), c(v11star, v12star, v22star) )
+ alphadata[i] = 1000;
+ betadata[i] = 100;
+ while (std::fabs(alphadata[i]) > 4 ||
+ std::fabs(betadata[i]) > 2 ||
+ std::fabs(alphadata[i] / betadata[i]) > 2){
+ alphadata[i] = stream.rnorm(m1star, s1star);
+ if (iter < 20 ){
+ alphadata[i] = 0.2*((m1star > 0) - (m1star < 0));
+ }
+ const double cond_mean = m2star - m1star * (v12star / v11star) +
+ alphadata[i] * (v12star / v11star);
+ const double cond_sd = std::sqrt(v22star * ( 1 - std::pow(rho, 2.0)));
+ betadata[i] = stream.rnorm(cond_mean, cond_sd);
+ if (iter < 20 ){
+ betadata[i] = 1.0*((cond_mean > 0) - (cond_mean < 0));
+ }
+ }
+ } // end i loop
+
+
+
+
+
+ // sample subject parameters (theta, tau2)
+ for (int s=0; s<*nsubj; ++s){
+ const int ntime_s = theta[s].size();
+ if (thetaeqdata[s] != -999){
+ for (int t=0; t<ntime_s; ++t){
+ theta[s][t] = thetaeqdata[s];
+ }
+ }
+ else{
+
+ // time period 0
+ vector<double> beta_s0;
+ beta_s0.reserve(STI[s][0].size());
+ vector<double> zalpha_s0;
+ zalpha_s0.reserve(STI[s][0].size());
+ for (int ii=0; ii<STI[s][0].size(); ++ii){
+ const int i = STI[s][0][ii];
+ beta_s0.push_back(betadata[i]);
+ zalpha_s0.push_back(Z[s][i] + alphadata[i]);
+ }
+ vector<double> a;
+ a.reserve(ntime_s);
+ a.push_back(e0[s]);
+ vector<double> R;
+ R.reserve(ntime_s);
+ R.push_back(E0inv[s] + tau2data[s]);
+ vector <double> f_0 = beta_s0;
+ for (int i=0; i<f_0.size(); ++i){
+ f_0[i] = f_0[i] * a[0];
+ }
+ Matrix <> Q_0(beta_s0.size(), beta_s0.size());
+ for (int i=0; i<beta_s0.size(); ++i){
+ for (int j=i; j<beta_s0.size(); ++j){
+ if (i!=j){
+ Q_0(i,j) = Q_0(j,i) = beta_s0[i] * beta_s0[j] * R[0];
+ }
+ else{
+ Q_0(i,j) = beta_s0[i] * beta_s0[j] * R[0] + 1.0;
+ }
+ }
+ }
+ vector <double> e_0 = zalpha_s0;
+ for (int i=0; i<e_0.size(); ++i){
+ e_0[i] = e_0[i] - f_0[i];
+ }
+ const Matrix <> Q_0_inv = invpd(Q_0);
+
+ vector <double> A_0;
+ A_0.reserve(beta_s0.size());
+ for (int i=0; i<beta_s0.size(); ++i){
+ double sumholder = 0.0;
+ for (int j=0; j<beta_s0.size(); ++j){
+ sumholder += beta_s0[j] * Q_0_inv(j,i);
+ }
+ A_0.push_back(sumholder * R[0]);
+ }
+ vector<double> m;
+ m.reserve(ntime_s);
+ double mhold = a[0];
+ for (int i=0; i<A_0.size(); ++i){
+ mhold += A_0[i] * e_0[i];
+ }
+ m.push_back(mhold);
+
+
+ vector<double> C;
+ C.reserve(ntime_s);
+ double Chold = 0.0;
+
+ for (int i=0; i<A_0.size(); ++i){
+ double hold2 = 0.0;
+ for (int j=0; j<A_0.size(); ++j){
+ //cout << "i = " << i << " j = " << j << endl;
+ hold2 += Q_0(i,j) * A_0[j];
+ }
+
+ Chold += hold2 * A_0[i];
+ }
+ C.push_back(R[0] - Chold);
+
+
+
+ // THE FORWARD RECURSION
+ // time periods 1 to T
+ for (int tt=1; tt<ntime_s; ++tt){
+ if (STI[s][tt].size() == 0){
+ a.push_back(m[tt-1]);
+ R.push_back(C[tt-1] + tau2data[s]);
+ m.push_back(a[tt]);
+ C.push_back(R[tt]);
+ }
+ else{
+ vector<double> beta_s;
+ beta_s.reserve(STI[s][tt].size());
+ vector<double> zalpha_s;
+ zalpha_s.reserve(STI[s][tt].size());
+ for (int ii=0; ii<STI[s][tt].size(); ++ii){
+ const int i = STI[s][tt][ii];
+ beta_s.push_back(betadata[i]);
+ zalpha_s.push_back(Z[s][i] + alphadata[i]);
+ }
+ // a
+ a.push_back(m[tt-1]);
+ // R
+ R.push_back(C[tt-1] + tau2data[s]);
+ vector <double> f = beta_s;
+ for (int i=0; i<f.size(); ++i){
+ f[i] = f[i] * a[tt];
+ }
+ Matrix <> Q(beta_s.size(), beta_s.size());
+ for (int i=0; i<beta_s.size(); ++i){
+ for (int j=i; j<beta_s.size(); ++j){
+ if (i!=j){
+ Q(i,j) = Q(j,i) = beta_s[i] * beta_s[j] * R[tt];
+ }
+ else{
+ Q(i,j) = beta_s[i] * beta_s[j] * R[tt] + 1.0;
+ }
+ }
+ }
+ vector <double> e = zalpha_s;
+ for (int i=0; i<e.size(); ++i){
+ e[i] = e[i] - f[i];
+ }
+ const Matrix <> Q_inv = invpd(Q);
+
+ vector <double> A;
+ A.reserve(beta_s.size());
+ for (int i=0; i<beta_s.size(); ++i){
+ double sumholder = 0.0;
+ for (int j=0; j<beta_s.size(); ++j){
+ sumholder += beta_s[j] * Q_inv(j,i);
+ }
+ A.push_back(sumholder * R[tt]);
+ }
+ mhold = a[tt];
+ for (int i=0; i<A.size(); ++i){
+ mhold += A[i] * e[i];
+ }
+ m.push_back(mhold);
+
+ Chold = 0.0;
+ for (int i=0; i<A.size(); ++i){
+ double hold2 = 0.0;
+ for (int j=0; j<A.size(); ++j){
+ hold2 += Q(i,j) * A[j];
+ }
+ Chold += hold2 * A[i];
+ }
+ C.push_back(R[tt] - Chold);
+ }
+ } // end tt loop
+
+ if (thetaineqdata[s] == 0){
+ theta[s][ntime_s-1] = stream.rnorm(m[ntime_s-1],
+ std::sqrt(C[ntime_s-1]));
+ }
+ else if (thetaineqdata[s] > 0){
+ theta[s][ntime_s-1] = stream.rtbnorm_combo(m[ntime_s-1],
+ C[ntime_s-1], 0.0);
+ }
+ else{
+ theta[s][ntime_s-1] = stream.rtanorm_combo(m[ntime_s-1],
+ C[ntime_s-1], 0.0);
+ }
+
+
+
+
+ for (int tt=(ntime_s-2); tt>=0; --tt){
+ const double B = C[tt] * (1.0 / R[tt+1]);
+ const double h = m[tt] + B * (theta[s][tt+1] - a[tt+1]);
+ const double H = C[tt] - B * R[tt+1] * B;
+ if (thetaineqdata[s] == 0){
+ theta[s][tt] = stream.rnorm(h, std::sqrt(H));
+ }
+ else if (thetaineqdata[s] > 0){
+ theta[s][tt] = stream.rtbnorm_combo(h, H, 0.0);
+ }
+ else {
+ theta[s][tt] = stream.rtanorm_combo(h, H, 0.0);
+ }
+ } // end backwards tt loop
+
+ } // end theteqdataa[s] == -999
+
+
+
+ // sample smoothing parameters (tau2)
+ if (c0[s] > 0 && d0[s] > 0){
+ double SSE = 0.0;
+ for (int t=1; t<ntime_s; ++t){
+ SSE += std::pow((theta[s][t] - theta[s][t-1]), 2.0);
+ }
+ const double param1 = (c0[s] + ntime_s - 1) * 0.5;
+ const double param2 = (d0[s] + SSE) * 0.5;
+ tau2data[s] = stream.rigamma(param1, param2);
+ }
+
+
+
+
+ } // end s loop (end of theta tau2 loop)
+
+
+
+
+
+
+
+
+
+
+
+
+
+ // store draws
+ if (iter >= *burnin && (iter % *thin == 0)) {
+ if (*storeability == 1){
+ int internalcount = 0;
+ for (int s=0; s<*nsubj; ++s){
+ for (int t=0; t<theta[s].size(); ++t){
+ thetadraws[M(sampcount, internalcount, *nrowthetadraws)] =
+ theta[s][t];
+ ++internalcount;
+ }
+ }
+ }
+ if (*storeitem == 1){
+ for (int i=0; i<*nitems; ++i){
+ alphadraws[M(sampcount, i, *nrowalphadraws)] = alphadata[i];
+ betadraws[M(sampcount, i, *nrowbetadraws)] = betadata[i];
+ }
+ }
+ for (int s=0; s<*nsubj; ++s){
+ tau2draws[M(sampcount, s, *nrowtau2draws)] = tau2data[s];
+ }
+ ++sampcount;
+ }
+
+
+ // print output to stdout
+ // print output to stdout
+ if(*verbose > 0 && iter % *verbose == 0){
+ Rprintf("\n\nMCMCdynamicIRT1d iteration %i of %i \n", (iter+1), tot_iter);
+ /*
+ for (int t=0; t<theta[1].size(); ++t){
+ Rprintf("%f ", theta[1][t]);
+ }
+ Rprintf("\n");
+ for (int t=0; t<theta[7].size(); ++t){
+ Rprintf("%f ", theta[7][t]);
+ }
+ Rprintf("\n");
+ for (int t=0; t<theta[5].size(); ++t){
+ Rprintf("%f ", theta[5][t]);
+ }
+ Rprintf("\n");
+ for (int t=0; t<theta[8].size(); ++t){
+ Rprintf("%f ", theta[8][t]);
+ }
+ Rprintf("\n");
+ for (int t=0; t<theta[2].size(); ++t){
+ Rprintf("%f ", theta[2][t]);
+ }
+ Rprintf("\n");
+ for (int t=0; t<theta[4].size(); ++t){
+ Rprintf("%f ", theta[4][t]);
+ }
+ Rprintf("\n");
+ for (int t=0; t<theta[0].size(); ++t){
+ Rprintf("%f ", theta[0][t]);
+ }
+ Rprintf("\n");
+ for (int t=0; t<theta[3].size(); ++t){
+ Rprintf("%f ", theta[3][t]);
+ }
+ Rprintf("\n");
+ for (int t=0; t<theta[6].size(); ++t){
+ Rprintf("%f ", theta[6][t]);
+ }
+ Rprintf("\n");
+ Rprintf("\n\n");
+ Rprintf("alpha[0] = %f \n", alphadata[0]);
+ Rprintf("beta[0] = %f \n", betadata[0]);
+ Rprintf("alpha[0]/beta[0] = %f \n", alphadata[0]/betadata[0]);
+ Rprintf("tau2[1] = %f \n", tau2data[1]);
+ Rprintf("tau2[7] = %f \n", tau2data[7]);
+ Rprintf("tau2[5] = %f \n", tau2data[5]);
+ Rprintf("tau2[8] = %f \n", tau2data[8]);
+ Rprintf("tau2[2] = %f \n", tau2data[2]);
+ Rprintf("tau2[4] = %f \n", tau2data[4]);
+ Rprintf("tau2[0] = %f \n", tau2data[0]);
+ Rprintf("tau2[3] = %f \n", tau2data[3]);
+ Rprintf("tau2[6] = %f \n", tau2data[6]);
+ Rprintf("\n\n");
+ */
+ }
+
+ R_CheckUserInterrupt(); // allow user interrupts
+
+
+
+ } // end iter loop
+
+
+
+} // end MCMCdynamicIRT1d_impl function
+
+
+
+
+
+
+
+
+extern "C" {
+
+ void MCMCdynamicIRT1d_b(double* thetadraws, const int* nrowthetadraws,
+ const int* ncolthetadraws,
+ double* alphadraws, const int* nrowalphadraws,
+ const int* ncolalphadraws,
+ double* betadraws, const int* nrowbetadraws,
+ const int* ncolbetadraws,
+ double* tau2draws, const int* nrowtau2draws,
+ const int* ncoltau2draws,
+ const int* nsubj, const int* nitems,
+ const int* ntime,
+ const int* Ydata, const int* nrowYdata,
+ const int* ncolYdata,
+ const int* ITdata, const int* lengthITdata,
+ const int* burnin, const int* mcmc, const int* thin,
+ const int* uselecuyer, const int* seedarray,
+ const int* lecuyerstream,
+ const int* verbose,
+ const double* thetastartdata,
+ const int* lengththetastart,
+ const int* thetainfodata,
+ const int* nrowthetainfo,
+ const int* ncolthetainfo,
+ double* alphastartdata,
+ const int* lengthalphastart,
+ double* betastartdata,
+ const int* lengthbetastart,
+ double* tau2startdata,
+ const int* lengthtau2start,
+ const double* c0,
+ const int* lengthc0,
+ const double* d0,
+ const int* lengthd0,
+ const double* a0,
+ double* A0,
+ const double* b0,
+ double* B0,
+ const double* e0,
+ const double* E0inv,
+ const double* thetaeqdata,
+ const int* nrowthetaeq,
+ const int* ncolthetaeq,
+ const double* thetaineqdata,
+ const int* nrowthetaineq,
+ const int* ncolthetaineq,
+ const int* storeitem,
+ const int* storeability
+ ){
+
+ *A0 = 0.0;
+ *B0 = 0.0;
+
+ MCMCPACK_PASSRNG2MODEL(MCMCdynamicIRT1d_b_impl, thetadraws, nrowthetadraws,
+ ncolthetadraws, alphadraws, nrowalphadraws,
+ ncolalphadraws, betadraws, nrowbetadraws,
+ ncolbetadraws, tau2draws, nrowtau2draws,
+ ncoltau2draws, nsubj, nitems, ntime,
+ Ydata, nrowYdata, ncolYdata,
+ ITdata, lengthITdata, burnin, mcmc, thin,
+ verbose, thetastartdata, lengththetastart,
+ thetainfodata, nrowthetainfo, ncolthetainfo,
+ alphastartdata, lengthalphastart, betastartdata,
+ lengthbetastart, tau2startdata, lengthtau2start,
+ c0, lengthc0, d0, lengthd0,
+ a0, A0, b0, B0, e0, E0inv,
+ thetaeqdata, nrowthetaeq,
+ ncolthetaeq, thetaineqdata, nrowthetaineq,
+ ncolthetaineq, storeitem, storeability );
+
+
+ } // end MCMCdynamicIRT1d
+
+} // end extern "C"
+
+
+#endif
diff --git a/src/MCMCdynamicIRT1d.cc b/src/MCMCdynamicIRT1d.cc
new file mode 100644
index 0000000..8768205
--- /dev/null
+++ b/src/MCMCdynamicIRT1d.cc
@@ -0,0 +1,699 @@
+//////////////////////////////////////////////////////////////////////////
+// MCMCdynamicIRT1d.cc is C++ code to estimate a dynamic 1d IRT model
+//
+// Kevin Quinn
+// 1/29/2008
+//
+// 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 MCMCDYNAMICIRT1D_CC
+#define MCMCDYNAMICIRT1D_CC
+
+#include<vector>
+#include<algorithm>
+
+#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 "mersenne.h"
+#include "lecuyer.h"
+
+#include <R.h> // needed to use Rprintf()
+#include <R_ext/Utils.h> // needed to allow user interrupts
+
+using namespace std;
+using namespace scythe;
+
+
+// used to access 1d arrays holding R matrices like a 2d array
+#define M(ROW,COL,NROWS) (COL*NROWS+ROW)
+
+
+// MCMCdynamicIRT1d implementation
+template <typename RNGTYPE>
+void MCMCdynamicIRT1d_impl(rng<RNGTYPE>& stream,
+ double* thetadraws, const int* nrowthetadraws,
+ const int* ncolthetadraws,
+ double* alphadraws, const int* nrowalphadraws,
+ const int* ncolalphadraws,
+ double* betadraws, const int* nrowbetadraws,
+ const int* ncolbetadraws,
+ double* tau2draws, const int* nrowtau2draws,
+ const int* ncoltau2draws,
+ const int* nsubj, const int* nitems,
+ const int* ntime,
+ const int* Ydata, const int* nrowYdata,
+ const int* ncolYdata,
+ const int* ITdata, const int* lengthITdata,
+ const int* burnin, const int* mcmc, const int* thin,
+ const int* verbose,
+ const double* thetadata,
+ const int* lengththeta,
+ const int* thetainfodata,
+ const int* nrowthetainfo,
+ const int* ncolthetainfo,
+ double* alphadata,
+ const int* lengthalpha,
+ double* betadata,
+ const int* lengthbeta,
+ double* tau2data,
+ const int* lengthtau2,
+ const double* c0,
+ const int* lengthc0,
+ const double* d0,
+ const int* lengthd0,
+ const double* a0,
+ const double* A0,
+ const double* b0,
+ const double* B0,
+ const double* e0,
+ const double* E0inv,
+ const double* thetaeqdata,
+ const int* nrowthetaeq,
+ const int* ncolthetaeq,
+ const double* thetaineqdata,
+ const int* nrowthetaineq,
+ const int* ncolthetaineq,
+ const int* storeitem,
+ const int* storeability){
+
+
+
+ const int tot_iter = *burnin + *mcmc;
+ const int nsamp = *mcmc / *thin;
+ //sparse matrix of latent outcome variables
+ double Z[*nsubj][*nitems];
+ for (int s=0; s<*nsubj; ++s){
+ for (int i=0; i<*nitems; ++i){
+ Z[s][i] = -999;
+ }
+ }
+ for (int j=0; j<*nrowYdata; ++j){
+ const int s = Ydata[M(j, 0, *nrowYdata)]; // subject id
+ const int i = Ydata[M(j, 1, *nrowYdata)]; // item id
+ const int y = Ydata[M(j, 2, *nrowYdata)]; // y value
+ Z[s][i] = 0.0;
+ }
+
+
+ // stuff to make working with theta easy
+ // theta[s][t] gives the tth theta for subject s
+ // the actual time period corresponding to theta[s][t] is theta_offset[s]+t
+ vector< vector<double> > theta;
+ vector<int> theta_offset;
+ int count = 0;
+ for (int s=0; s<*nsubj; ++s){
+ vector<double> holder;
+ int ntime_s = thetainfodata[M(s, 1, *nrowthetainfo)];
+ theta_offset.push_back(thetainfodata[M(s, 2, *nrowthetainfo)]);
+ for (int t=0; t<ntime_s; ++t){
+ holder.push_back(thetadata[count]);
+ ++count;
+ }
+ theta.push_back(holder);
+ }
+
+
+
+ // define constants used in the sampling
+ const double A0a0 = *A0 * *a0;
+ const double B0b0 = *B0 * *b0;
+
+
+ // set up mappings
+ //
+ // IS gives the mapping from items to subjects
+ // IS[i] provides a vector of integers corresponding to the subjects
+ // who voted on item i.
+ // IS[i][s] gets the subject index of the sth subject to vote on item i
+ vector< vector<int> > IS;
+ for (int i=0; i<*nitems; ++i){
+ vector<int> subjholder;
+ for (int j=0; j<*nrowYdata; ++j){
+ if (Ydata[M(j,1,*nrowYdata)] == i){
+ subjholder.push_back(Ydata[M(j,0,*nrowYdata)]);
+ }
+ }
+ sort(subjholder.begin(), subjholder.end());
+ IS.push_back(subjholder);
+ }
+
+ // SI gives the mapping from subjects to items
+ // SI[s] provides a vector of integers corresponding to the items
+ // voted on by subject s.
+ // SI[s][i] gets the item index of the ith item voted on by subject s
+ vector< vector<int> > SI;
+ for (int i=0; i<*nsubj; ++i){
+ vector<int> itemholder;
+ for (int j=0; j<*nrowYdata; ++j){
+ if (Ydata[M(j,0,*nrowYdata)] == i){
+ itemholder.push_back(Ydata[M(j,1,*nrowYdata)]);
+ }
+ }
+ sort(itemholder.begin(), itemholder.end());
+ SI.push_back(itemholder);
+ }
+
+ // TI gives the mapping from times to items
+ // TI[t] provides a vector of integers corresponding to the items
+ // voted on in time period t.
+ // TI[t][i] gets the item index of the ith item voted on in time period t
+ vector< vector<int> > TI;
+ for (int t=0; t<*ntime; ++t){
+ vector<int> itemholder;
+ for (int i=0; i<*lengthITdata; ++i){
+ if (ITdata[i] == t){
+ itemholder.push_back(i);
+ }
+ }
+ sort(itemholder.begin(), itemholder.end());
+ TI.push_back(itemholder);
+ }
+
+ // IT gives the mapping from items to times and is just fancy
+ // way of holding the stuff in *ITdata
+ vector<int> IT;
+ for (int i=0; i<*nitems; ++i){
+ IT.push_back(ITdata[i]);
+ }
+
+
+ // ST gives the mapping from subjects to times
+ // ST[s] provides a vector of integers corresponding to the times
+ // in which items were voted on by subject s.
+ // ST[s][t] gets the time index of the tth time period served by subject s
+ /*
+ vector <vector <int> > ST;
+ for (int s=0; s<*nsubj; ++s){
+ vector <int> timeholder;
+ for (int iind=0; iind<SI[s].size(); ++iind){
+ const int i = SI[s][iind];
+ const int t = IT[i];
+ if (timeholder.empty()){
+ timeholder.push_back(t);
+ }
+ vector<int>::iterator myiter;
+ myiter = find(timeholder.begin(), timeholder.end(), t);
+ if (myiter == timeholder.end()){ // t not currently in timeholder
+ timeholder.push_back(t);
+ }
+ }
+ sort(timeholder.begin(), timeholder.end());
+ ST.push_back(timeholder);
+ }
+ */
+
+ // STI gives the list of item indices of items voted on by by subject s
+ // in the tth time period served by s
+ vector< vector < vector <int> > > STI;
+ for (int s=0; s<*nsubj; ++s){
+ vector < vector <int> > timeitemholder;
+ for (int tt=0; tt<theta[s].size(); ++tt){
+ const int t = tt + theta_offset[s];
+ vector <int> itemholder;
+ for (int ii=0; ii<TI[t].size(); ++ii){
+ const int i = TI[t][ii];
+ if (Z[s][i] != -999){
+ itemholder.push_back(i);
+ }
+ }
+ timeitemholder.push_back(itemholder);
+ }
+ STI.push_back(timeitemholder);
+ }
+
+
+
+
+
+ // MCMC iterations start here
+ int sampcount = 0;
+ for (int iter=0; iter < tot_iter; ++iter){
+
+
+
+ // sample latent Z
+ for (int j=0; j<*nrowYdata; ++j){
+ const int s = Ydata[M(j, 0, *nrowYdata)]; // subject id
+ const int i = Ydata[M(j, 1, *nrowYdata)]; // item id
+ const int y = Ydata[M(j, 2, *nrowYdata)]; // y value
+ const int t = IT[i]; // time period
+ const double mu = -alphadata[i] +
+ betadata[i] * theta[s][t-theta_offset[s]];
+ if (y == 1.0){
+ Z[s][i] = stream.rtbnorm_combo(mu, 1.0, 0);
+ }
+ if (y == 0.0){
+ Z[s][i] = stream.rtanorm_combo(mu, 1.0, 0);
+ }
+ }
+
+
+
+
+
+ // sample item parameters (alpha, beta)
+ for (int i=0; i<*nitems; ++i){
+ const int nsubj_i = IS[i].size();
+ double sumthetasq = 0.0;
+ double sumtheta = 0.0;
+ double sumz = 0.0;
+ double sumthetaz = 0.0;
+ const int t = IT[i];
+
+ for (int ss=0; ss<nsubj_i; ++ss){
+ const int s = IS[i][ss];
+ const double theta_st = theta[s][t - theta_offset[s]];
+ sumthetasq += std::pow(theta_st, 2.0);
+ sumtheta += -1*theta_st;
+ sumz += -1*Z[s][i];
+ sumthetaz += Z[s][i] * theta_st;
+ } // end ss loop
+
+ const double sumthetaallsq = std::pow(sumtheta, 2.0);
+ const double invdet = 1.0 / ((nsubj_i + *A0) *
+ (sumthetasq + *B0) -
+ sumthetaallsq);
+ const double v11star = invdet * (sumthetasq + *B0);
+ const double v12star = invdet * (-1.0 * sumtheta);
+ const double v22star = invdet * (nsubj_i + *B0);
+ const double s1star = std::sqrt(v11star);
+ const double s2star = std::sqrt(v22star);
+ const double rho = v12star / (s1star * s2star);
+ const double holder1 = sumz + A0a0;
+ const double holder2 = sumthetaz + B0b0;
+ const double m1star = v11star * holder1 + v12star * holder2;
+ const double m2star = v12star * holder1 + v22star * holder2;
+ // (alpha[i], beta[i]) ~
+ // N((m1star, m2star), c(v11star, v12star, v22star) )
+ alphadata[i] = stream.rnorm(m1star, s1star);
+ const double cond_mean = m2star - m1star * (v12star / v11star) +
+ alphadata[i] * (v12star / v11star);
+ const double cond_sd = std::sqrt(v22star * ( 1 - std::pow(rho, 2.0)));
+ betadata[i] = stream.rnorm(cond_mean, cond_sd);
+ } // end i loop
+
+
+
+
+
+ // sample subject parameters (theta, tau2)
+ for (int s=0; s<*nsubj; ++s){
+ const int ntime_s = theta[s].size();
+ if (thetaeqdata[s] != -999){
+ for (int t=0; t<ntime_s; ++t){
+ theta[s][t] = thetaeqdata[s];
+ }
+ }
+ else{
+
+ // time period 0
+ vector<double> beta_s0;
+ beta_s0.reserve(STI[s][0].size());
+ vector<double> zalpha_s0;
+ zalpha_s0.reserve(STI[s][0].size());
+ for (int ii=0; ii<STI[s][0].size(); ++ii){
+ const int i = STI[s][0][ii];
+ beta_s0.push_back(betadata[i]);
+ zalpha_s0.push_back(Z[s][i] + alphadata[i]);
+ }
+ vector<double> a;
+ a.reserve(ntime_s);
+ a.push_back(e0[s]);
+ vector<double> R;
+ R.reserve(ntime_s);
+ R.push_back(E0inv[s] + tau2data[s]);
+ vector <double> f_0 = beta_s0;
+ for (int i=0; i<f_0.size(); ++i){
+ f_0[i] = f_0[i] * a[0];
+ }
+ Matrix <> Q_0(beta_s0.size(), beta_s0.size());
+ for (int i=0; i<beta_s0.size(); ++i){
+ for (int j=i; j<beta_s0.size(); ++j){
+ if (i!=j){
+ Q_0(i,j) = Q_0(j,i) = beta_s0[i] * beta_s0[j] * R[0];
+ }
+ else{
+ Q_0(i,j) = beta_s0[i] * beta_s0[j] * R[0] + 1.0;
+ }
+ }
+ }
+ vector <double> e_0 = zalpha_s0;
+ for (int i=0; i<e_0.size(); ++i){
+ e_0[i] = e_0[i] - f_0[i];
+ }
+ const Matrix <> Q_0_inv = invpd(Q_0);
+
+ vector <double> A_0;
+ A_0.reserve(beta_s0.size());
+ for (int i=0; i<beta_s0.size(); ++i){
+ double sumholder = 0.0;
+ for (int j=0; j<beta_s0.size(); ++j){
+ sumholder += beta_s0[j] * Q_0_inv(j,i);
+ }
+ A_0.push_back(sumholder * R[0]);
+ }
+ vector<double> m;
+ m.reserve(ntime_s);
+ double mhold = a[0];
+ for (int i=0; i<A_0.size(); ++i){
+ mhold += A_0[i] * e_0[i];
+ }
+ m.push_back(mhold);
+
+
+ vector<double> C;
+ C.reserve(ntime_s);
+ double Chold = 0.0;
+
+ for (int i=0; i<A_0.size(); ++i){
+ double hold2 = 0.0;
+ for (int j=0; j<A_0.size(); ++j){
+ //cout << "i = " << i << " j = " << j << endl;
+ hold2 += Q_0(i,j) * A_0[j];
+ }
+
+ Chold += hold2 * A_0[i];
+ }
+ C.push_back(R[0] - Chold);
+
+
+
+ // THE FORWARD RECURSION
+ // time periods 1 to T
+ for (int tt=1; tt<ntime_s; ++tt){
+ if (STI[s][tt].size() == 0){
+ a.push_back(m[tt-1]);
+ R.push_back(C[tt-1] + tau2data[s]);
+ m.push_back(a[tt]);
+ C.push_back(R[tt]);
+ }
+ else{
+ vector<double> beta_s;
+ beta_s.reserve(STI[s][tt].size());
+ vector<double> zalpha_s;
+ zalpha_s.reserve(STI[s][tt].size());
+ for (int ii=0; ii<STI[s][tt].size(); ++ii){
+ const int i = STI[s][tt][ii];
+ beta_s.push_back(betadata[i]);
+ zalpha_s.push_back(Z[s][i] + alphadata[i]);
+ }
+ // a
+ a.push_back(m[tt-1]);
+ // R
+ R.push_back(C[tt-1] + tau2data[s]);
+ vector <double> f = beta_s;
+ for (int i=0; i<f.size(); ++i){
+ f[i] = f[i] * a[tt];
+ }
+ Matrix <> Q(beta_s.size(), beta_s.size());
+ for (int i=0; i<beta_s.size(); ++i){
+ for (int j=i; j<beta_s.size(); ++j){
+ if (i!=j){
+ Q(i,j) = Q(j,i) = beta_s[i] * beta_s[j] * R[tt];
+ }
+ else{
+ Q(i,j) = beta_s[i] * beta_s[j] * R[tt] + 1.0;
+ }
+ }
+ }
+ vector <double> e = zalpha_s;
+ for (int i=0; i<e.size(); ++i){
+ e[i] = e[i] - f[i];
+ }
+ const Matrix <> Q_inv = invpd(Q);
+
+ vector <double> A;
+ A.reserve(beta_s.size());
+ for (int i=0; i<beta_s.size(); ++i){
+ double sumholder = 0.0;
+ for (int j=0; j<beta_s.size(); ++j){
+ sumholder += beta_s[j] * Q_inv(j,i);
+ }
+ A.push_back(sumholder * R[tt]);
+ }
+ mhold = a[tt];
+ for (int i=0; i<A.size(); ++i){
+ mhold += A[i] * e[i];
+ }
+ m.push_back(mhold);
+
+ Chold = 0.0;
+ for (int i=0; i<A.size(); ++i){
+ double hold2 = 0.0;
+ for (int j=0; j<A.size(); ++j){
+ hold2 += Q(i,j) * A[j];
+ }
+ Chold += hold2 * A[i];
+ }
+ C.push_back(R[tt] - Chold);
+ }
+ } // end tt loop
+
+ if (thetaineqdata[s] == 0){
+ theta[s][ntime_s-1] = stream.rnorm(m[ntime_s-1],
+ std::sqrt(C[ntime_s-1]));
+ }
+ else if (thetaineqdata[s] > 0){
+ theta[s][ntime_s-1] = stream.rtbnorm_combo(m[ntime_s-1],
+ C[ntime_s-1], 0.0);
+ }
+ else{
+ theta[s][ntime_s-1] = stream.rtanorm_combo(m[ntime_s-1],
+ C[ntime_s-1], 0.0);
+ }
+
+
+
+
+ for (int tt=(ntime_s-2); tt>=0; --tt){
+ const double B = C[tt] * (1.0 / R[tt+1]);
+ const double h = m[tt] + B * (theta[s][tt+1] - a[tt+1]);
+ const double H = C[tt] - B * R[tt+1] * B;
+ if (thetaineqdata[s] == 0){
+ theta[s][tt] = stream.rnorm(h, std::sqrt(H));
+ }
+ else if (thetaineqdata[s] > 0){
+ theta[s][tt] = stream.rtbnorm_combo(h, H, 0.0);
+ }
+ else {
+ theta[s][tt] = stream.rtanorm_combo(h, H, 0.0);
+ }
+ } // end backwards tt loop
+
+ } // end theteqdataa[s] == -999
+
+
+
+ // sample smoothing parameters (tau2)
+ if (c0[s] > 0 && d0[s] > 0){
+ double SSE = 0.0;
+ for (int t=1; t<ntime_s; ++t){
+ SSE += std::pow((theta[s][t] - theta[s][t-1]), 2.0);
+ }
+ const double param1 = (c0[s] + ntime_s - 1) * 0.5;
+ const double param2 = (d0[s] + SSE) * 0.5;
+ tau2data[s] = stream.rigamma(param1, param2);
+ }
+
+
+
+
+ } // end s loop (end of theta tau2 loop)
+
+
+
+
+
+
+
+
+
+
+
+
+
+ // store draws
+ if (iter >= *burnin && (iter % *thin == 0)) {
+ if (*storeability == 1){
+ int internalcount = 0;
+ for (int s=0; s<*nsubj; ++s){
+ for (int t=0; t<theta[s].size(); ++t){
+ thetadraws[M(sampcount, internalcount, *nrowthetadraws)] =
+ theta[s][t];
+ ++internalcount;
+ }
+ }
+ }
+ if (*storeitem == 1){
+ for (int i=0; i<*nitems; ++i){
+ alphadraws[M(sampcount, i, *nrowalphadraws)] = alphadata[i];
+ betadraws[M(sampcount, i, *nrowbetadraws)] = betadata[i];
+ }
+ }
+ for (int s=0; s<*nsubj; ++s){
+ tau2draws[M(sampcount, s, *nrowtau2draws)] = tau2data[s];
+ }
+ ++sampcount;
+ }
+
+
+ // print output to stdout
+ // print output to stdout
+ if(*verbose > 0 && iter % *verbose == 0){
+ Rprintf("\n\nMCMCdynamicIRT1d iteration %i of %i \n", (iter+1), tot_iter);
+ /*
+ for (int t=0; t<theta[1].size(); ++t){
+ Rprintf("%f ", theta[1][t]);
+ }
+ Rprintf("\n");
+ for (int t=0; t<theta[7].size(); ++t){
+ Rprintf("%f ", theta[7][t]);
+ }
+ Rprintf("\n");
+ for (int t=0; t<theta[5].size(); ++t){
+ Rprintf("%f ", theta[5][t]);
+ }
+ Rprintf("\n");
+ for (int t=0; t<theta[8].size(); ++t){
+ Rprintf("%f ", theta[8][t]);
+ }
+ Rprintf("\n");
+ for (int t=0; t<theta[2].size(); ++t){
+ Rprintf("%f ", theta[2][t]);
+ }
+ Rprintf("\n");
+ for (int t=0; t<theta[4].size(); ++t){
+ Rprintf("%f ", theta[4][t]);
+ }
+ Rprintf("\n");
+ for (int t=0; t<theta[0].size(); ++t){
+ Rprintf("%f ", theta[0][t]);
+ }
+ Rprintf("\n");
+ for (int t=0; t<theta[3].size(); ++t){
+ Rprintf("%f ", theta[3][t]);
+ }
+ Rprintf("\n");
+ for (int t=0; t<theta[6].size(); ++t){
+ Rprintf("%f ", theta[6][t]);
+ }
+ Rprintf("\n");
+ Rprintf("\n\n");
+ Rprintf("alpha[0] = %f \n", alphadata[0]);
+ Rprintf("beta[0] = %f \n", betadata[0]);
+ Rprintf("alpha[0]/beta[0] = %f \n", alphadata[0]/betadata[0]);
+ Rprintf("tau2[1] = %f \n", tau2data[1]);
+ Rprintf("tau2[7] = %f \n", tau2data[7]);
+ Rprintf("tau2[5] = %f \n", tau2data[5]);
+ Rprintf("tau2[8] = %f \n", tau2data[8]);
+ Rprintf("tau2[2] = %f \n", tau2data[2]);
+ Rprintf("tau2[4] = %f \n", tau2data[4]);
+ Rprintf("tau2[0] = %f \n", tau2data[0]);
+ Rprintf("tau2[3] = %f \n", tau2data[3]);
+ Rprintf("tau2[6] = %f \n", tau2data[6]);
+ Rprintf("\n\n");
+ */
+ }
+
+ R_CheckUserInterrupt(); // allow user interrupts
+
+
+
+ } // end iter loop
+
+
+
+} // end MCMCdynamicIRT1d_impl function
+
+
+
+
+
+
+
+
+extern "C" {
+
+ void MCMCdynamicIRT1d(double* thetadraws, const int* nrowthetadraws,
+ const int* ncolthetadraws,
+ double* alphadraws, const int* nrowalphadraws,
+ const int* ncolalphadraws,
+ double* betadraws, const int* nrowbetadraws,
+ const int* ncolbetadraws,
+ double* tau2draws, const int* nrowtau2draws,
+ const int* ncoltau2draws,
+ const int* nsubj, const int* nitems,
+ const int* ntime,
+ const int* Ydata, const int* nrowYdata,
+ const int* ncolYdata,
+ const int* ITdata, const int* lengthITdata,
+ const int* burnin, const int* mcmc, const int* thin,
+ const int* uselecuyer, const int* seedarray,
+ const int* lecuyerstream,
+ const int* verbose,
+ const double* thetastartdata,
+ const int* lengththetastart,
+ const int* thetainfodata,
+ const int* nrowthetainfo,
+ const int* ncolthetainfo,
+ double* alphastartdata,
+ const int* lengthalphastart,
+ double* betastartdata,
+ const int* lengthbetastart,
+ double* tau2startdata,
+ const int* lengthtau2start,
+ const double* c0,
+ const int* lengthc0,
+ const double* d0,
+ const int* lengthd0,
+ const double* a0,
+ const double* A0,
+ const double* b0,
+ const double* B0,
+ const double* e0,
+ const double* E0inv,
+ const double* thetaeqdata,
+ const int* nrowthetaeq,
+ const int* ncolthetaeq,
+ const double* thetaineqdata,
+ const int* nrowthetaineq,
+ const int* ncolthetaineq,
+ const int* storeitem,
+ const int* storeability
+ ){
+
+
+ MCMCPACK_PASSRNG2MODEL(MCMCdynamicIRT1d_impl, thetadraws, nrowthetadraws,
+ ncolthetadraws, alphadraws, nrowalphadraws,
+ ncolalphadraws, betadraws, nrowbetadraws,
+ ncolbetadraws, tau2draws, nrowtau2draws,
+ ncoltau2draws, nsubj, nitems, ntime,
+ Ydata, nrowYdata, ncolYdata,
+ ITdata, lengthITdata, burnin, mcmc, thin,
+ verbose, thetastartdata, lengththetastart,
+ thetainfodata, nrowthetainfo, ncolthetainfo,
+ alphastartdata, lengthalphastart, betastartdata,
+ lengthbetastart, tau2startdata, lengthtau2start,
+ c0, lengthc0, d0, lengthd0,
+ a0, A0, b0, B0, e0, E0inv, thetaeqdata, nrowthetaeq,
+ ncolthetaeq, thetaineqdata, nrowthetaineq,
+ ncolthetaineq, storeitem, storeability );
+
+
+ } // end MCMCdynamicIRT1d
+
+} // end extern "C"
+
+
+#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