[r-cran-mnp] 32/51: Import Upstream version 2.5-1
Andreas Tille
tille at debian.org
Fri Sep 8 14:14:47 UTC 2017
This is an automated email from the git hooks/post-receive script.
tille pushed a commit to branch master
in repository r-cran-mnp.
commit 488605ffa949aa62badf37718274dcb6ddcef29a
Author: Andreas Tille <tille at debian.org>
Date: Fri Sep 8 15:54:55 2017 +0200
Import Upstream version 2.5-1
---
DESCRIPTION | 6 +-
R/predict.mnp.R | 85 ++++++++++---------------
man/mnp.Rd | 8 ++-
man/predict.mnp.Rd | 57 +++++++----------
src/MNP.c | 179 +++++++++++++++++++++++++++++++++++++++++++++++++++++
5 files changed, 244 insertions(+), 91 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index 925d80c..80897b6 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
Package: MNP
-Version: 2.4-2
-Date: 2006-10-17
+Version: 2.5-1
+Date: 2006-11-21
Title: R Package for Fitting the Multinomial Probit Model
Author: Kosuke Imai <kimai at princeton.edu>,
David A. van Dyk <dvd at uci.edu>.
@@ -24,4 +24,4 @@ LazyLoad: yes
LazyData: yes
License: GPL (version 2 or later)
URL: http://imai.princeton.edu/research/MNP.html
-Packaged: Tue Oct 17 20:57:21 2006; kimai
+Packaged: Tue Nov 21 11:06:30 2006; kimai
diff --git a/R/predict.mnp.R b/R/predict.mnp.R
index 2f17ff5..a57fddb 100644
--- a/R/predict.mnp.R
+++ b/R/predict.mnp.R
@@ -1,10 +1,12 @@
-predict.mnp <- function(object, newdata = NULL, newdraw = NULL, moredraw = 1,
- type = c("prob", "choice", "order", "latent"),
+predict.mnp <- function(object, newdata = NULL, newdraw = NULL, n.draws = 1,
+ type = c("prob", "choice", "order"),
verbose = FALSE, ...){
- if (NA %in% match(type, c("prob", "choice", "order", "latent")))
+ if (NA %in% match(type, c("prob", "choice", "order")))
stop("Invalid input for `type'.")
-
+ if (n.draws < 1)
+ stop("Invalid input for `n.draws'.")
+
p <- object$n.alt
if (is.null(newdraw))
param <- object$param
@@ -12,7 +14,7 @@ predict.mnp <- function(object, newdata = NULL, newdraw = NULL, moredraw = 1,
param <- newdraw
coef <- coef(object)
n.cov <- ncol(coef)
- n.draws <- nrow(param)
+ n.mcmc <- nrow(coef)
cov <- param[,(n.cov+1):ncol(param)]
## get X matrix ready
@@ -31,7 +33,6 @@ predict.mnp <- function(object, newdata = NULL, newdraw = NULL, moredraw = 1,
else if (sum(is.na(x))>0)
stop("Invalid input for `newdata'.")
}
-
n.obs <- nrow(x)
if (verbose) {
if (n.obs == 1)
@@ -39,64 +40,44 @@ predict.mnp <- function(object, newdata = NULL, newdraw = NULL, moredraw = 1,
else
cat("There are", n.obs, "observations to predict. Please wait...\n")
}
-
+
alt <- object$alt
if (object$base != alt[1])
alt <- c(object$base, alt[1:(length(alt)-1)])
-
- ## computing W
- W <- array(NA, dim=c(p-1, n.obs, n.draws, moredraw),
- dimnames=c(alt[2:p], NULL, 1:n.draws, 1:moredraw))
- tmp <- floor(n.draws/10)
- inc <- 1
- Sigma <- cov.mnp(object)
- for (i in 1:n.draws) {
- for (j in 1:n.obs)
- W[,j,i,] <- c(matrix(x[j,], ncol=n.cov) %*% matrix(coef[i,], ncol = 1)) +
- mvrnorm(moredraw, mu = rep(0, p-1), Sigma = Sigma[,,i])
- if (i == inc*tmp & verbose) {
- cat("", inc*10, "percent done.\n")
- inc <- inc + 1
- }
- }
- ans <- list()
- if ("latent" %in% type)
- ans$w <- W
- else
- ans$w <- NULL
- ## computing Y
- Y <- array(NA, dim = c(n.obs, n.draws, moredraw),
- dimnames=list(NULL, 1:n.draws, 1:moredraw))
- O <- array(NA, dim=c(p, n.obs, n.draws, moredraw),
- dimnames=list(alt, NULL, 1:n.draws, 1:moredraw))
- for (i in 1:n.obs)
- for (j in 1:n.draws) {
- for (k in 1:moredraw) {
- Y[i,j,k] <- alt[match(max(c(0, W[,i,j,k])), c(0,W[,i,j,k]))]
- O[,i,j,k] <- rank(c(0, -W[,i,j,k]))
- }
- }
+ res <- .C("predict", as.double(x), as.integer(n.obs),
+ as.double(coef), as.double(cov), as.integer(p-1),
+ as.integer(n.cov), as.integer(n.mcmc),
+ as.integer(n.draws), as.integer(verbose),
+ prob = if (n.draws > 1) double(n.obs*n.mcmc*p)
+ else double(n.obs*p),
+ choice = double(n.obs*n.mcmc),
+ order = double(n.obs*n.mcmc*p),
+ PACKAGE = "MNP")
+
+ ans <- list()
if ("choice" %in% type)
- ans$y <- Y
+ ans$y <- matrix(res$choice, ncol=n.mcmc, nrow = n.obs,
+ byrow = TRUE, dimnames = list(rownames(x), 1:n.mcmc))
else
ans$y <- NULL
if ("order" %in% type)
- ans$o <- O
+ ans$o <- aperm(array(res$order, dim=c(p, n.mcmc, n.obs),
+ dimnames = c(alt, 1:n.mcmc,
+ rownames(x))), c(3,1,2))
else
ans$o <- NULL
-
- ## computing P
- if ("prob" %in% type) {
- P <- array(NA, dim = c(n.obs, p, moredraw),
- dimnames = c(NULL, alt, 1:moredraw))
- for (i in 1:p)
- for (j in 1:moredraw)
- P[,i,j] <- apply(Y[,,j,drop=FALSE]==alt[i], 1, mean)
- ans$p <- P
- }
+ if ("prob" %in% type)
+ if (n.draws > 1)
+ ans$p <- aperm(array(res$prob, dim = c(p, n.mcmc, n.obs),
+ dimnames = c(alt, 1:nrow(coef),
+ rownames(x))), c(3,1,2))
+ else
+ ans$p <- matrix(res$prob, nrow = n.obs, ncol = p, byrow = TRUE,
+ dimnames = list(rownames(x), alt))
else
ans$p <- NULL
+ class(ans) <- "predict.mnp"
return(ans)
}
diff --git a/man/mnp.Rd b/man/mnp.Rd
index a9457e6..5240a3e 100644
--- a/man/mnp.Rd
+++ b/man/mnp.Rd
@@ -161,8 +161,8 @@ res1 <- mnp(choice ~ 1, choiceX = list(Surf=SurfPrice, Tide=TidePrice,
thin = 3, verbose = TRUE)
## summarize the results
summary(res1)
-## calculate the predicted probabilities for the first 5 observations
-predict(res1, newdata = detergent[1:3,], type="prob", verbose = TRUE)
+## calculate the quantities of interest for the first 3 observations
+pre1 <- predict(res1, newdata = detergent[1:3,])
## load the Japanese election data
data(japan)
@@ -172,7 +172,9 @@ res2 <- mnp(cbind(LDP, NFP, SKG, JCP) ~ gender + education + age, data = japan,
## summarize the results
summary(res2)
## calculate the predicted probabilities for the 10th observation
-predict(res2, newdata = japan[10,], type = "prob", verbose = TRUE)
+## averaging over 100 additional Monte Carlo draws given each of MCMC draw.
+pre2 <- predict(res2, newdata = japan[10,], type = "prob", n.draws = 100,
+ verbose = TRUE)
}
\value{
diff --git a/man/predict.mnp.Rd b/man/predict.mnp.Rd
index 6c5f108..3d7c627 100644
--- a/man/predict.mnp.Rd
+++ b/man/predict.mnp.Rd
@@ -10,8 +10,8 @@
}
\usage{
- \method{predict}{mnp}(object, newdata = NULL, newdraw = NULL, moredraw = 1,
- type = c("prob", "choice", "order", "latent"), verbose = FALSE, ...)
+ \method{predict}{mnp}(object, newdata = NULL, newdraw = NULL, n.draws = 1,
+ type = c("prob", "choice", "order"), verbose = FALSE, ...)
}
\arguments{
@@ -26,10 +26,11 @@
posterior predictions. The default is the original MCMC draws stored
in \code{object}.
}
- \item{moredraw}{The number of additional draws of latent variables for
- each of MCMC draws. Given a particular MCMC draw of coefficients and
- covariance matrix, the specified number of latent variables will be
- sampled from the multivariate normal distribution. This will be
+ \item{n.draws}{The number of additional Monte Carlo draws given each
+ MCMC draw of coefficients and covariance matrix. The specified
+ number of latent variables will be sampled from the multivariate
+ normal distribution, and the quantities of interest will be
+ calculated by averaging over these draws. This will be
particularly useful calculating the uncertainty of predicted
probabilities. The default is \code{1}.
}
@@ -39,11 +40,8 @@
the most preferred choice among the choice set.
\code{type = "choice"} returns the Monte Carlo sample of the
most preferred choice,
- \code{type = "order"} returns the Monte Carlo sample of the
+ and \code{type = "order"} returns the Monte Carlo sample of the
ordered preferences,
- and \code{type = "latent"} returns the Monte Carlo sample of the
- predictive values of the latent variable. The default is to return
- all four types of posterior predictions.
}
\item{verbose}{logical. If \code{TRUE}, helpful messages along with a
progress report on the Monte Carlo sampling from the posterior
@@ -73,40 +71,33 @@
}
\value{
- \code{predict.mnp} yields a list containing at least one of the
+ \code{predict.mnp} yields a list of class
+ \code{predict.mnp} containing at least one of the
following elements:
- \item{o}{A four dimensional array of the Monte Carlo sample from the
+ \item{o}{A three dimensional array of the Monte Carlo sample from the
posterior predictive distribution of the ordered preferences. The
first dimension corresponds to the alternatives in the choice set,
the second dimension corresponds to the rows of \code{newdata} (or
the original data set if \code{newdata} is left unspecified),
- the third dimension indexes the Monte Carlo sample, and the fourth
- dimension is the number of additional draws given by
- \code{moredraw}.}
- \item{p}{A four dimensional array of the posterior predictive
+ and the third dimension indexes the Monte Carlo sample. If
+ \code{n.draws} is greater than 1, then each entry will be an average
+ over these additional draws.
+ }
+ \item{p}{A two or three dimensional array of the posterior predictive
probabilities for each alternative in the choice set being most
preferred. The first demension corresponds to the rows of
\code{newdata} (or the original data set if \code{newdata} is left
- unspecified), the second dimension corresponds
- to the alternatives in the choice set, and the third diemsion
- represents the number of additional draws given by \code{moredraw}.
+ unspecified), and the second dimension corresponds
+ to the alternatives in the choice set. If \code{n.draws} is greater
+ than 1, then the third diemsion exists and indexes the Monte Carlo
+ sample.
}
- \item{y}{A three dimensional array of the Monte Carlo sample from the
+ \item{y}{A matrix of the Monte Carlo sample from the
posterior predictive distribution of the most preferred choice. The
first dimension correspond to the rows of \code{newdata} (or the
- original data set if \code{newdata} is left unspecified), the
- second dimension indexes the Monte Carlo sample, and the third
- dimension represents the number of additional draws given by
- \code{moredraw}.
- }
- \item{w}{A four dimensional array of the Monte Carlo sample from the
- posterior predictive distribution of the latent
- variable. The first dimension corresponds to the alternatives in the
- choice set, the second dimension corresponds to the rows of
- \code{newdata} (or the original data set if \code{newdata} is left
- unspecified), the third dimension indexes the Monte Carlo sample,
- and the four dimension represents the number of additional draws
- given by \code{moredraw}.
+ original data set if \code{newdata} is left unspecified), and the
+ second dimension indexes the Monte Carlo sample. \code{n.draws} will
+ be set to 1 when computing this quantity of interest.
}
}
diff --git a/src/MNP.c b/src/MNP.c
index 8b3f389..b8dddea 100644
--- a/src/MNP.c
+++ b/src/MNP.c
@@ -419,4 +419,183 @@ void cMNPgibbs(int *piNDim, int *piNCov, int *piNSamp, int *piNGen,
Free3DMatrix(PerSig, n_dim, n_dim);
} /* main */
+
+
+/* unitility function */
+void R_max_col2(double **matrix,
+ int nr,
+ int nc,
+ int *maxes,
+ int ties_meth) {
+
+ int *ncol = intArray(1);
+ int *nrow = intArray(1);
+ int *ties = intArray(1);
+ int *itmp = intArray(1);
+ double *tmp = doubleArray(nr*nc);
+ int i, j, k;
+
+ ncol[0] = nc; nrow[0] = nr; ties[0] = ties_meth;
+ i = 0;
+ for (k = 0; k < nc; k++)
+ for (j = 0; j < nr; j++)
+ tmp[i++] = matrix[j][k];
+
+ R_max_col(tmp, nrow, ncol, maxes, ties);
+
+ free(ncol);
+ free(nrow);
+ free(itmp);
+ free(tmp);
+
+}
+
+void predict(double *dX, /* X matrix */
+ int *nobs, /* number of observations */
+ double *dcoef, /* coefficients */
+ double *dSigma, /* covariances */
+ int *ndims, /* number of dimensions */
+ int *ncovs, /* number of covariates */
+ int *ndraws, /* number of MCMC draws */
+ int *moredraws, /* number of extra draws */
+ int *verbose,
+ double *prob, /* probability output */
+ double *choice, /* choice output */
+ double *order /* order output */
+ ) {
+
+ int n_samp = *nobs;
+ int n_draw = *ndraws;
+ int n_dim = *ndims;
+ int n_cov = *ncovs;
+ int n_extra = *moredraws;
+
+ double **X = doubleMatrix(n_samp*n_dim, n_cov);
+ double *Xbeta = doubleArray(n_cov);
+ double **W = doubleMatrix(n_extra, n_dim+1);
+ double **beta = doubleMatrix(n_draw, n_cov);
+ double ***Sigma = doubleMatrix3D(n_draw, n_dim, n_dim);
+
+ int i, j, k, main_loop, itemp, itempP, itempO, itempC;
+ int total = n_extra*n_samp*n_draw;
+ int count, progress = 1;
+ int itempQ = ftrunc((double) total/10);
+ int *maxdim = intArray(n_extra);
+ int *ind = intArray(n_dim+1);
+ int *sumorder = intArray(n_dim+1);
+ int *probTemp = intArray(n_dim+1);
+
+ /** reading the data */
+ itemp = 0;
+ for (k = 0; k < n_cov; k++)
+ for (j = 0; j < n_dim; j++)
+ for (i = 0; i < n_samp; i++)
+ X[i*n_dim+j][k] = dX[itemp++];
+
+ /** reading the MCMC draws **/
+ itemp = 0;
+ for (k = 0; k < n_cov; k++)
+ for (j = 0; j < n_draw; j++)
+ beta[j][k] = dcoef[itemp++];
+
+ itemp = 0;
+ for (k = 0; k < n_dim; k++)
+ for (j = k; j < n_dim; j++)
+ for (i = 0; i < n_draw; i++) {
+ Sigma[i][j][k] = dSigma[itemp++];
+ Sigma[i][k][j] = Sigma[i][j][k];
+ }
+
+ /** get random seed **/
+ GetRNGstate();
+
+ /** Posterior predictive simulations **/
+ itempC = 0; itempO = 0; itempP = 0; count = 0;
+ itempQ = 0;
+
+ /* loop for observations */
+ for (i = 0; i < n_samp; i++) {
+ if (n_extra == 1) {
+ for (j = 0; j <= n_dim; j++)
+ probTemp[j] = 0;
+ }
+ /* loop for MCMC draws */
+ for (main_loop = 0; main_loop < n_draw; main_loop++) {
+ if (n_extra > 1) {
+ for (j = 0; j <= n_dim; j++)
+ probTemp[j] = 0;
+ }
+ /* compute the mean for each dimension */
+ for (j = 0; j < n_dim; j++) {
+ Xbeta[j] = 0;
+ for (k = 0; k < n_cov; k++)
+ Xbeta[j] += X[i*n_dim+j][k] * beta[main_loop][k];
+ }
+ /* sample W */
+ for (j = 0; j < n_extra; j++) {
+ rMVN(W[j], Xbeta, Sigma[main_loop], n_dim);
+ W[j][n_dim] = 0;
+ }
+ /* which dimension is max for each of n_extra W draws? */
+ /* PdoubleMatrix(W, n_extra, n_dim+1); */
+ R_max_col2(W, n_extra, n_dim+1, maxdim, 1);
+ /* PintArray(maxdim, n_extra); */
+ /* order */
+ for (j = 0; j < n_extra; j++) {
+ for (k = 0; k <= n_dim; k++) {
+ ind[k] = k; sumorder[k] = 0;
+ }
+ rsort_with_index(W[j], ind, n_dim+1);
+ /* PintArray(ind, n_dim+1); */
+ for (k = 0; k <= n_dim; k++)
+ sumorder[ind[k]] += (n_dim-k);
+ if(*verbose) {
+ if(count == itempQ) {
+ Rprintf("%3d percent done.\n", progress*10);
+ itempQ += ftrunc((double) total/10); progress++;
+ R_FlushConsole();
+ }
+ count++;
+ }
+ }
+ if (n_extra > 1) {
+ /* store probability and mean order */
+ for (j = 0; j <= n_dim; j++) {
+ itemp = 0;
+ for (k = 0; k < n_extra; k++)
+ if (maxdim[k] == (j+1))
+ itemp++;
+ prob[itempP++] = ((double) itemp / (double) n_extra);
+ order[itempO++] = ((double) sumorder[j] / (double) n_extra);
+ }
+ } else {
+ /* store choice */
+ for (j = 0; j <= n_dim; j++) {
+ if (maxdim[0] == (j+1)) {
+ choice[itempC++] = j;
+ probTemp[j]++;
+ }
+ order[itempO++] = sumorder[j];
+ }
+ }
+ }
+ if (n_extra == 1)
+ for (j = 0; j <= n_dim; j++)
+ prob[itempP++] = ((double) probTemp[j] / (double) n_draw);
+ }
+
+ /** write out the random seed **/
+ PutRNGstate();
+
+ /* freeing memory */
+ FreeMatrix(X, n_samp*n_dim);
+ free(Xbeta);
+ FreeMatrix(W, n_extra);
+ FreeMatrix(beta, n_draw);
+ Free3DMatrix(Sigma, n_draw, n_dim);
+ free(maxdim);
+ free(ind);
+ free(sumorder);
+ free(probTemp);
+}
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/r-cran-mnp.git
More information about the debian-science-commits
mailing list