Let $X$ be a random variable taking values in the set of natural numbers. For each $i >0$ let $p_i = P[X = i]$. Say I don't know $p_1,p_2,...$ but I do know $q_i = ap_i$ for every $i > 0$, where $a$ is a unknown constant. If I were able to compute $a = \sum_i q_i$, then I known how to sample from the distribution of $X$. Otherwise, if $a$ is difficult to compute, is there a method to obtain samples from the distribution of $X$?

  • $\begingroup$ This is an interesting question. It is usually not computationally prohibitive to sum a series of numbers (particularly since they decreased to zero), so can you give any explanation of the context where this problem arises? If summing a series of numbers is computationally prohibitive, then presumably there will be other computational limitations on the kinds of methods that can be applied. $\endgroup$ – Ben Nov 7 '18 at 0:37
  • $\begingroup$ @MGF: Okay, that seems like a very different question to me (and one that has an analytic solution that is well-known). Anyway, thanks for elaborating. $\endgroup$ – Ben Nov 7 '18 at 0:51
  • $\begingroup$ Well, that is just the context in which I came up with the question, I was just wondering, in the general scenario (described in the question) is there a general method to sample from the desired distribution? $\endgroup$ – MFG Nov 7 '18 at 1:03
  • $\begingroup$ Do you know where the maximum of the $q_i$ values is? Alternatively, is it feasible to idenfify it? Do you know if it's unimodal (where for a discrete with an interval where there's a flat region with $q_i=q_{i+1}=...=q_{i+k}$ which are all the highest q-value, we still call that unimodal). Do you know anything else about it? (e.g. is it strictly monotonic-decreasing away from the mode? Is it decreasing at least as fast as exponential?) $\endgroup$ – Glen_b Nov 7 '18 at 5:16

An exact simulation is possible if you can find another sequence from which you can draw simulations $\{p_i\}_{i>0}$ such that $\max\left\{\frac{q_i}{p_i}\right\}=M<\infty$ by means of rejection sampling.

The algorithm goes as follows:

  • Obtain a sample $j$ from distribution $Q$ and a sample $u$ from $\mathrm{Unif}(0,M)$.
  • Check if $u<\frac{q_j}{p_j}$.
    • If this holds, accept $j$ as a sample drawn from $Q$;
    • if not, repeat from begining.

That being said, if $Q$ is very far from $P$ (in the Lévy-Prokhorov metric (e.g. if $Q$ is an exponential variable shifted to 1e10 and $P$ is a discretised Log-Cauchy) getting a sample with this method will be more computationally intensive that just crawling the sequence $q_i$ up to a certain value, cap it and normalize it.

Here's a sample code using Python

import numpy as np
from scipy.stats import poisson
import matplotlib.pyplot as plt

lam = 10
def seq(n, lam):
    return poisson.pmf((np.ceil(n/5)*5), lam)

def sup(n, lam):
    return poisson.pmf(n, lam)

Y = np.random.poisson(lam=lam, size=sim_size)
U = np.random.uniform(size=sim_size)
X = Y[U<seq(Y, lam)/sup(Y, lam)/84]

gp = np.unique(X, return_counts=True)

gp2 = np.unique(Y, return_counts=True)

plt.plot(gp[0], gp[1]/len(X))
plt.plot(gp2[0], gp2[1]/len(Y))

enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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