[r-cran-mnp] 17/51: Import Upstream version 2.3-4
Andreas Tille
tille at debian.org
Fri Sep 8 14:14:45 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 727a82e5443c69c40aa7e97f7105c5878c1dcdcd
Author: Andreas Tille <tille at debian.org>
Date: Fri Sep 8 15:54:48 2017 +0200
Import Upstream version 2.3-4
---
DESCRIPTION | 8 +++---
R/mnp.R | 15 ++++++-----
R/print.mnp.R | 4 +--
man/coef.mnp.Rd | 2 +-
man/cov.mnp.Rd | 2 +-
man/detergent.Rd | 5 ++--
man/mnp.Rd | 15 +++++++----
man/predict.mnp.Rd | 2 +-
man/summary.mnp.Rd | 2 +-
src/MNP.c | 9 ++++---
src/rand.c | 79 +++++++++++++++++++++++++++++++-----------------------
src/rand.h | 2 +-
12 files changed, 82 insertions(+), 63 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index ae483d3..404a67c 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
Package: MNP
-Version: 2.3-3
-Date: 2005-06-23
+Version: 2.3-4
+Date: 2005-09-06
Title: R Package for Fitting the Multinomial Probit Model
Author: Kosuke Imai <kimai at princeton.edu>,
David A. van Dyk <dvd at uci.edu>.
@@ -23,5 +23,5 @@ Description: MNP is a publicly available R package that fits the Bayesian
LazyLoad: yes
LazyData: yes
License: GPL (version 2 or later)
-URL: http://www.princeton.edu/~kimai/research/MNP.html
-Packaged: Thu Jun 23 14:11:21 2005; kimai
+URL: http://imai.princeton.edu/research/MNP.html
+Packaged: Tue Sep 6 21:11:41 2005; kimai
diff --git a/R/mnp.R b/R/mnp.R
index f5bb133..e17c1d2 100644
--- a/R/mnp.R
+++ b/R/mnp.R
@@ -1,12 +1,12 @@
mnp <- function(formula, data = parent.frame(), choiceX = NULL,
cXnames = NULL, base = NULL, latent = FALSE,
- n.draws = 5000, p.var = "Inf", p.df = n.dim+1,
- p.scale = 1, coef.start = 0, cov.start = 1,
- burnin = 0, thin = 0, verbose = FALSE) {
+ invcdf = FALSE, n.draws = 5000, p.var = "Inf",
+ p.df = n.dim+1, p.scale = 1, coef.start = 0,
+ cov.start = 1, burnin = 0, thin = 0, verbose = FALSE) {
call <- match.call()
mf <- match.call(expand = FALSE)
mf$choiceX <- mf$cXnames <- mf$base <- mf$n.draws <- mf$latent <-
- mf$p.var <- mf$p.df <- mf$p.scale <- mf$coef.start <-
+ mf$p.var <- mf$p.df <- mf$p.scale <- mf$coef.start <- mf$invcdf <-
mf$cov.start <- mf$verbose <- mf$burnin <- mf$thin <- NULL
mf[[1]] <- as.name("model.frame")
mf$na.action <- 'na.pass'
@@ -147,7 +147,8 @@ mnp <- function(formula, data = parent.frame(), choiceX = NULL,
as.double(p.mean), as.double(p.prec), as.integer(p.df),
as.double(p.scale*p.alpha0), as.double(X), as.integer(Y),
as.double(coef.start), as.double(cov.start),
- as.integer(p.imp), as.integer(burnin), as.integer(keep),
+ as.integer(p.imp), as.integer(invcdf),
+ as.integer(burnin), as.integer(keep),
as.integer(verbose), as.integer(MoP), as.integer(latent),
pdStore = double(n.par*floor((n.draws-burnin)/keep)),
PACKAGE="MNP")$pdStore
@@ -168,8 +169,8 @@ mnp <- function(formula, data = parent.frame(), choiceX = NULL,
## returning the object
res <- list(param = param, x = X, y = Y, w = W, call = call, alt = lev,
- n.alt = p, base = base, p.mean = if(p.imp) NULL else p.mean,
- p.var = p.var,
+ n.alt = p, base = base, invcdf = invcdf,
+ p.mean = if(p.imp) NULL else p.mean, p.var = p.var,
p.df = p.df, p.scale = p.scale, burnin = burnin, thin = thin)
class(res) <- "mnp"
return(res)
diff --git a/R/print.mnp.R b/R/print.mnp.R
index e71a382..1d04d14 100644
--- a/R/print.mnp.R
+++ b/R/print.mnp.R
@@ -3,11 +3,11 @@ print.mnp <- function (x, digits = max(3, getOption("digits") - 3), ...)
cat("\nCall:\n", deparse(x$call), "\n\n", sep = "")
param <- apply(x$param, 2, mean)
if (length(param)) {
- cat("Parameters:\n")
+ cat("Parameter estimates (posterior means):\n")
print.default(format(param, digits = digits), print.gap = 2,
quote = FALSE)
}
- else cat("No parameters\n")
+ else cat("No parameter estimates\n")
cat("\n")
invisible(x)
}
diff --git a/man/coef.mnp.Rd b/man/coef.mnp.Rd
index 882f308..68cc010 100644
--- a/man/coef.mnp.Rd
+++ b/man/coef.mnp.Rd
@@ -35,7 +35,7 @@
}
\seealso{\code{mnp}, \code{cov.mnp}; MNP home page at
- \url{http://www.princeton.edu/~kimai/research/MNP.html}}
+ \url{http://imai.princeton.edu/research/MNP.html}}
\author{
Kosuke Imai, Department of Politics, Princeton University
diff --git a/man/cov.mnp.Rd b/man/cov.mnp.Rd
index a86053f..04ed5ae 100644
--- a/man/cov.mnp.Rd
+++ b/man/cov.mnp.Rd
@@ -32,7 +32,7 @@
}
\seealso{\code{mnp}, \code{coef.mnp}; MNP home page at
- \url{http://www.princeton.edu/~kimai/research/MNP.html}}
+ \url{http://imai.princeton.edu/research/MNP.html}}
\author{
Kosuke Imai, Department of Politics, Princeton University
diff --git a/man/detergent.Rd b/man/detergent.Rd
index 7c260a9..d3f812a 100644
--- a/man/detergent.Rd
+++ b/man/detergent.Rd
@@ -26,6 +26,7 @@
\references{Chintagunta, P. K. and Prasad, A. R. (1998) \dQuote{An Empirical
Investigation of the `Dynamic McFadden' Model of Purchase Timing and
- Brand Choice: Implications for Market Structure}. \emph{Journal of Business and
- Economic Statistics} vol. 16 no. 1 pp.2-12.}
+ Brand Choice: Implications for Market Structure}. \emph{Journal of
+ Business and Economic Statistics} vol. 16 no. 1 pp.2-12.
+}
\keyword{datasets}
diff --git a/man/mnp.Rd b/man/mnp.Rd
index f58a378..25d2aa1 100644
--- a/man/mnp.Rd
+++ b/man/mnp.Rd
@@ -15,9 +15,9 @@
\usage{
mnp(formula, data = parent.frame(), choiceX = NULL, cXnames = NULL,
- base = NULL, latent = FALSE, n.draws = 5000, p.var = "Inf",
- p.df = n.dim+1, p.scale = 1, coef.start = 0, cov.start = 1,
- burnin = 0, thin = 0, verbose = FALSE)
+ base = NULL, latent = FALSE, invcdf = FALSE, n.draws = 5000,
+ p.var = "Inf", p.df = n.dim+1, p.scale = 1, coef.start = 0,
+ cov.start = 1, burnin = 0, thin = 0, verbose = FALSE)
}
\arguments{
@@ -48,6 +48,10 @@ mnp(formula, data = parent.frame(), choiceX = NULL, cXnames = NULL,
be returned. See Imai and van Dyk (2005) for the notation. The
default is \code{FALSE}.
}
+ \item{invcdf}{logical. If \code{TRUE}, then the inverse cdf method is
+ used for truncated normal sampling. If \code{FALSE}, then the
+ rejection sampling method is used. The default is \code{FALSE}.
+ }
\item{n.draws}{A positive integer. The number of MCMC draws. The
default is \code{5000}.
}
@@ -187,6 +191,7 @@ predict(res2, newdata = japan[10,], type = "prob")
\item{alt}{The names of alternatives.}
\item{n.alt}{The total number of alternatives.}
\item{base}{The base category used for fitting.}
+ \item{invcdf}{The value of \code{invcdf}.}
\item{p.var}{The prior variance for the coefficients.}
\item{p.df}{The prior degrees of freedom parameter for the covariance matrix.}
\item{p.scale}{The prior scale matrix for the covariance matrix.}
@@ -208,7 +213,7 @@ predict(res2, newdata = japan[10,], type = "prob")
\author{
Kosuke Imai, Department of Politics, Princeton University
- \email{kimai at Princeton.Edu}, \url{http://www.princeton.edu/~kimai};
+ \email{kimai at Princeton.Edu}, \url{http://imai.princeton.edu};
Jordan R. Vance, Princeton University; David A. van Dyk, Department of
Statistics, University of California, Irvine \email{dvd at uci.edu},
\url{http://www.ics.uci.edu/~dvd}.
@@ -216,7 +221,7 @@ predict(res2, newdata = japan[10,], type = "prob")
\seealso{\code{coef.mnp}, \code{cov.mnp}, \code{predict.mnp},
\code{summary.mnp}; MNP home page at
- \url{http://www.princeton.edu/~kimai/research/MNP.html}}
+ \url{http://imai.princeton.edu/research/MNP.html}}
\keyword{models}
diff --git a/man/predict.mnp.Rd b/man/predict.mnp.Rd
index 96e020b..0018dcd 100644
--- a/man/predict.mnp.Rd
+++ b/man/predict.mnp.Rd
@@ -96,7 +96,7 @@
}
\seealso{\code{mnp}; MNP home page at
- \url{http://www.princeton.edu/~kimai/research/MNP.html}}
+ \url{http://imai.princeton.edu/research/MNP.html}}
\author{
Kosuke Imai, Department of Politics, Princeton University
diff --git a/man/summary.mnp.Rd b/man/summary.mnp.Rd
index 1f31b16..9bd63cf 100644
--- a/man/summary.mnp.Rd
+++ b/man/summary.mnp.Rd
@@ -42,7 +42,7 @@
}
\seealso{\code{mnp}; MNP home page at
- \url{http://www.princeton.edu/~kimai/research/MNP.html}}
+ \url{http://imai.princeton.edu/research/MNP.html}}
\author{
Kosuke Imai, Department of Politics, Princeton University
diff --git a/src/MNP.c b/src/MNP.c
index aeec2e1..d530780 100644
--- a/src/MNP.c
+++ b/src/MNP.c
@@ -18,11 +18,12 @@ void cMNPgibbs(int *piNDim, int *piNCov, int *piNSamp, int *piNGen,
double *pdA0, int *piNu0, double *pdS, double *pdX,
int *y, /* response variable: -1 for missing */
double *pdbeta, double *pdSigma, int *piImp,
+ int *invcdf, /* use inverse cdf for TruncNorm? */
int *piBurnin, /* the number of burnin */
int *piKeep,
int *verbose, /* 1 if extra print is needed */
int *piMoP, /* 1 if Multinomial ordered Probit */
- int *latent, /* 1 if W is stored */
+ int *latent, /* 1 if W is stored */
double *pdStore){
/* paramerters from R */
@@ -245,7 +246,7 @@ void cMNPgibbs(int *piNDim, int *piNCov, int *piNSamp, int *piNGen,
else {
if(Y[i][j]==0) minw = cmean - 1000*sqrt(cvar);
if(Y[i][j]==maxy[i]) maxw = cmean + 1000*sqrt(cvar);
- W[i][j]=TruncNorm(minw,maxw,cmean,cvar);
+ W[i][j]=TruncNorm(minw,maxw,cmean,cvar,*invcdf);
}
/* printf("%14g\n", W[i][j]); */
X[i*n_dim+j][n_cov]=W[i][j];
@@ -277,9 +278,9 @@ void cMNPgibbs(int *piNDim, int *piNCov, int *piNSamp, int *piNGen,
if(y[i]==-1)
W[i][j]=cmean+norm_rand()*sqrt(cvar);
else if(y[i]==(j+1))
- W[i][j]=TruncNorm(maxw,cmean+1000*sqrt(cvar),cmean,cvar);
+ W[i][j]=TruncNorm(maxw,cmean+1000*sqrt(cvar),cmean,cvar,*invcdf);
else
- W[i][j]=TruncNorm(cmean-1000*sqrt(cvar),maxw,cmean,cvar);
+ W[i][j]=TruncNorm(cmean-1000*sqrt(cvar),maxw,cmean,cvar,*invcdf);
X[i*n_dim+j][n_cov]=W[i][j];
X[i*n_dim+j][n_cov]*=sqrt(alpha2);
}
diff --git a/src/rand.c b/src/rand.c
index bb35be4..103dd75 100644
--- a/src/rand.c
+++ b/src/rand.c
@@ -15,51 +15,62 @@
#include "subroutines.h"
#include "rand.h"
+
/* Sample from a univariate truncated Normal distribution
- (truncated both from above and below):
+ (truncated both from above and below): choose either inverse cdf
+ method or rejection sampling method. For rejection sampling,
if the range is too far from mu, it uses standard rejection
sampling algorithm with exponential envelope function. */
double TruncNorm(
double lb, /* lower bound */
double ub, /* upper bound */
double mu, /* mean */
- double var /* variance */
+ double var, /* variance */
+ int invcdf /* use inverse cdf method? */
) {
-
- double stlb, stub, temp, M, u, z, exp_par;
- int flag=0; /* 1 if stlb, stub <-2 */
- stlb = (lb-mu)/sqrt(var); /* standardized lower bound */
- stub = (ub-mu)/sqrt(var); /* standardized upper bound */
+ double z;
+ double sigma = sqrt(var);
+ double stlb = (lb-mu)/sigma; /* standardized lower bound */
+ double stub = (ub-mu)/sigma; /* standardized upper bound */
if(stlb >= stub)
- error("TurncNorm: lower bound is greater than upper bound\n");
- if(stub<=-2){
- flag=1;
- temp=stub;
- stub=-stlb;
- stlb=-temp;
+ error("TruncNorm: lower bound is greater than upper bound\n");
+ if (invcdf) { /* inverse cdf method */
+ z = qnorm(runif(pnorm(stlb, 0, 1, 1, 0), pnorm(stub, 0, 1, 1, 0)),
+ 0, 1, 1, 0);
+ }
+ else { /* rejection sampling method */
+ double tol=2.0;
+ double temp, M, u, exp_par;
+ int flag=0; /* 1 if stlb, stub <-tol */
+ if(stub<=-tol){
+ flag=1;
+ temp=stub;
+ stub=-stlb;
+ stlb=-temp;
+ }
+ if(stlb>=tol){
+ exp_par=stlb;
+ while(pexp(stub,1/exp_par,1,0) - pexp(stlb,1/exp_par,1,0) < 0.000001)
+ exp_par/=2.0;
+ if(dnorm(stlb,0,1,1) - dexp(stlb,1/exp_par,1) >=
+ dnorm(stub,0,1,1) - dexp(stub,1/exp_par,1))
+ M=exp(dnorm(stlb,0,1,1) - dexp(stlb,1/exp_par,1));
+ else
+ M=exp(dnorm(stub,0,1,1) - dexp(stub,1/exp_par,1));
+ do{
+ u=unif_rand();
+ z=-log(1-u*(pexp(stub,1/exp_par,1,0)-pexp(stlb,1/exp_par,1,0))
+ -pexp(stlb,1/exp_par,1,0))/exp_par;
+ }while(unif_rand() > exp(dnorm(z,0,1,1)-dexp(z,1/exp_par,1))/M );
+ if(flag==1) z=-z;
+ }
+ else{
+ do z=norm_rand();
+ while( z<stlb || z>stub );
+ }
}
- if(stlb>=2){
- exp_par=stlb;
- while(pexp(stub,1/exp_par,1,0) - pexp(stlb,1/exp_par,1,0) < 0.000001)
- exp_par/=2.0;
- if(dnorm(stlb,0,1,1) - dexp(stlb,1/exp_par,1) >=
- dnorm(stub,0,1,1) - dexp(stub,1/exp_par,1))
- M=exp(dnorm(stlb,0,1,1) - dexp(stlb,1/exp_par,1));
- else
- M=exp(dnorm(stub,0,1,1) - dexp(stub,1/exp_par,1));
- do{ /* sample from Exponential by inverse cdf method */
- u=unif_rand();
- z=-log(1-u*(pexp(stub,1/exp_par,1,0)-pexp(stlb,1/exp_par,1,0))
- -pexp(stlb,1/exp_par,1,0))/exp_par;
- }while(unif_rand() > exp(dnorm(z,0,1,1)-dexp(z,1/exp_par,1))/M );
- if(flag==1) z=-z;
- } /* if(stlb>=2) */
- else{ /* if the range is not extreme */
- do z=norm_rand();
- while( z<stlb || z>stub );
- }/* else */
- return(z*sqrt(var) + mu);
+ return(z*sigma + mu);
}
diff --git a/src/rand.h b/src/rand.h
index 857d596..b8d793e 100644
--- a/src/rand.h
+++ b/src/rand.h
@@ -5,7 +5,7 @@
Copyright: GPL version 2 or later.
*******************************************************************/
-double TruncNorm(double lb, double ub, double mu, double var);
+double TruncNorm(double lb, double ub, double mu, double var, int invcdf);
void rMVN(double *Sample, double *mean, double **inv_Var, int size);
void rWish(double **Sample, double **S, int df, int size);
--
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