correct discrete spatial resids

andeek · Apr 28, 2017 · cf5b87e · cf5b87e
1 parent 34fb4e7
commit cf5b87e
Show file tree

Hide file tree

Showing 6 changed files with 70 additions and 36 deletions.
diff --git a/R/RcppExports.R b/R/RcppExports.R
@@ -125,9 +125,10 @@ binary_two_param_sampler <- function(data, params) {
 #'        which are vectors that contain the horizontal and vertical location of each point in space.
 #'        The input "params" is a list of parameter values. This function returns the inverse cdf at a value between 0 and 1 from the conditional distribution
 #' @param params A list of parameters to be passed to the conditional_density function 
+#' @param discrete If you are dealing with a discrete distribution, the string name of the pmf function. Otherwise, leave as an empty string.
 #' @param ncols A integer of the number of columns in the original grid. Only necessary if dealing with "u" and "v" in cdf function.
 #' @export
-spatial_residuals <- function(data, neighbors, conditional_cdf, params, ncols = 0L) {
-    .Call('conclique_spatial_residuals', PACKAGE = 'conclique', data, neighbors, conditional_cdf, params, ncols)
+spatial_residuals <- function(data, neighbors, conditional_cdf, params, discrete = "", ncols = 0L) {
+    .Call('conclique_spatial_residuals', PACKAGE = 'conclique', data, neighbors, conditional_cdf, params, discrete, ncols)
 }
 
diff --git a/man/spatial_residuals.Rd b/man/spatial_residuals.Rd
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
@@ -86,17 +86,18 @@ BEGIN_RCPP
 END_RCPP
 }
 // spatial_residuals
-arma::vec spatial_residuals(arma::vec data, List neighbors, std::string conditional_cdf, List params, int ncols);
-RcppExport SEXP conclique_spatial_residuals(SEXP dataSEXP, SEXP neighborsSEXP, SEXP conditional_cdfSEXP, SEXP paramsSEXP, SEXP ncolsSEXP) {
+arma::vec spatial_residuals(arma::vec data, List neighbors, std::string conditional_cdf, List params, std::string discrete, int ncols);
+RcppExport SEXP conclique_spatial_residuals(SEXP dataSEXP, SEXP neighborsSEXP, SEXP conditional_cdfSEXP, SEXP paramsSEXP, SEXP discreteSEXP, SEXP ncolsSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
     Rcpp::traits::input_parameter< arma::vec >::type data(dataSEXP);
     Rcpp::traits::input_parameter< List >::type neighbors(neighborsSEXP);
     Rcpp::traits::input_parameter< std::string >::type conditional_cdf(conditional_cdfSEXP);
     Rcpp::traits::input_parameter< List >::type params(paramsSEXP);
+    Rcpp::traits::input_parameter< std::string >::type discrete(discreteSEXP);
     Rcpp::traits::input_parameter< int >::type ncols(ncolsSEXP);
-    rcpp_result_gen = Rcpp::wrap(spatial_residuals(data, neighbors, conditional_cdf, params, ncols));
+    rcpp_result_gen = Rcpp::wrap(spatial_residuals(data, neighbors, conditional_cdf, params, discrete, ncols));
     return rcpp_result_gen;
 END_RCPP
 }
diff --git a/src/cdf.cpp b/src/cdf.cpp
@@ -14,6 +14,10 @@ using namespace Rcpp; using namespace arma;
 //' @export
 // [[Rcpp::export]]
 arma::vec gaussian_single_param_cdf(List data, List params) {
+
+  NumericVector datum = data["data"];
+  vec res(datum.length());
+
   RNGScope scope;
 
   double rho = params["rho"];
@@ -26,10 +30,8 @@ arma::vec gaussian_single_param_cdf(List data, List params) {
   vec sums_vec = sums[0];
   vec nums_vec = nums[0];
 
-  NumericVector datum = data["data"];
-
   vec mean_structure = kappa + eta * (sums_vec - nums_vec * kappa);
-  vec res(datum.length());
+
   for(int i = 0; i < mean_structure.n_elem; ++i) {
     NumericVector prob = pnorm(datum, mean_structure(i), rho);
     res(i) = prob(i);

diff --git a/src/sampler.cpp b/src/sampler.cpp
@@ -15,20 +15,23 @@ using namespace Rcpp; using namespace arma;
 //' @export
 // [[Rcpp::export]]
 arma::vec gaussian_single_param_sampler(List data, List params) {
-  RNGScope scope;
-
-  double rho = params["rho"];
-  double kappa = params["kappa"];
-  double eta = params["eta"];
 
   List sums = data["sums"];
   List nums = data["nums"];
 
   vec sums_vec = sums[0];
   vec nums_vec = nums[0];
 
+  vec res(sums_vec.n_elem);
+
+  RNGScope scope;
+
+  double rho = params["rho"];
+  double kappa = params["kappa"];
+  double eta = params["eta"];
+
   vec mean_structure = kappa + eta * (sums_vec - nums_vec * kappa);
-  vec res(rnorm(mean_structure.n_elem));
+  res = rnorm(mean_structure.n_elem);
   res = res * rho + mean_structure;
 
   return(res);
@@ -38,19 +41,21 @@ arma::vec gaussian_single_param_sampler(List data, List params) {
 //' @export
 // [[Rcpp::export]]
 arma::vec binary_single_param_sampler(List data, List params) {
-  RNGScope scope;
-
-  double kappa = params["kappa"];
-  double eta = params["eta"];
-
   List sums = data["sums"];
   List nums = data["nums"];
 
   vec sums_vec = sums[0];
   vec nums_vec = nums[0];
 
+  vec res(sums_vec.n_elem);
+
+  RNGScope scope;
+
+  double kappa = params["kappa"];
+  double eta = params["eta"];
+
   vec mean_structure = log(kappa) - log(1 - kappa) + eta * (sums_vec - nums_vec * kappa);
-  vec res(mean_structure.n_elem);
+
   int i;
 
   for(i = 0; i < res.n_elem; ++i) {
@@ -63,11 +68,6 @@ arma::vec binary_single_param_sampler(List data, List params) {
 //' @export
 // [[Rcpp::export]]
 arma::vec binary_two_param_sampler(List data, List params) {
-  RNGScope scope;
-
-  double kappa = params["kappa"];
-  double eta_1 = params["eta_1"];
-  double eta_2 = params["eta_2"];
 
   List sums = data["sums"];
   List nums = data["nums"];
@@ -77,8 +77,16 @@ arma::vec binary_two_param_sampler(List data, List params) {
   vec sums_vec_2 = sums[1];
   vec nums_vec_2 = nums[1];
 
+  vec res(sums_vec_1.n_elem);
+
+  RNGScope scope;
+
+  double kappa = params["kappa"];
+  double eta_1 = params["eta_1"];
+  double eta_2 = params["eta_2"];
+
   vec mean_structure = log(kappa) - log(1 - kappa) + eta_1 * (sums_vec_1 - nums_vec_1 * kappa) + eta_2 * (sums_vec_2 - nums_vec_2 * kappa);
-  vec res(mean_structure.n_elem);
+
   int i;
 
   for(i = 0; i < res.n_elem; ++i) {

diff --git a/src/spatial_residuals.cpp b/src/spatial_residuals.cpp
@@ -27,15 +27,22 @@ typedef arma::vec (*cdfPtr)(List, List);
 //'        which are vectors that contain the horizontal and vertical location of each point in space.
 //'        The input "params" is a list of parameter values. This function returns the inverse cdf at a value between 0 and 1 from the conditional distribution
 //' @param params A list of parameters to be passed to the conditional_density function 
+//' @param discrete If you are dealing with a discrete distribution, the string name of the pmf function. Otherwise, leave as an empty string.
 //' @param ncols A integer of the number of columns in the original grid. Only necessary if dealing with "u" and "v" in cdf function.
 //' @export
 // [[Rcpp::export]]
-arma::vec spatial_residuals(arma::vec data, List neighbors, std::string conditional_cdf, List params, int ncols = 0) {
+arma::vec spatial_residuals(arma::vec data, List neighbors, std::string conditional_cdf, List params, std::string discrete = "", int ncols = 0) {
+
   int N = neighbors.length();
+  vec res(N);
+
+  RNGScope rngScope;
+
+
   List sums(N);
   List nums(N);
   int n, m;
-  vec loc_u(N), loc_v(N);
+  vec loc_u(N), loc_v(N), lim(N), A_vec(N);
 
   if(ncols > 0) {
     // location information
@@ -46,14 +53,13 @@ arma::vec spatial_residuals(arma::vec data, List neighbors, std::string conditio
     loc_u = s - row*data.n_cols + 1;
   }
 
-
   bool r_func = true;
   Environment global = Environment::global_env();
 
   std::map<std::string, cdfPtr> cdf_map;
   cdf_map["gaussian_single_param"] = gaussian_single_param_cdf;
 
-  // defined samplers in the package
+  // defined cdf in the package
   std::map<std::string, cdfPtr>::iterator it;
   it = cdf_map.find(conditional_cdf);
   if (it != cdf_map.end())
@@ -86,15 +92,28 @@ arma::vec spatial_residuals(arma::vec data, List neighbors, std::string conditio
     sums_nums_loc.push_back(loc_u, "u");
     sums_nums_loc.push_back(loc_v, "v");
   }
-
+
+
   if(r_func) {
     Function cdf = global[conditional_cdf];
     NumericVector resid = cdf(sums_nums_loc, params);
-    vec res(resid.begin(), resid.length());
-    return(res);
+    res = resid;
   } else {
-    vec res = cdf_map[conditional_cdf](sums_nums_loc, params);
-    return(res);
+    res = cdf_map[conditional_cdf](sums_nums_loc, params);
+  }
+
+  if(discrete.length() > 0) {
+    Function pmf = global[discrete];
+    NumericVector limit = pmf(sums_nums_loc, params);
+    NumericVector A = runif(limit.length());
+
+    vec lim(limit.begin(), limit.length());
+    vec A_vec(A.begin(), A.length());
+    vec one(lim.n_elem); one.ones(); 
+
+    res = (one - A_vec) % res + A_vec % (res - lim);
   }
 
+  return(res);
+
 }