# Sampling from a recursively defined distribution

I'd like to know if there are techniques for sampling from a recursively defined probability distribution, assuming that solving the recursion for a formula for the distribution is too difficult.

As an example of what I have in mind by a recursively defined distribution take $p(t) = p(t-1)^2 + p(t-2)^3$, $t \in \mathbb{N}$, and choose the initial values that make it a probability distribution.

-
For the example you have, this sequence would converge to zero fairly quickly. Why not store the first probabilities until reaching the underflow bound, numerical precision-wise, and then do accept reject with some other distribution (say Negative binomial) whose parameters you would have optimized for the accept-reject algorithm to have good properties? – an12 Oct 22 '12 at 4:18

All you need to do is solve your recursion and then sample using the resulting probabilities.

So for example in R you might assume $p(0)=0$ and use the following

max_n         <- 20            # number of probabilities to calculate
p             <- rep(0, max_n)
rec_p         <- function(x1, steps=20){
p[1]   <- x1
p[2]   <- p[1]^2
for (i in 3:steps){ p[i] <- p[i-1]^2 + p[i-2]^3 }
p
}
err_rec_p     <- function(x){ abs(sum(rec_p(x, max_n)) - 1) }
best_p1       <- optimize( f=err_rec_p, interval=c(0,1), tol = 1e-9 )$minimum probabilities <- rec_p(best_p1, max_n) / sum(rec_p(best_p1, max_n))  to get the probabilities  [1] 5.005662e-01 2.505666e-01 1.882088e-01 5.115401e-02 9.283564e-03 [6] 2.200409e-04 8.485179e-07 1.137392e-11 6.110475e-19 1.471774e-33 [11] 2.281524e-55 3.188038e-99 1.187613e-164 3.240190e-296 0.000000e+00 [16] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00  and then sample using these probabilities samplesize <- 75 sample(1:max_n, samplesize, replace=TRUE, prob=probabilities )  which might for example give the following values  [1] 1 1 2 3 1 3 4 2 2 1 1 1 2 1 3 1 2 5 1 3 3 1 2 1 1 1 1 1 3 1 1 2 1 1 3 2 3 1 [39] 2 1 3 2 3 2 2 3 1 1 2 2 1 3 1 1 1 1 1 2 2 1 3 1 1 1 2 1 1 3 1 3 1 3 1 1 1  - Since the distribution you gave will go to zero very rapidly, an efficient sampling method might simply be the obvious one given as follows (assuming your initial probabilities$p(1)$and$p(2)$are known): 1. Set$k=1$. 2. Sample$u$uniformly from$[0,1]$. 3. If$u < p(k)$, choose$t_k$and exit. 4. Otherwise, set:$u \leftarrow u - p(k)p(k+2) = p(k+1)^2 + p(k)^3k \leftarrow k+1$And repeat from step 3. I hope it's clear what is being done here. We're simply using the CDF method and counting on the fact that on average$k\$ will be incremented very few times before accepting an answer.

-