init header only rgen

.Rbuildignore

@@ -0,0 +1,4 @@
+^.*\.Rproj$
+^\.Rproj\.user$
+^README\.Rmd$
+^README-.*\.png$

.gitignore

@@ -0,0 +1,3 @@
+.Rproj.user
+.Rhistory
+.RData

DESCRIPTION

@@ -0,0 +1,15 @@
+Package: rgen
+Title: Random Sampling Distribution C++ Routines
+Version: 0.0.0.9000
+Authors@R: c(person("James", "Balamuta", email = "balamut2@illinois.edu", role = c("aut", "cre")),
+             person("Steven", "Culpepper", email = "sculpepp@illinois.edu", role = c("ctb")))
+Description: Provides C++ armadillo headers for popular sampling distributions
+Depends:
+    R (>= 3.3.1)
+License: MIT + file LICENSE
+Encoding: UTF-8
+LazyData: true
+LinkingTo: Rcpp, RcppArmadillo
+Imports:
+    Rcpp
+RoxygenNote: 5.0.1

LICENSE

@@ -0,0 +1,2 @@
+YEAR: 2016
+COPYRIGHT HOLDER: James Joseph Balamuta

NAMESPACE

@@ -0,0 +1,4 @@
+# Generated by roxygen2: do not edit by hand
+
+importFrom(Rcpp,sourceCpp)
+useDynLib(rgen)

R/rgen-package.R

@@ -0,0 +1,3 @@
+#' @useDynLib rgen
+#' @importFrom Rcpp sourceCpp
+"_PACKAGE"

README.Rmd

@@ -0,0 +1,17 @@
+---
+output: github_document
+---
+
+<!-- README.md is generated from README.Rmd. Please edit that file -->
+
+```{r, echo = FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>",
+  fig.path = "README-"
+)
+```
+
+# `rgen` - C++ Headers for Sampling Distributions
+The repository houses a few random sampling routines that connect into *R*'s seed generation to be used with [`RcppArmadillo`](https://github.com/RcppCore/RcppArmadillo).
+

README.md

@@ -1,2 +1,6 @@
-# rgen
-Random generation library for distributions
+
+<!-- README.md is generated from README.Rmd. Please edit that file -->
+`rgen` - C++ Headers for Sampling Distributions
+===============================================
+
+The repository houses a few random sampling routines that connect into *R*'s seed generation to be used with [`RcppArmadillo`](https://github.com/RcppCore/RcppArmadillo).

inst/include/dist/rdirichlet.h

@@ -0,0 +1,24 @@
+#ifndef RGEN_RDIRICHLET_H
+#define RGEN_RDIRICHLET_H
+
+#include <RcppArmadillo.h>
+
+// @title Generate Dirichlet Random Variable
+// @description Sample a Dirichlet random variable.
+// @usage rDirichlet(deltas)
+// @param deltas A \code{vector} of Dirichlet parameters.
+// @return A \code{vector} from a Dirichlet.
+// @author Steven Andrew Culpepper
+arma::vec rDirichlet(const arma::vec& deltas){
+    size_t C = deltas.n_elem;
+    arma::vec Xgamma(C);
+
+    //generating gamma(deltac,1)
+    for(size_t c = 0; c < C; c++){
+        Xgamma(c) = R::rgamma(deltas(c),1.0);
+    }
+
+    return Xgamma/sum(Xgamma);
+}
+
+#endif

inst/include/dist/riwishart.h

@@ -0,0 +1,22 @@
+#ifndef RGEN_RIWISHART_H
+#define RGEN_RIWISHART_H
+
+#include <RcppArmadillo.h>
+#include <rwishart.h>
+
+// @title Generate Random Inverse Wishart Distribution
+// @description Creates a random inverse wishart distribution when given degrees of freedom and a sigma matrix.
+// @param df An \code{int} that represents the degrees of freedom.  (> 0)
+// @param S A \code{matrix} with dimensions m x m that provides Sigma, the covariance matrix.
+// @return A \code{matrix} that is an inverse wishart distribution.
+// @seealso \code{\link{rwishart}}
+// @author James J Balamuta
+// @examples
+// #Call with the following data:
+// riwishart(3, diag(2))
+// [[Rcpp::export]]
+arma::mat riwishart(unsigned int df, const arma::mat& S){
+    return rwishart(df, S.i()).i();
+}
+
+#endif

inst/include/dist/rmatrixnorm.h

