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
$\begingroup$
$\endgroup$
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
sim_size=10000000
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))
