|
| 1 | +/* ARS.C - Procedure for performing Adaptive Rejection Sampling. */ |
| 2 | + |
| 3 | +/* Copyright (c) 1995-2003 by Carl Edward Rasmussen and Radford M. Neal |
| 4 | + * |
| 5 | + * Permission is granted for anyone to copy, use, modify, or distribute this |
| 6 | + * program and accompanying programs and documents for any purpose, provided |
| 7 | + * this copyright notice is retained and prominently displayed, along with |
| 8 | + * a note saying that the original programs are available from Radford Neal's |
| 9 | + * web page, and note is made of any changes made to the programs. The |
| 10 | + * programs and documents are distributed without any warranty, express or |
| 11 | + * implied. As the programs were written for research purposes only, they have |
| 12 | + * not been tested to the degree that would be advisable in any important |
| 13 | + * application. All use of these programs is entirely at the user's own risk. |
| 14 | + */ |
| 15 | + |
| 16 | + |
| 17 | +/* This module implements the "Adaptive Rejection Sampling" scheme due to |
| 18 | + * Gilks and Wild. See "Adaptive Rejection Sampling for Gibbs Sampling", |
| 19 | + * Applied Statistics, vol. 41, no. 2, pp. 337-348 (1992). This is not |
| 20 | + * the most sophisticated possible implementation of the method, however. */ |
| 21 | + |
| 22 | + |
| 23 | +#include <math.h> |
| 24 | +#include <stdio.h> |
| 25 | +#include <stdlib.h> |
| 26 | + |
| 27 | +#include "rand.h" |
| 28 | +#include "ars.h" |
| 29 | + |
| 30 | + |
| 31 | +#define MIN_DROP 0.1 /* minimum drop in log probability for initial points */ |
| 32 | +#define TINY 1e-9 /* minimum tolerance for inserting new points */ |
| 33 | +#define MAX_LIST 100 /* max number of segments in piecewise exponential */ |
| 34 | + |
| 35 | +struct segtype { double a; /* in the log domain, the piece wise... */ |
| 36 | + double b; /* exponential is a piece wise linear: y=ax+b */ |
| 37 | + double x; /* interior point in piece */ |
| 38 | + double xmax; /* upper limit of this piece */ |
| 39 | + double mass; /* the probability mass of this piece */ |
| 40 | + struct segtype *prv; /* ptr to previous piece */ |
| 41 | + struct segtype *nxt; }; /* ptr to next piece */ |
| 42 | + |
| 43 | +static struct segtype *root; |
| 44 | + |
| 45 | +static double rpwed(struct segtype **); |
| 46 | + |
| 47 | + |
| 48 | +/* ADAPTIVE REJECTION SAMPLING PROCEDURE. |
| 49 | + * |
| 50 | + * The ars function returns a single point sampled at random from the univariate |
| 51 | + * distribution specified by the function logp, provided by the user. The |
| 52 | + * logp function takes a point as its first argument, and returns the log |
| 53 | + * of the probability density at that point, plus any arbitrary constant (not |
| 54 | + * depending on the point). The logp function takes two additional arguments, |
| 55 | + * the first a pointer to a place where it must store the derivative of the |
| 56 | + * log probability density, the second a pointer to additional information |
| 57 | + * describing the distribution, which is passed on unchanged from the last |
| 58 | + * argument of ars. |
| 59 | + * |
| 60 | + * The logp function passed MUST be log-concave. It is assumed that any real |
| 61 | + * number is legal as input to logp. |
| 62 | + * |
| 63 | + * The user must also supply an initial guess, "init", and a typical "scale" |
| 64 | + * of variation. It is not essential that these values be very accurate, but |
| 65 | + * performance will generally depend on their accuracy. |
| 66 | + * |
| 67 | + * The ars function first tries to locate points on either side of the mode; |
| 68 | + * the derivative must have the right sign and be non-negligible, and the drop |
| 69 | + * from the max. probability seen must be of at least moderate size to qualify. |
| 70 | + * Then a piece-wise exponential distribution is iteratively improved using |
| 71 | + * knowledge from rejected points, until a sample is accepted. At most MAX_LIST |
| 72 | + * pieces are introduced in the approximation. Before pieces are inserted |
| 73 | + * various checks are made in order to prevent numerical problems. If new points |
| 74 | + * don't qualify, the piece-wise exponential is simply not updated. A warning |
| 75 | + * will be issued (one time only) when 10000 rejections are exceeded */ |
| 76 | + |
| 77 | +double ars( |
| 78 | + double init, /* initial guess */ |
| 79 | + double scale, /* guess for scale of variation in x */ |
| 80 | + double (*logp)(double, double *, void *), /* function to sample from */ |
| 81 | + void *extra /* any extra information to pass to logp() */ |
| 82 | +) |
| 83 | +{ |
| 84 | + struct segtype seg[MAX_LIST], *prv, *cur, *nxt; |
| 85 | + int i, no_seg = 2; |
| 86 | + static int warning = 0; |
| 87 | + double x, max, f, df; |
| 88 | + |
| 89 | + root = &seg[0]; nxt = &seg[1]; |
| 90 | + root->prv = NULL; root->nxt = nxt; nxt->prv = root; nxt->nxt = NULL; |
| 91 | + |
| 92 | + x = init; f = logp(x, &df, extra); /* find point to the left of mode */ |
| 93 | + max = f; |
| 94 | + while (df<TINY || f>max-MIN_DROP) { |
| 95 | + x -= scale+(init-x); |
| 96 | + f = logp(x, &df, extra); |
| 97 | + if (f>max) max = f; |
| 98 | + } |
| 99 | + root->x = x; root->a = df; root->b = f-x*df; |
| 100 | + |
| 101 | + x = init; f = logp(x, &df, extra); /* find point to the right of mode */ |
| 102 | + while (df>-TINY || f>max-MIN_DROP) { |
| 103 | + x += scale+(x-init); |
| 104 | + f = logp(x, &df, extra); |
| 105 | + if (f>max) max = f; |
| 106 | + } |
| 107 | + nxt->x = x; nxt->a = df; nxt->b = f-x*df; |
| 108 | + root->xmax = (nxt->b-root->b)/(root->a-nxt->a); |
| 109 | + |
| 110 | + for (i=0; ; i++) { /* repeat until a point is accepted */ |
| 111 | + |
| 112 | + if (i==10000 && !(warning)) { |
| 113 | + fprintf(stderr, "WARNING: More than 10000 rejections in ars\n"); |
| 114 | + warning = 1; |
| 115 | + } |
| 116 | + |
| 117 | + cur = root; /* find max y-value; needed to avoid numerical problems */ |
| 118 | + max = cur->a*cur->xmax+cur->b; |
| 119 | + while (cur = cur->nxt, cur->nxt) |
| 120 | + if ((x = cur->a*cur->xmax+cur->b) > max) max = x; |
| 121 | + |
| 122 | + cur = root; /* compute masses */ |
| 123 | + cur->mass = exp(cur->a*cur->xmax+cur->b-max)/cur->a; |
| 124 | + while (prv = cur, cur = cur->nxt, cur->nxt) |
| 125 | + cur->mass = (exp(cur->a*cur->xmax+cur->b-max)- |
| 126 | + exp(cur->a*prv->xmax+cur->b-max))/cur->a; |
| 127 | + cur->mass = -exp(cur->a*prv->xmax+cur->b-max)/cur->a; |
| 128 | + |
| 129 | + x = rpwed(&cur); /* this is the new sample */ |
| 130 | + f = logp(x, &df, extra); |
| 131 | + if (rand_uniform() <= exp(f-cur->a*x-cur->b)) return x; /* success! */ |
| 132 | + |
| 133 | +/* Now, insert a new piece in the piece-wise approximation if the situation is |
| 134 | + * appropriate. Eg, if we have enough memory, if the slope at the new x isn't |
| 135 | + * too small (the exponential distribution will degenerate), and if the slope |
| 136 | + * isn't too close to the slope of current, previous and next segment (or |
| 137 | + * which ever of these may exist), since this may cause numerical problems. */ |
| 138 | + |
| 139 | + if (no_seg < MAX_LIST && fabs(df) > TINY && fabs(df-cur->a) > TINY && |
| 140 | + (!(cur->prv) || fabs(df-cur->prv->a) > TINY) && |
| 141 | + (!(cur->nxt) || fabs(df-cur->nxt->a) > TINY)) { |
| 142 | + |
| 143 | + if (x<cur->x) cur = cur->prv; /* now, insert *after* cur */ |
| 144 | + prv = cur; cur = &seg[no_seg++]; cur->prv = prv; |
| 145 | + if (prv) |
| 146 | + { cur->nxt = prv->nxt; prv->nxt = cur; } |
| 147 | + else |
| 148 | + { cur->nxt = root; root = cur; } |
| 149 | + nxt = cur->nxt; if (nxt) nxt->prv = cur; |
| 150 | + cur->x = x; cur->a = df; cur->b = f-x*df; |
| 151 | + |
| 152 | + if (prv) prv->xmax = (cur->b-prv->b)/(prv->a-cur->a); |
| 153 | + if (nxt) cur->xmax = (nxt->b-cur->b)/(cur->a-nxt->a); |
| 154 | + } |
| 155 | + } |
| 156 | +} |
| 157 | + |
| 158 | + |
| 159 | +/* Private function to sample from piece-wise exponential distribution. First a |
| 160 | + * piece under the piece-wise distribution is sampled at random. Then a random |
| 161 | + * sample is drawn from this piece. A pointer to the segment, or piece which |
| 162 | + * was used is returned in q, and the function returns the random sample. Care |
| 163 | + * is taken to avoid numerical over and underflow. */ |
| 164 | + |
| 165 | +static double rpwed(struct segtype **q) |
| 166 | +{ |
| 167 | + double mass = 0.0, t, u; |
| 168 | + |
| 169 | + *q = root; while (*q) { mass += (*q)->mass; *q = (*q)->nxt; } |
| 170 | + t = mass*rand_uniform(); |
| 171 | + *q = root; while ((*q)->nxt && ((t -= (*q)->mass) >= 0.0)) *q = (*q)->nxt; |
| 172 | + |
| 173 | + u = rand_uniopen(); |
| 174 | + if ((*q)->prv == NULL) |
| 175 | + return (*q)->xmax+log(u)/(*q)->a; |
| 176 | + if ((*q)->nxt == NULL) |
| 177 | + return (*q)->prv->xmax+log(u)/(*q)->a; |
| 178 | + t = log(u+(1.0-u)*exp(fabs((*q)->a)*((*q)->prv->xmax-(*q)->xmax)))/(*q)->a; |
| 179 | + return ((*q)->a > 0) ? (*q)->xmax+t : (*q)->prv->xmax+t; |
| 180 | +} |
0 commit comments