# Sample from an infinite discrete distribution where probabilities are known up to a normalization constant

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$$?

• 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. – Ben Nov 7 '18 at 0:37
• @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. – Ben Nov 7 '18 at 0:51
• 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? – MFG Nov 7 '18 at 1:03
• 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?) – 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

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))