org.knowceans.arms
Class ArmSampler

java.lang.Object
  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