@@ -0,0 +1,85 @@
+#ifndef RGEN_RMATRIXNORM_H
+#define RGEN_RMATRIXNORM_H
+
+#include <RcppArmadillo.h>
+
+// @title Sample From Matrix Normal Distribution
+// @description Provides the ability to draw one sample from a Matrix Normal Distribution
+// @param mu        A \code{matrix} with dimensions \eqn{N \times P}{N x P} that contains the mean of the Matrix Normal Distribution
+// @param Sigma_row A \code{matrix} with dimensions \eqn{N \times N}{N x N} that contains the row covariance matrix (positive definite) of the Matrix Normal Distribution
+// @param Sigma_col A \code{matrix} with dimensions \eqn{P \times P}{P x P} that contains the column covariance matrix (positive definite) of the Matrix Normal Distribution
+// @rdname rmatnormal
+// @export
+// @details
+// The implementation is based off of the notation that
+// \deqn{X \sim MN_{N\times P}\left({0,I,I}\right)}{X ~ MN[NxP](0,I,I)}
+// Then,
+// \deqn{Y = M + AXB}
+// so that:
+// \deqn{Y \sim MN_{N\times P}\left({M,U=AA^{T},V=B^{T}B}\right)}{Y ~ MN[NxP](M,U=AA^T,V=B^{T}B)}
+//
+// There are two ways to proceed in this generation. The first is to obtain a cholesky decomposition and
+// the second is to use a symmetric eigen decomposition
+arma::mat rmatnormal_chol(const arma::mat& mu, const arma::mat& Sigma_row, const arma::mat& Sigma_col) {
+
+    // Construct X by repeatively sampling from N(0,1)
+    arma::mat X(Sigma_row.n_rows, Sigma_col.n_cols);
+    X.imbue(norm_rand); // Use R's PNG seed instead of Armadillo
+
+    // N x P + N x N * N x P * P x P
+    return mu + arma::chol(Sigma_row) * X * arma::chol(Sigma_col);
+}
+
+// helper function
+inline arma::mat make_mat(const arma::mat& X){
+    // Eigen values
+    arma::vec eigval;
+    arma::mat eigvec;
+
+    eig_sym(eigval, eigvec, X);
+
+    return eigvec * diagmat(sqrt(eigval)) * eigvec.t();
+}
+
+arma::mat rmatnormal_eigen(const arma::mat& mu, const arma::mat& Sigma_row, const arma::mat& Sigma_col) {
+
+    // Construct X by repeatively sampling from N(0,1)
+    arma::mat X(Sigma_row.n_rows, Sigma_col.n_cols);
+    X.imbue(norm_rand); // Use R's PNG seed instead of Armadillo
+
+    // N x P + N x N * N x P * P x P
+    return mu + make_mat(Sigma_row) * X * make_mat(Sigma_col);
+}
+
+
+/*
+# N == P
+set.seed(234)
+mu = matrix(rnorm(5*5), nrow=5)
+sigma_r = matrix(runif(5*5),nrow=5)
+sigma_c = crossprod(matrix(runif(5*5),nrow=5))
+
+rmatnormal(mu,sigma_r, sigma_c)
+
+# N != P
+set.seed(234)
+mu = matrix(rnorm(4*6), nrow=4,ncol = 6)
+sigma_r = crossprod(matrix(runif(4*4),nrow=4))
+sigma_c = crossprod(matrix(runif(6*6),nrow=6))
+
+rmatnormal(mu,sigma_r, sigma_c)
+
+# Test repeated sampling matches distribution
+set.seed(234)
+mu = matrix(rnorm(5*4), nrow=4,ncol = 5)
+sigma_r = crossprod(matrix(runif(4*4),nrow=4))
+sigma_c = crossprod(matrix(runif(5*5),nrow=5))
+
+# Multiple draws
+large_sample = replicate(20,  rmatnormal(mu,sigma_r, sigma_c))
+
+# Find the mean
+apply(large_sample, c(1,2), mean)
+*/
+
+#endif

inst/include/dist/rmultinomial.h

@@ -0,0 +1,23 @@
+#ifndef RGEN_RMULTINOMIAL_H
+#define RGEN_RMULTINOMIAL_H
+
+#include <RcppArmadillo.h>
+
+// @title Generate Multinomial Random Variable
+// @description Sample a multinomial random variable for given probabilities.
+// @usage rmultinomial(ps)
+// @param ps A \code{vector} for the probability of each category.
+// @return A \code{vector} from a multinomial with probability ps.
+// @author Steven Andrew Culpepper
+double rmultinomial(const arma::vec& ps){
+    unsigned int C = ps.n_elem;
+    double u = R::runif(0,1);
+    arma::vec cps = cumsum(ps);
+    arma::vec Ips = arma::zeros<arma::vec>(C);
+
+    Ips.elem(arma::find(cps < u) ).fill(1.0);
+
+    return sum(Ips);
+}
+
+#endif

inst/include/dist/rmvnorm.h

