View source: R/sampling_nonlinear.R
sampling_nonlinear | R Documentation |
A Gibbs sampler that draws posterior samples of probability parameters conditional on a (possibly nonlinear) indicator function defining a restricted parameter space that is convex.
sampling_nonlinear( k, options, inside, prior = rep(1, sum(options)), M = 1000, start, burnin = 10, eps = 1e-06, progress = TRUE, cpu = 1 )
k |
vector of observed response frequencies. |
options |
number of observable categories/probabilities for each item
type/multinomial distribution, e.g., |
inside |
an indicator function that takes a vector with probabilities
|
prior |
a vector with two positive numbers defining the shape parameters of the beta prior distributions for each binomial rate parameter. |
M |
number of posterior samples drawn from the encompassing model |
start |
only relevant if |
burnin |
number of burnin samples that are discarded. Can be chosen to be small if the maxmimum-a-posteriori estimate is used as the (default) starting value. |
eps |
precision of the bisection algorithm |
progress |
whether a progress bar should be shown (if |
cpu |
either the number of CPUs used for parallel sampling, or a parallel
cluster (e.g., |
Inequality constraints are defined via an indicator function inside
which returns inside(x)=1
(or 0
) if the vector of free parameters
x
is inside (or outside) the model space. Since the vector x
must include only free (!) parameters, the last probability for each
multinomial must not be used in the function inside(x)
!
Efficiency can be improved greatly if the indicator function is defined as C++ code via the function cppXPtr in the package RcppXPtrUtils (see below for examples). In this case, please keep in mind that indexing in C++ starts with 0,1,2... (not with 1,2,3,... as in R)!
For each parameter, the Gibbs sampler draws a sample from the conditional posterior distribution (a scaled, truncated beta). The conditional truncation boundaries are computed with a bisection algorithm. This requires that the restricted parameteter space defined by the indicator function is convex.
# two binomial success probabilities: x = c(x1, x2) # restriction to a circle: model <- function(x) { (x[1] - .50)^2 + (x[2] - .50)^2 <= .15 } # draw prior samples mcmc <- sampling_nonlinear( k = 0, options = c(2, 2), inside = model, M = 1000 ) head(mcmc) plot(c(mcmc[, 1]), c(mcmc[, 2]), xlim = 0:1, ylim = 0:1) ##### Using a C++ indicator function (much faster) cpp_code <- "SEXP inside(NumericVector x){ return wrap( sum(pow(x-.50, 2)) <= .15);}" # NOTE: Uses Rcpp sugar syntax (vectorized sum & pow) # define function via C++ pointer: model_cpp <- RcppXPtrUtils::cppXPtr(cpp_code) mcmc <- sampling_nonlinear( k = 0, options = c(2, 2), inside = model_cpp ) head(mcmc) plot(c(mcmc[, 1]), c(mcmc[, 2]), xlim = 0:1, ylim = 0:1)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.