org.knowceans.arms
Class ArmSampler

java.lang.Object
  extended by org.knowceans.arms.ArmSampler
Direct Known Subclasses:
GammaArms, GaussArms, GmmArms, IgmmGibbsSampler.AlphaArms, IgmmGibbsSampler.BetaArms, InvGammaArms

public abstract class ArmSampler
extends java.lang.Object

ArmSampler implements an adaptive rejection Metropolis sampler (ARMS) that can sample from virtually any univariate distribution. The method performs best if a log-concave density is being sampled from, but it also works for other densities, for which an additional Metropolis step is inserted. The greater the difference to log-concave shape, the more Metropolis rejections must be expected.

This implementation is a port of the original c / fortran implementation by Wally Gilks available at http://www.mrc-bsu.cam.ac.uk/BSUsite/Research/ars.shtml.

Please acknowledge this work by referencing the relevant scientific literature and program code (Web: http://www.arbylon.net/projects).

References:

Gilks, W. R. (1992) Derivative-free adaptive rejection sampling for Gibbs sampling. Bayesian Statistics 4, (eds. Bernardo, J., Berger, J., Dawid, A. P., and Smith, A. F. M.) Oxford University Press.

Gilks, W. R., Best, N. G. and Tan, K. K. C. (1995) Adaptive rejection Metropolis sampling. Applied Statistics, 44, 455-472.

Gilks, W. R. and Wild, P. (1992) Adaptive rejection sampling for Gibbs sampling. Applied Statistics 41, pp 337-348.

Author:
gregor heinrich

Nested Class Summary
(package private)  class ArmSampler.Envelope
          attributes of the entire rejection envelope
(package private)  class ArmSampler.Metropolis
          for metropolis step
(package private)  class ArmSampler.Point
          a point in the x,y plane
 
Field Summary
static int DEREF
          dereference for ported c pointers
static double EYEPS
          critical relative exp(y) difference
(package private)  java.lang.Object params
           
static double XEPS
          critical relative x-value difference
static double YCEIL
          maximum y avoiding overflow in exp(y)
static double YEPS
          critical y-value difference
 
Constructor Summary
ArmSampler()
          init (nothing to do)
 
Method Summary
(package private)  double area(ArmSampler.Point q)
          To integrate piece of exponentiated envelope to left of POINT q
 double[] arms(java.lang.Object params, double[] xinit, int ninit, double[] xl, double[] xr, double[] convex, int npoint, boolean dometrop, double[] xprev, double[] xsamp, int nsamp, double[] qcent, double[] xcent, int ncent, int[] neval)
          to perform derivative-free adaptive rejection sampling with metropolis step
 double armsSimple(java.lang.Object params, int ninit, double[] xl, double[] xr, boolean dometrop, double[] xprev)
          adaptive rejection metropolis sampling - simplified argument list
(package private)  void cumulate(ArmSampler.Envelope env)
          to exponentiate and integrate envelope
(package private)  double expshift(double y, double y0)
          to exponentiate shifted y without underflow
(package private)  int initial(double[] xinit, int ninit, double xl, double xr, int npoint, ArmSampler.Envelope env, double[] convex, int[] neval, ArmSampler.Metropolis metrop)
          to set up initial envelope
(package private)  void invert(double prob, ArmSampler.Envelope env, ArmSampler.Point p)
          to obtain a point corresponding to a qiven cumulative probability
abstract  double logpdf(double x, java.lang.Object params)
          Abstract function to implement the log pdf.
(package private)  double logshift(double y, double y0)
          inverse of function expshift
(package private)  int meet(ArmSampler.Point q, ArmSampler.Envelope env, ArmSampler.Metropolis metrop)
          To find where two chords intersect
(package private)  double perfunc(ArmSampler.Envelope env, double x)
          to evaluate log density and increment count of evaluations
(package private)  void sample(ArmSampler.Envelope env, ArmSampler.Point p)
          To sample from piecewise exponential envelope
(package private)  int test(ArmSampler.Envelope env, ArmSampler.Point p, ArmSampler.Metropolis metrop)
          to perform rejection, squeezing, and metropolis tests
(package private)  int update(ArmSampler.Envelope env, ArmSampler.Point p, ArmSampler.Metropolis metrop)
          to update envelope to incorporate new point on log density
 
Methods inherited from class java.lang.Object
clone, equals, finalize, getClass, hashCode, notify, notifyAll, toString, wait, wait, wait
 

Field Detail

params

java.lang.Object params

DEREF

public static final int DEREF
dereference for ported c pointers

See Also:
Constant Field Values

XEPS

public static final double XEPS
critical relative x-value difference

See Also:
Constant Field Values

YEPS

public static final double YEPS
critical y-value difference

See Also:
Constant Field Values

EYEPS

public static final double EYEPS
critical relative exp(y) difference

See Also:
Constant Field Values

YCEIL

public static final double YCEIL
maximum y avoiding overflow in exp(y)

See Also:
Constant Field Values
Constructor Detail

ArmSampler

public ArmSampler()
init (nothing to do)

Method Detail

logpdf

public abstract double logpdf(double x,
                              java.lang.Object params)
Abstract function to implement the log pdf.

Parameters:
x -
params -
Returns:

armsSimple

public double armsSimple(java.lang.Object params,
                         int ninit,
                         double[] xl,
                         double[] xr,
                         boolean dometrop,
                         double[] xprev)
                  throws java.lang.Exception
adaptive rejection metropolis sampling - simplified argument list

Parameters:
params[] - parameters of the pdf to be sampled
ninit - : number of starting values to be used
xl[] - : left bound ([] for pointer)
xr[] - : right bound ([] for pointer)
dometrop - : whether metropolis step is required
xprev[] - : current value from markov chain ([] for pointer)
Returns:
: sampled value
Throws:
java.lang.Exception

arms

public double[] arms(java.lang.Object params,
                     double[] xinit,
                     int ninit,
                     double[] xl,
                     double[] xr,
                     double[] convex,
                     int npoint,
                     boolean dometrop,
                     double[] xprev,
                     double[] xsamp,
                     int nsamp,
                     double[] qcent,
                     double[] xcent,
                     int ncent,
                     int[] neval)
              throws java.lang.Exception
to perform derivative-free adaptive rejection sampling with metropolis step

Parameters:
params - parameters of the pdf to be sampled
xinit[] - : starting values for x in ascending order
ninit - : number of starting values supplied
xl[] - : left bound ([] for pointer)
xr[] - : right bound ([] for pointer)
convex[] - : adjustment for convexity ([] for pointer)
npoint - : maximum number of envelope points
dometrop - : whether metropolis step is required
xprev[] - : previous value from markov chain ([] for pointer)
xsamp[] - : to store sampled values
nsamp - : number of sampled values to be obtained
qcent[] - : percentages for envelope centiles
xcent[] - : to store requested centiles
ncent - : number of centiles requested
neval[] - : on exit, the number of function evaluations performed ([] for pointer)
Returns:
xsamp[] : sampled values
Throws:
java.lang.Exception

initial

int initial(double[] xinit,
            int ninit,
            double xl,
            double xr,
            int npoint,
            ArmSampler.Envelope env,
            double[] convex,
            int[] neval,
            ArmSampler.Metropolis metrop)
      throws java.lang.Exception
to set up initial envelope

Parameters:
xinit - : initial x-values
ninit - : number of initial x-values
xl,xr - : lower and upper x-bounds
npoint - : maximum number of POINTs allowed in envelope
env[] - : rejection envelope attributes
convex[] - : adjustment for convexity ([] for pointer)
neval[] - : current number of function evaluations ([] for pointer)
metrop - : for metropolis step
Throws:
java.lang.Exception

sample

void sample(ArmSampler.Envelope env,
            ArmSampler.Point p)
      throws java.lang.Exception
To sample from piecewise exponential envelope

Parameters:
env - : envelope attributes
p - : a working POINT to hold the sampled value
Throws:
java.lang.Exception

invert

void invert(double prob,
            ArmSampler.Envelope env,
            ArmSampler.Point p)
      throws java.lang.Exception
to obtain a point corresponding to a qiven cumulative probability

Parameters:
prob - : cumulative probability under envelope
env - : envelope attributes
p - : a working POINT to hold the sampled value
Throws:
java.lang.Exception

test

int test(ArmSampler.Envelope env,
         ArmSampler.Point p,
         ArmSampler.Metropolis metrop)
   throws java.lang.Exception
to perform rejection, squeezing, and metropolis tests

Parameters:
env - : envelope attributes
p - : point to be tested
metrop - : data required for metropolis step
Throws:
java.lang.Exception

update

int update(ArmSampler.Envelope env,
           ArmSampler.Point p,
           ArmSampler.Metropolis metrop)
     throws java.lang.Exception
to update envelope to incorporate new point on log density

Parameters:
env - : envelope attributes
p - : point to be incorporated
metrop - : for metropolis step
Throws:
java.lang.Exception

cumulate

void cumulate(ArmSampler.Envelope env)
        throws java.lang.Exception
to exponentiate and integrate envelope

Parameters:
env - : envelope attributes
Throws:
java.lang.Exception

meet

int meet(ArmSampler.Point q,
         ArmSampler.Envelope env,
         ArmSampler.Metropolis metrop)
   throws java.lang.Exception
To find where two chords intersect

Parameters:
q - : to store point of intersection
env - : envelope attributes
metrop - : for metropolis step
Throws:
java.lang.Exception

area

double area(ArmSampler.Point q)
      throws java.lang.Exception
To integrate piece of exponentiated envelope to left of POINT q

Parameters:
Point - q
Throws:
java.lang.Exception

expshift

double expshift(double y,
                double y0)
to exponentiate shifted y without underflow


logshift

double logshift(double y,
                double y0)
inverse of function expshift


perfunc

double perfunc(ArmSampler.Envelope env,
               double x)
to evaluate log density and increment count of evaluations

Parameters:
*lpdf - : structure containing pointers to log-density function and data
*env - : envelope attributes
x - : point at which to evaluate log density