[r-cran-mcmcpack] 38/90: Imported Upstream version 1.0-3
Andreas Tille
tille at debian.org
Fri Dec 16 09:07:38 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 cddd909c0df50ba7718abf0a1f900ee08de99cc6
Author: Andreas Tille <tille at debian.org>
Date: Fri Dec 16 08:07:15 2016 +0100
Imported Upstream version 1.0-3
---
DESCRIPTION | 8 +--
HISTORY | 6 +-
NAMESPACE | 1 +
R/MCMCquantreg.R | 84 +++++++++++++++++++++++++
man/MCMCirtHier1d.Rd | 4 +-
man/MCMCquantreg.Rd | 156 ++++++++++++++++++++++++++++++++++++++++++++++
src/MCMCfcds.h | 172 ++++++++++++++++++++++++++++++++++++++++++++++++++-
src/MCMCmedreg.cc | 138 +++++++++++++++++++++++++++++++++++++++++
src/MCMCquantreg.cc | 140 +++++++++++++++++++++++++++++++++++++++++
9 files changed, 698 insertions(+), 11 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index 63ef4be..88533af 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
Package: MCMCpack
-Version: 1.0-2
-Date: 2009-7-17
+Version: 1.0-3
+Date: 2009-7-27
Title: Markov chain Monte Carlo (MCMC) Package
Author: Andrew D. Martin <admartin at wustl.edu>, Kevin M. Quinn
<kevin_quinn at harvard.edu>, Jong Hee Park <jhp at uchicago.edu>
@@ -19,6 +19,6 @@ Description: This package contains functions to perform Bayesian
License: GPL-2
SystemRequirements: gcc (>= 4.0)
URL: http://mcmcpack.wustl.edu
-Packaged: 2009-07-17 18:52:45 UTC; adm
+Packaged: 2009-07-28 17:41:00 UTC; adm
Repository: CRAN
-Date/Publication: 2009-07-18 15:22:20
+Date/Publication: 2009-07-30 12:51:08
diff --git a/HISTORY b/HISTORY
index 698091d..984aa75 100644
--- a/HISTORY
+++ b/HISTORY
@@ -2,6 +2,10 @@
// Changes and Bug Fixes
//
+1.0-2 to 1.0-3
+ * Mix fix to hierarchical IRT documentation [written by Mike Malecki]
+ * added Bayesian quantile regression [contributed by Craig Reed]
+
1.0-1 to 1.0-2
* Added Poisson regression changepoint analysis MCMCpoissonChange() [JHP]
* Old MCMCpoissonChangepoint() is deprecated [JHP]
@@ -10,8 +14,6 @@
* plotChangepoint() function modified to support new models
* Added a number of new helper functions in btsutil.R [JHP]
-
-
0.9-6 to 1.0-1
* added one-dimensional dynamic IRT model
* added Rehnquist Court data to illustrate dynamic IRT model
diff --git a/NAMESPACE b/NAMESPACE
index cda65b0..fb9cbc9 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -34,6 +34,7 @@ export(
MCMCpoisson,
MCMCpoissonChange,
MCMCprobit,
+ MCMCquantreg,
MCMCregress,
MCMCSVDreg,
MCMCtobit,
diff --git a/R/MCMCquantreg.R b/R/MCMCquantreg.R
new file mode 100644
index 0000000..865af80
--- /dev/null
+++ b/R/MCMCquantreg.R
@@ -0,0 +1,84 @@
+"MCMCquantreg" <-
+ function(formula, p=0.5, data=NULL, burnin = 1000, mcmc = 10000,
+ thin=1, verbose = 0, seed = NA, beta.start = NA,
+ b0 = 0, B0 = 0, c0 = 0.001, d0 = 0.001,
+ ...) {
+
+ ## checks
+ check.offset(list(...))
+ check.mcmc.parameters(burnin, mcmc, thin)
+
+ cl <- match.call()
+ if (p<=0 || p>=1){
+ stop("p must be in (0,1).\n Please respecify and call again.\n")
+ }
+
+ ## seeds
+ seeds <- form.seeds(seed)
+ lecuyer <- seeds[[1]]
+ seed.array <- seeds[[2]]
+ lecuyer.stream <- seeds[[3]]
+
+ ## form response and model matrices
+ holder <- parse.formula(formula, data=data)
+ Y <- holder[[1]]
+ X <- holder[[2]]
+ xnames <- holder[[3]]
+ K <- ncol(X) # number of covariates
+
+ ## starting values and priors
+ ols.fit <- lm(formula)
+ defaults <- matrix(coef(ols.fit),K,1)
+ defaults[1] <- defaults[1]+summary(ols.fit)$sigma*qnorm(p)
+ beta.start <- coef.start(beta.start, K, formula, family=gaussian, data, defaults=defaults)
+ mvn.prior <- form.mvn.prior(b0, B0, K)
+ b0 <- mvn.prior[[1]]
+ B0 <- mvn.prior[[2]]
+ check.ig.prior(c0, d0)
+
+
+ B0.eigenvalues <- eigen(B0)$values
+ if (min(B0.eigenvalues) < 0){
+ stop("B0 is not positive semi-definite.\nPlease respecify and call again.\n")
+ }
+
+
+
+ ## define holder for posterior sample
+ sample <- matrix(data=0, mcmc/thin, K+1)
+ posterior <- NULL
+
+ if (p==0.5) {
+ ## call C++/Scythe function "MCMCmedreg" to draw samples
+ auto.Scythe.call(output.object="posterior", cc.fun.name="MCMCmedreg",
+ sample.nonconst=sample, Y=Y, X=X, burnin=as.integer(burnin),
+ mcmc=as.integer(mcmc), thin=as.integer(thin),
+ lecuyer=as.integer(lecuyer),
+ seedarray=as.integer(seed.array),
+ lecuyerstream=as.integer(lecuyer.stream),
+ verbose=as.integer(verbose), betastart=beta.start,
+ b0=b0, B0=B0, c0=as.double(c0), d0=as.double(d0), package="MCMCpack")
+}
+
+ else{
+
+ ## call C++/Scythe function "MCMCquantreg" to draw samples
+ auto.Scythe.call(output.object="posterior", cc.fun.name="MCMCquantreg",
+ sample.nonconst=sample, p=as.double(p), Y=Y, X=X,
+ burnin=as.integer(burnin),
+ mcmc=as.integer(mcmc), thin=as.integer(thin),
+ lecuyer=as.integer(lecuyer),
+ seedarray=as.integer(seed.array),
+ lecuyerstream=as.integer(lecuyer.stream),
+ verbose=as.integer(verbose), betastart=beta.start,
+ b0=b0, B0=B0, c0=as.double(c0), d0=as.double(d0), package="MCMCpack")
+}
+
+## pull together matrix and build MCMC object to return
+ output <- form.mcmc.object(posterior,
+ names=c(xnames,"sigma"),
+ title="MCMCquantreg Posterior Sample",
+ y=Y, call=cl)
+
+ return(output)
+}
diff --git a/man/MCMCirtHier1d.Rd b/man/MCMCirtHier1d.Rd
index 048b30f..7a42fbf 100644
--- a/man/MCMCirtHier1d.Rd
+++ b/man/MCMCirtHier1d.Rd
@@ -52,7 +52,7 @@ MCMCirtHier1d(datamatrix, Xjdata,
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
+ 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.}
@@ -286,7 +286,7 @@ Xjdata <- data.frame(presparty= c(1,1,0,1,1,1,1,0,0),
## the variance of the (unidentified) latent parameter alpha.
scMiss <- SupremeCourt
-scMiss[matrix(as.logical(rbinom( nrow(SupremeCourt)*ncol(SupremeCourt), 1, .1)), dim(SupremeCourt))] <- NA
+scMiss[matrix(as.logical(rbinom(nrow(SupremeCourt)*ncol(SupremeCourt), 1, .1)), dim(SupremeCourt))] <- NA
posterior1.miss <- MCMCirtHier1d(t(scMiss),
burnin=80000, mcmc=10000, thin=20,
diff --git a/man/MCMCquantreg.Rd b/man/MCMCquantreg.Rd
new file mode 100644
index 0000000..ca7a8cf
--- /dev/null
+++ b/man/MCMCquantreg.Rd
@@ -0,0 +1,156 @@
+\name{MCMCquantreg}
+
+\alias{MCMCquantreg}
+
+\title{Bayesian quantile regression using Gibbs sampling}
+
+\description{
+ This function fits quantile regression models under Bayesian inference.
+ The function samples from the posterior distribution using Gibbs sampling with data augmentation.
+ A multivariate normal prior is assumed for \eqn{\beta}{beta} and an inverse Gamma prior
+ is assumed for \eqn{\sigma}{sigma}. The user supplies the prior parameters.
+ A sample of the posterior distribution is returned as an mcmc object,
+ which can then be analysed by functions in the coda package.
+}
+
+\usage{
+MCMCquantreg(formula, p=0.5, data = NULL, burnin = 1000,
+ mcmc = 10000, thin = 1, verbose = 0, seed = NA,
+ beta.start = NA, b0 = 0, B0 = 0, c0 = 0.001, d0 = 0.001, ...)
+}
+
+\arguments{
+
+ \item{formula}{ Model formula. }
+
+ \item{p}{The quantile of interest. Must be between 0 and 1. The default value of 0.5 corresponds to median regression.}
+
+
+ \item{data}{ Data frame. }
+
+ \item{burnin}{ The number of burn-in iterations for the sampler. }
+
+ \item{mcmc}{ The number of MCMC iterations after burnin. }
+
+ \item{thin}{ The thinning interval used in the simulation. The number of
+ MCMC iterations must be divisible by this value. }
+
+ \item{verbose}{ A switch which determines whether or not the progress of
+ the sampler is printed to the screen. If \code{verbose} is greater
+ than 0 the iteration number and the most recently sampled values of
+ \eqn{\beta}{beta} and \eqn{\sigma}{sigma} are printed to
+ the screen every \code{verbose}th iteration. }
+
+ \item{seed}{ The seed for the random number generator. If NA, the Mersenne
+ 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{beta.start}{ The starting values for \eqn{\beta}{beta}.
+ This can either be a scalar or a
+ column vector with dimension equal to the dimension of \eqn{\beta}{beta}.
+ The default value of NA will use the OLS
+ estimate \eqn{\hat{\beta}}{beta^hat} with
+ \eqn{\hat{\sigma}\Phi^{-1}(p)}{sigma^hat*Phi^(-1)(p)} added on to the first element of \eqn{\hat{\beta}}{beta^hat} as the starting value.
+ (\eqn{\hat{\sigma}^2}{(sigma^hat)^2} denotes the usual unbiased estimator of
+ \eqn{\sigma^2}{sigma^2} under ordinary mean regression and
+ \eqn{\Phi^{-1}(p)}{Phi^(-1)(p)} denotes the inverse of the
+ cumulative density function of the standard normal distribution.)
+ Note that the default value assume that an intercept is included in the model.
+ If a scalar is given, that value will serve as the starting value
+ for all \eqn{\beta}{beta}. }
+
+ \item{b0}{ The prior mean of \eqn{\beta}{beta}. This can either be a
+ scalar or a
+ column vector with dimension equal to the dimension of \eqn{\beta}{beta}. If this
+ takes a scalar value, then that value will serve as the prior
+ mean for all \eqn{\beta}{beta}. }
+
+ \item{B0}{ The prior precision of \eqn{\beta}{beta}. This can either be a
+ scalar or a square matrix with dimensions equal to the number of betas.
+ If this
+ takes a scalar value, then that value times an identity matrix serves
+ as the prior precision of \eqn{\beta}{beta}. Default value of 0 is equivalent to
+ an improper uniform prior for \eqn{\beta}{beta}. }
+
+ \item{c0}{ \eqn{c_0/2}{c0/2} is the shape parameter for the inverse
+ Gamma prior on \eqn{\sigma}{sigma}. }
+
+ \item{d0}{ \eqn{d_0/2}{d0/2} is the scale parameter for the
+ inverse Gamma prior on \eqn{\sigma}{sigma}.}
+
+ \item{\dots}{ further arguments to be passed }
+}
+
+\details{
+ \code{MCMCquantreg} simulates from the posterior distribution using
+ a partially collapsed Gibbs sampler with data augmentation (see \url{http://people.brunel.ac.uk/~mastkky/}).
+ This involves a multivariate Normal draw for \eqn{\beta}{beta}, and an
+ inverse Gamma draw for \eqn{\sigma}{sigma}.
+ The augmented data are drawn conditionally from the inverse Gaussian distribution. The simulation
+ is carried out in compiled C++ code to maximise efficiency. Please consult
+ the coda documentation for a comprehensive list of functions that can be
+ used to analyse the posterior sample.
+
+ The model takes the following form:
+ \deqn{y_i = x_i ' \beta + \varepsilon_i}{y_i = x_i'beta + epsilon_i}
+ The errors are assumed to have an Asymmetric Laplace distribution with pth quantile equal to zero
+ and shape parameter equal to \eqn{\sigma}{sigma}:
+ \deqn{\varepsilon_i(p) \sim \mathcal{AL}(0, \sigma, p)}{epsilon_i(p) ~ AL(0,
+ sigma, p),}
+ where \eqn{\beta}{beta} and \eqn{\sigma}{sigma} depend on p.
+ We assume standard, semi-conjugate priors:
+ \deqn{\beta \sim \mathcal{N}(b_0,B_0^{-1})}{beta ~ N(b0,B0^(-1)),}
+ and
+ \deqn{\sigma^{-1} \sim \mathcal{G}amma(c_0/2, d_0/2)}{sigma^(-1) ~
+ Gamma(c0/2, d0/2),}
+ where \eqn{\beta}{beta} and \eqn{\sigma^{-1}}{sigma^(-1)} are assumed
+ \emph{a priori} independent of each other and the parameters associated with the other quantiles. Only starting values for
+ \eqn{\beta}{beta} are allowed for this sampler.
+}
+
+\value{
+ An mcmc object that contains the posterior sample. This
+ object can be summarised by functions provided by the coda package.
+}
+
+\references{ Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007.
+ \emph{Scythe Statistical Library 1.2.} \url{http://scythe.wustl.edu}.
+
+ Craig Reed and Keming Yu. 2009. ``A Partially Collapsed Gibbs Sampler for Bayesian Quantile Regression." Technical Report.
+
+ Keming Yu and Jin Zhang. 2005. ``A Three Parameter Asymmetric Laplace Distribution and it's extensions."
+\emph{Communications in Statistics - Theory and Methods}, 34, 1867-1879.
+
+ Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2002.
+ \emph{Output Analysis and Diagnostics for MCMC (CODA)}.
+ \url{http://www-fis.iarc.fr/coda/}.}
+
+\author{Craig Reed}
+
+\seealso{ \code{\link[MCMCpack]{MCMCregress}}, \code{\link[coda]{plot.mcmc}},
+ \code{\link[coda]{summary.mcmc}}, \code{\link[stats]{lm}}, \code{\link[quantreg]{rq}} }
+
+\examples{
+\dontrun{
+
+x<-rep(1:10,5)
+y<-rnorm(50,mean=x)
+posterior_50 <- MCMCquantreg(y~x)
+posterior_95 <- MCMCquantreg(y~x, p=0.95, verbose=10000,
+ mcmc=50000, thin=10, seed=2)
+plot(posterior_50)
+plot(posterior_95)
+raftery.diag(posterior_50)
+autocorr.plot(posterior_95)
+summary(posterior_50)
+summary(posterior_95)
+}
+}
+
+\keyword{models}
diff --git a/src/MCMCfcds.h b/src/MCMCfcds.h
index e77e9c8..b55ea66 100644
--- a/src/MCMCfcds.h
+++ b/src/MCMCfcds.h
@@ -19,14 +19,14 @@
//
// KQ 6/10/2004
// DBP 7/01/2007 [ported to scythe 1.0.x (partial)]
+// ADM 7/28/2009 [added some functions from Craig Reed for quantile
+// regression]
//
// 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 MCMCFCDS_H
#define MCMCFCDS_H
@@ -39,8 +39,8 @@
#include "distributions.h"
#include <iostream>
-using namespace std;
+using namespace std;
using namespace scythe;
// linear regression with Gaussian errors beta draw
@@ -87,6 +87,172 @@ NormIGregress_sigma2_draw (const Matrix <> &X, const Matrix <> &Y,
return stream.rigamma (c_post, d_post);
}
+////////////// New Additions ////////////////////////
+
+// linear regression with Laplace errors beta draw
+// (multivariate Normal prior)
+// regression model is y = X * beta + epsilon, epsilon ~ Laplace(0,sigma)
+// b0 is the prior mean of beta
+// B0 is the prior precision (the inverse variance) of beta
+template <typename RNGTYPE>
+Matrix<double>
+LaplaceNormregress_beta_draw (const Matrix<>& X, const Matrix<>& Y, const Matrix<>& weights,
+ const Matrix<>& b0, const Matrix<>& B0, double sigma,
+ rng<RNGTYPE>& stream)
+{
+
+ const unsigned int k = X.cols();
+ const unsigned int n_obs = X.rows();
+ const double one_over_two_sigma = 1.0/(2.0*sigma);
+ Matrix<> XtwX(k,k,false);
+ Matrix<> XtwY(k,1,false);
+ double temp_x = 0.0;
+ double temp_y = 0.0;
+
+//Calculate XtwY, where w denotes a diagonal matrix with the augmented data (weights) on the diagonal
+ for (unsigned int i=0; i<k; ++i){
+ for (unsigned int m=0; m<n_obs; ++m){
+ temp_y = temp_y + weights(m)*X(m,i)*Y(m);
+ }
+ XtwY(i) = temp_y;
+ temp_y = 0.0;
+ }
+//Calculate XtwX
+ for (unsigned int i=0; i<k; ++i){
+ for (unsigned int j=0; j<(i+1); ++j){
+ for (unsigned int m=0; m<n_obs; ++m){
+ temp_x = temp_x + weights(m)*X(m,i)*X(m,j);
+ }
+ XtwX(i,j) = temp_x;
+ XtwX(j,i) = temp_x;
+ temp_x = 0.0;
+ }
+ }
+
+ const Matrix<> var_matrix_beta = invpd(B0+one_over_two_sigma*XtwX);
+ const Matrix<> C = cholesky(var_matrix_beta);
+ const Matrix<> betahat = var_matrix_beta*gaxpy(B0,b0,one_over_two_sigma*XtwY);
+
+ return( gaxpy(C, stream.rnorm(k,1, 0, 1), betahat) );
+}
+
+
+// linear regression with Laplace errors sigma draw
+// (inverse-Gamma prior)
+// regression model is y = X * beta + epsilon, epsilon ~ Laplace(0,sigma)
+// c0/2 is the prior shape parameter for sigma
+// d0/2 is the prior scale parameter for sigma
+template <typename RNGTYPE>
+double
+LaplaceIGammaregress_sigma_draw (const Matrix <> &abse, double c0, double d0,
+ rng<RNGTYPE>& stream)
+
+{
+ const double c_post = 0.5*c0 + abse.rows();
+ const double d_post = 0.5*(d0 + sum(abse));
+
+ return stream.rigamma (c_post, d_post);
+
+}
+
+// linear regression with Asymmetric Laplace errors beta draw
+// (multivariate Normal prior)
+// regression model is y = X * beta + epsilon, epsilon ~ ALaplace(0,sigma,p)
+// b0 is the prior mean of beta
+// B0 is the prior precision (the inverse variance) of beta
+template <typename RNGTYPE>
+Matrix<double>
+ALaplaceNormregress_beta_draw (double p, const Matrix<>& X, const Matrix<>& Y, const Matrix<>& weights,
+ const Matrix<>& b0, const Matrix<>& B0, double sigma,
+ rng<RNGTYPE>& stream)
+{
+
+ const unsigned int k = X.cols();
+ const unsigned int n_obs = X.rows();
+ const double one_over_two_sigma = 1.0/(2.0*sigma);
+ const Matrix<> U = Y - (1.0-2.0*p)*pow(weights,-1.0);
+ Matrix<> XtwX(k,k,false);
+ Matrix<> XtwU(k,1,false);
+ double temp_x = 0.0;
+ double temp_u = 0.0;
+
+//Calculate XtwU where w denotes a diagonal matrix with the augmented data (weights) on the diagonal
+ for (unsigned int i=0; i<k; ++i){
+ for (unsigned int m=0; m<n_obs; ++m){
+ temp_u = temp_u + weights(m)*X(m,i)*U(m);
+ }
+ XtwU(i) = temp_u;
+ temp_u = 0.0;
+ }
+//Calculate XtwX
+ for (unsigned int i=0; i<k; ++i){
+ for (unsigned int j=0; j<(i+1); ++j){
+ for (unsigned int m=0; m<n_obs; ++m){
+ temp_x = temp_x + weights(m)*X(m,i)*X(m,j);
+ }
+ XtwX(i,j) = temp_x;
+ XtwX(j,i) = temp_x;
+ temp_x = 0.0;
+ }
+ }
+
+ const Matrix<> var_matrix_beta = invpd(B0+one_over_two_sigma*XtwX);
+ const Matrix<> C = cholesky(var_matrix_beta);
+ const Matrix<> betahat = var_matrix_beta*gaxpy(B0,b0,one_over_two_sigma*XtwU);
+
+ return( gaxpy(C, stream.rnorm(k,1, 0, 1), betahat) );
+}
+
+// linear regression with Asymmetric Laplace errors sigma draw
+// (inverse-Gamma prior)
+// regression model is y = X * beta + epsilon, epsilon ~ ALaplace(0,sigma, p)
+// c0/2 is the prior shape parameter for sigma
+// d0/2 is the prior scale parameter for sigma
+template <typename RNGTYPE>
+double
+ALaplaceIGammaregress_sigma_draw (double p, const Matrix<> &e, const Matrix <> &abse, double c0, double d0,
+ rng<RNGTYPE>& stream)
+
+{
+ const double c_post = 0.5*c0 + abse.rows();
+ const double d_post = 0.5*(d0 + sum(abse)+(2.0*p-1.0)*sum(e));
+
+ return stream.rigamma (c_post, d_post);
+
+}
+
+// This function draws from the full conditional distribution of the latent random variables (weights) under quantile regression (including median regression) and returns a column vector of those weights.
+
+template <typename RNGTYPE>
+Matrix<double>
+ALaplaceIGaussregress_weights_draw (const Matrix <> &abse, double sigma,
+ rng<RNGTYPE>& stream)
+
+{
+ const double lambda = 1.0/(2.0*sigma);
+ const Matrix<double> nu_params = pow(abse,-1.0);
+ Matrix<> w(abse);
+ unsigned int n_obs = abse.rows();
+
+ // The inverse Gaussian distribution
+
+ for (unsigned int i=0; i<n_obs; ++i){
+ double chisq = stream.rchisq(1);
+ double nu = nu_params(i);
+ double smallroot = nu/(2.0*lambda)*(nu*chisq+2.0*lambda-std::sqrt(nu*nu*chisq*chisq+4.0*nu*lambda*chisq));
+ unsigned int q = stream.rbern(nu/(nu+smallroot));
+ if (q == 1){
+ w(i) = smallroot;
+ }
+ else{
+ w(i) = nu*nu/smallroot;
+ }
+ }
+ return w;
+}
+
+////////////////////////////////////////////////////
+
// update latent data for standard item response models
// only works for 1 dimensional case
template <typename RNGTYPE>
diff --git a/src/MCMCmedreg.cc b/src/MCMCmedreg.cc
new file mode 100644
index 0000000..81c52fb
--- /dev/null
+++ b/src/MCMCmedreg.cc
@@ -0,0 +1,138 @@
+// MCMCmedreg.cc is a program that simulates draws from the
+// posterior distribution of a linear regression model with
+// errors that come from a Laplace (double exponential) distribution.
+//
+// The initial version of this file was generated by the
+// auto.Scythe.call() function in the MCMCpack R package
+// written by:
+//
+// Andrew D. Martin
+// Dept. of Political Science
+// Washington University in St. Louis
+// admartin at wustl.edu
+//
+// Kevin M. Quinn
+// Dept. of Government
+// Harvard University
+// kevin_quinn at harvard.edu
+//
+// This software is distributed under the terms of the GNU GENERAL
+// PUBLIC LICENSE Version 2, June 1991. See the package LICENSE
+// file for more information.
+//
+// Copyright (C) 2009 Andrew D. Martin and Kevin M. Quinn
+//
+// This file was initially generated on Wed Apr 1 11:39:12 2009
+// The function was rewritten by:
+//
+// Craig Reed
+// Department of Mathematical Sciences
+// Brunel University
+// craig.reed at brunel.ac.uk
+//
+// CR 7/4/09 [Rewritten using templates]
+
+#ifndef MCMCMEDREG_CC
+#define MCMCMEDREG_CC
+
+#include "MCMCrng.h"
+#include "MCMCfcds.h"
+#include "matrix.h"
+#include "distributions.h"
+#include "stat.h"
+#include "la.h"
+#include "ide.h"
+#include "smath.h"
+#include "rng.h"
+
+#include <R.h> // needed to use Rprintf()
+#include <R_ext/Utils.h> // needed to allow user interrupts
+
+using namespace std;
+using namespace scythe;
+
+/* MCMCmedreg implementation. Takes Matrix<> reference which it
+ * fills with the posterior.
+ */
+template <typename RNGTYPE>
+void MCMCmedreg_impl (rng<RNGTYPE>& stream, const Matrix<>& Y,
+ const Matrix<>& X, Matrix<>& beta, const Matrix<>& b0,
+ const Matrix<>& B0, double c0, double d0,
+ unsigned int burnin, unsigned int mcmc, unsigned int thin,
+ unsigned int verbose,
+ Matrix<>& result)
+{
+ // define constants
+ const unsigned int tot_iter = burnin + mcmc; //total iterations
+ const unsigned int nstore = mcmc / thin; // number of draws to store
+ const unsigned int k = X.cols ();
+ const unsigned int n_obs = X.rows();
+
+ // storage matrices
+ Matrix<> betamatrix (k, nstore);
+ Matrix<> sigmamatrix (1, nstore);
+
+ // Gibbs sampler
+ unsigned int count = 0;
+ for (unsigned int iter = 0; iter < tot_iter; ++iter) {
+ Matrix<> abse = fabs(gaxpy(X, (-1*beta), Y));
+ double sigma = LaplaceIGammaregress_sigma_draw (abse, c0, d0, stream);
+ Matrix<> weights = ALaplaceIGaussregress_weights_draw (abse, sigma, stream);
+ beta = LaplaceNormregress_beta_draw (X, Y, weights, b0, B0, sigma, stream);
+
+ // store draws in storage matrix (or matrices)
+ if (iter >= burnin && (iter % thin == 0)) {
+ sigmamatrix (0, count) = sigma;
+ betamatrix(_, count) = beta;
+ ++count;
+ }
+
+ // print output to stdout
+ if(verbose > 0 && iter % verbose == 0) {
+ Rprintf("\n\nMCMCmedreg iteration %i of %i \n",
+ (iter+1), tot_iter);
+ Rprintf("beta = \n");
+ for (unsigned int j=0; j<k; ++j)
+ Rprintf("%10.5f\n", beta(j));
+ Rprintf("sigma = %10.5f\n", sigma);
+ }
+
+ R_CheckUserInterrupt(); // allow user interrupts
+
+ } // end MCMC loop
+
+ result = cbind(t(betamatrix), t(sigmamatrix));
+ } // end MCMCmedreg_impl
+
+extern "C" {
+ void MCMCmedreg(double *sampledata, const int *samplerow,
+ const int *samplecol, const double *Ydata, const int *Yrow,
+ const int *Ycol, const double *Xdata, const int *Xrow,
+ const int *Xcol, const int *burnin, const int *mcmc,
+ const int *thin, const int *uselecuyer, const int *seedarray,
+ const int *lecuyerstream, const int *verbose,
+ const double *betastartdata, const int *betastartrow,
+ const int *betastartcol, const double *b0data,
+ const int *b0row, const int *b0col,
+ const double *B0data, const int *B0row,
+ const int *B0col, const double *c0, const double *d0)
+ {
+ // pull together Matrix objects
+ Matrix<> Y(*Yrow, *Ycol, Ydata);
+ Matrix<> X(*Xrow, *Xcol, Xdata);
+ Matrix<> betastart(*betastartrow, *betastartcol, betastartdata);
+ Matrix<> b0(*b0row, *b0col, b0data);
+ Matrix<> B0(*B0row, *B0col, B0data);
+
+ Matrix<> storagematrix;
+ MCMCPACK_PASSRNG2MODEL(MCMCmedreg_impl, Y, X, betastart, b0, B0,
+ *c0, *d0, *burnin, *mcmc, *thin, *verbose,
+ storagematrix);
+
+ const unsigned int size = *samplerow * *samplecol;
+ for (unsigned int i = 0; i < size; ++i)
+ sampledata[i] = storagematrix(i);
+ }
+}
+
+#endif
diff --git a/src/MCMCquantreg.cc b/src/MCMCquantreg.cc
new file mode 100644
index 0000000..1b81884
--- /dev/null
+++ b/src/MCMCquantreg.cc
@@ -0,0 +1,140 @@
+// MCMCquantreg.cc is a function that draws from the posterior
+// distribution of a linear regression model with errors that
+// come from an Asymmetric Laplace distribution having pth quantile equal to zero.
+//
+// The initial version of this file was generated by the
+// auto.Scythe.call() function in the MCMCpack R package
+// written by:
+//
+// Andrew D. Martin
+// Dept. of Political Science
+// Washington University in St. Louis
+// admartin at wustl.edu
+//
+// Kevin M. Quinn
+// Dept. of Government
+// Harvard University
+// kevin_quinn at harvard.edu
+//
+// This software is distributed under the terms of the GNU GENERAL
+// PUBLIC LICENSE Version 2, June 1991. See the package LICENSE
+// file for more information.
+//
+// Copyright (C) 2009 Andrew D. Martin and Kevin M. Quinn
+//
+// This file was initially generated on Wed Apr 1 11:39:12 2009
+//
+// The function was rewritten by:
+//
+// Craig Reed
+// Department of Mathematical Sciences
+// Brunel University
+// craig.reed at brunel.ac.uk
+//
+// CR 9/4/09 [Rewritten function using templates]
+
+#ifndef MCMCQUANTREG_CC
+#define MCMCQUANTREG_CC
+
+#include "MCMCrng.h"
+#include "MCMCfcds.h"
+#include "matrix.h"
+#include "distributions.h"
+#include "stat.h"
+#include "la.h"
+#include "ide.h"
+#include "smath.h"
+#include "rng.h"
+
+#include <R.h> // needed to use Rprintf()
+#include <R_ext/Utils.h> // needed to allow user interrupts
+
+using namespace std;
+using namespace scythe;
+
+/* MCMCquantreg implementation. Takes Matrix<> reference which it
+ * fills with the posterior.
+ */
+template <typename RNGTYPE>
+void MCMCquantreg_impl (rng<RNGTYPE>& stream, double p, const Matrix<>& Y,
+ const Matrix<>& X, Matrix<>& beta, const Matrix<>& b0,
+ const Matrix<>& B0, double c0, double d0,
+ unsigned int burnin, unsigned int mcmc, unsigned int thin,
+ unsigned int verbose,
+ Matrix<>& result)
+{
+ // define constants
+ const unsigned int tot_iter = burnin + mcmc; //total iterations
+ const unsigned int nstore = mcmc / thin; // number of draws to store
+ const unsigned int k = X.cols ();
+ const unsigned int n_obs = X.rows();
+
+ // storage matrices
+ Matrix<> betamatrix (k, nstore);
+ Matrix<> sigmamatrix (1, nstore);
+
+ // Gibbs sampler
+ unsigned int count = 0;
+ for (unsigned int iter = 0; iter < tot_iter; ++iter) {
+ Matrix<> e = gaxpy(X, (-1*beta), Y);
+ Matrix<> abse = fabs(e);
+ double sigma = ALaplaceIGammaregress_sigma_draw (p, e, abse, c0, d0, stream);
+ Matrix<> weights = ALaplaceIGaussregress_weights_draw (abse, sigma, stream);
+ beta = ALaplaceNormregress_beta_draw (p, X, Y, weights, b0, B0, sigma, stream);
+
+ // store draws in storage matrix (or matrices)
+ if (iter >= burnin && (iter % thin == 0)) {
+ sigmamatrix (0, count) = sigma;
+ betamatrix(_, count) = beta;
+ ++count;
+ }
+
+ // print output to stdout
+ if(verbose > 0 && iter % verbose == 0) {
+ Rprintf("\n\nMCMCquantreg iteration %i of %i \n",
+ (iter+1), tot_iter);
+ Rprintf("beta = \n");
+ for (unsigned int j=0; j<k; ++j)
+ Rprintf("%10.5f\n", beta(j));
+ Rprintf("sigma = %10.5f\n", sigma);
+ }
+
+ R_CheckUserInterrupt(); // allow user interrupts
+
+ } // end MCMC loop
+
+ result = cbind(t(betamatrix), t(sigmamatrix));
+ } // end MCMCquantreg_impl
+
+extern "C" {
+ void MCMCquantreg(double *sampledata, const int *samplerow,
+ const int *samplecol, const double *p, const double *Ydata, const int *Yrow,
+ const int *Ycol, const double *Xdata, const int *Xrow,
+ const int *Xcol, const int *burnin, const int *mcmc,
+ const int *thin, const int *uselecuyer, const int *seedarray,
+ const int *lecuyerstream, const int *verbose,
+ const double *betastartdata, const int *betastartrow,
+ const int *betastartcol, const double *b0data,
+ const int *b0row, const int *b0col,
+ const double *B0data, const int *B0row,
+ const int *B0col, const double *c0, const double *d0)
+ {
+ // pull together Matrix objects
+ Matrix<> Y(*Yrow, *Ycol, Ydata);
+ Matrix<> X(*Xrow, *Xcol, Xdata);
+ Matrix<> betastart(*betastartrow, *betastartcol, betastartdata);
+ Matrix<> b0(*b0row, *b0col, b0data);
+ Matrix<> B0(*B0row, *B0col, B0data);
+
+ Matrix<> storagematrix;
+ MCMCPACK_PASSRNG2MODEL(MCMCquantreg_impl, *p, Y, X, betastart, b0, B0,
+ *c0, *d0, *burnin, *mcmc, *thin, *verbose,
+ storagematrix);
+
+ const unsigned int size = *samplerow * *samplecol;
+ for (unsigned int i = 0; i < size; ++i)
+ sampledata[i] = storagematrix(i);
+ }
+}
+
+#endif
--
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