@@ -0,0 +1,24 @@
+#ifndef RGEN_RMVNORM_H
+#define RGEN_RMVNORM_H
+
+#include <RcppArmadillo.h>
+
+// @title Generate Random Multivariate Normal Distribution
+// @description Creates a random Multivariate Normal when given number of obs, mean, and sigma.
+// @param n  An \code{int}, which gives the number of observations.  (> 0)
+// @param mu A \code{vector} length m that represents the means of the normals.
+// @param S  A \code{matrix} with dimensions m x m that provides Sigma, the covariance matrix.
+// @return A \code{matrix} that is a Multivariate Normal distribution
+// @seealso \code{\link{TwoPLChoicemcmc}} and \code{\link{probitHLM}}
+// @author James J Balamuta
+// @examples
+// #Call with the following data:
+// rmvnorm(2, c(0,0), diag(2))
+arma::mat rmvnorm(unsigned int n, const arma::vec& mu, const arma::mat& S){
+    unsigned int ncols = S.n_cols;
+    arma::mat Y(n, ncols);
+    Y.imbue( norm_rand ) ;
+    return arma::repmat(mu, 1, n).t() + Y * arma::chol(S);
+}
+
+#endif

inst/include/dist/rwishart.h

@@ -0,0 +1,60 @@
+#ifndef RGEN_RWISHART_H
+#define RGEN_RWISHART_H
+
+#include <RcppArmadillo.h>
+
+// @title Generate Random Wishart Distribution
+// @description Creates a random wishart distribution when given degrees of freedom and a sigma matrix.
+// @param df An \code{int}, which gives the degrees of freedom of the Wishart.  (> 0)
+// @param S A \code{matrix} with dimensions m x m that provides Sigma, the covariance matrix.
+// @return A \code{matrix} that is a Wishart distribution, aka the sample covariance matrix of a Multivariate Normal Distribution
+// @seealso \code{\link{riwishart}}
+// @author James J Balamuta
+// @examples
+// #Call with the following data:
+// rwishart(3, diag(2))
+//
+// # Validation
+// set.seed(1337)
+// S = toeplitz((10:1)/10)
+// n = 10000
+// o = array(dim = c(10,10,n))
+// for(i in 1:n){
+// o[,,i] = rwishart(20, S)
+// }
+// mR = apply(o, 1:2, mean)
+// Va = 20*(S^2 + tcrossprod(diag(S)))
+// vR = apply(o, 1:2, var)
+// stopifnot(all.equal(vR, Va, tolerance = 1/16))
+//
+// [[Rcpp::export]]
+arma::mat rwishart(unsigned int df, const arma::mat& S){
+    // Dimension of returned wishart
+    unsigned int m = S.n_rows;
+
+    // Z composition:
+    // sqrt chisqs on diagonal
+    // random normals below diagonal
+    // misc above diagonal
+    arma::mat Z(m,m);
+
+    // Fill the diagonal
+    for(unsigned int i = 0; i < m; i++){
+        Z(i,i) = sqrt(R::rchisq(df-i));
+    }
+
+    // Fill the lower matrix with random guesses
+    for(unsigned int j = 0; j < m; j++){
+        for(unsigned int i = j+1; i < m; i++){
+            Z(i,j) = R::rnorm(0,1);
+        }
+    }
+
+    // Lower triangle * chol decomp
+    arma::mat C = arma::trimatl(Z).t() * arma::chol(S);
+
+    // Return random wishart
+    return C.t()*C;
+}
+
+#endif

inst/include/rgen.h

@@ -0,0 +1,11 @@
+#ifndef RGEN_PACKAGE_H
+#define RGEN_PACKAGE_H
+
+#include <dist/rdirichlet.h>
+#include <dist/rmultinomial.h>
+#include <dist/rwishart.h>
+#include <dist/riwishart.h>
+#include <dist/rmvnorm.h>
+#include <dist/rmatrixnorm.h>
+
+#endif

rgen.Rproj

@@ -0,0 +1,21 @@
+Version: 1.0
+
+RestoreWorkspace: No
+SaveWorkspace: No
+AlwaysSaveHistory: Default
+
+EnableCodeIndexing: Yes
+UseSpacesForTab: Yes
+NumSpacesForTab: 4
+Encoding: UTF-8
+
+RnwWeave: knitr
+LaTeX: pdfLaTeX
+
+AutoAppendNewline: Yes
+StripTrailingWhitespace: Yes
+
+BuildType: Package
+PackageUseDevtools: Yes
+PackageInstallArgs: --no-multiarch --with-keep.source
+PackageRoxygenize: rd,collate,namespace

src/.gitignore

@@ -0,0 +1,3 @@
+*.o
+*.so
+*.dll