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$?
 A: 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))


