Take the 2-minute tour ×
MathOverflow is a question and answer site for professional mathematicians. It's 100% free, no registration required.

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.

share|improve this question
1  
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

2 Answers 2

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
share|improve this answer

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)^3$

    $k \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.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.