From cf5b87e3a9c1ce556fe798f2df5ee5f05b48c75c Mon Sep 17 00:00:00 2001 From: Andee Kaplan Date: Fri, 28 Apr 2017 16:01:39 -0500 Subject: [PATCH] correct discrete spatial resids --- R/RcppExports.R | 5 +++-- man/spatial_residuals.Rd | 5 ++++- src/RcppExports.cpp | 7 ++++--- src/cdf.cpp | 8 ++++--- src/sampler.cpp | 44 +++++++++++++++++++++++---------------- src/spatial_residuals.cpp | 37 ++++++++++++++++++++++++-------- 6 files changed, 70 insertions(+), 36 deletions(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index 75b9405..dbc6f0a 100644 --- 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 index 3901f45..6a0ffc3 100644 --- a/man/spatial_residuals.Rd +++ b/man/spatial_residuals.Rd @@ -4,7 +4,8 @@ \alias{spatial_residuals} \title{Get spatial residuals from data given a neighborhood structure and MRF model.} \usage{ -spatial_residuals(data, neighbors, conditional_cdf, params, ncols = 0L) +spatial_residuals(data, neighbors, conditional_cdf, params, discrete = "", + ncols = 0L) } \arguments{ \item{data}{A vector containing data values for each location in the lattice.} @@ -31,6 +32,8 @@ The input "params" is a list of parameter values. This function returns the inve \item{params}{A list of parameters to be passed to the conditional_density function} +\item{discrete}{If you are dealing with a discrete distribution, the string name of the pmf function. Otherwise, leave as an empty string.} + \item{ncols}{A integer of the number of columns in the original grid. Only necessary if dealing with "u" and "v" in cdf function.} } \description{ diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 57a69d5..6659c8b 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -86,8 +86,8 @@ 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; @@ -95,8 +95,9 @@ BEGIN_RCPP 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 index 2ca51c8..dc57b99 100644 --- 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 index 5edfe06..14c20d6 100644 --- a/src/sampler.cpp +++ b/src/sampler.cpp @@ -15,11 +15,6 @@ 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"]; @@ -27,8 +22,16 @@ arma::vec gaussian_single_param_sampler(List data, List params) { 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 index ed0a77d..4edd09c 100644 --- 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 cdf_map; cdf_map["gaussian_single_param"] = gaussian_single_param_cdf; - // defined samplers in the package + // defined cdf in the package std::map::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); + }