如何用Python实现从c^a/(-a log(1-c))分布采样自然数a？

你的目标分布是对数级数分布（Logarithmic Series Distribution），其概率质量函数为：

P(a) = -cᵃ/(a·ln(1-c)) ，其中a取1,2,3,...，且0<c<1

下面提供两种Python实现方案，解决你无法用scipy.stats.rv_discrete处理无限取值的问题：

方案1：用Scipy内置的对数级数分布（最推荐） #

Scipy已经内置了这个分布的实现——scipy.stats.logser，无需手动处理无限取值的问题，直接调用采样方法即可：

import scipy.stats as stats

def sample_log_series(c, sample_size=1):
    # logser的参数p对应你的c，分布定义完全匹配
    return stats.logser.rvs(p=c, size=sample_size)

测试示例：

# 采样10个样本，c=0.3
samples = sample_log_series(0.3, sample_size=10)
print(samples)

这个方法简洁高效，Scipy内部已经优化了采样逻辑，是最优选择。

方案2：手动实现逆变换采样（适合理解底层逻辑） #

如果不想依赖Scipy，或者需要自定义采样逻辑，可以用逆变换采样。对数级数分布的CDF没有解析表达式，但可以通过迭代累加概率找到满足条件的a：

1. 生成[0,1)区间的均匀随机数U
2. 从a=1开始，逐步累加每个a的概率，直到累积概率超过U，此时的a就是采样结果

代码实现：

import numpy as np

def sample_log_series_manual(c, sample_size=1):
    samples = []
    log_term = np.log(1 - c)  # 预计算，减少重复计算
    for _ in range(sample_size):
        u = np.random.uniform(0, 1)
        cum_prob = 0.0
        a = 1
        while True:
            prob_a = - (c ** a) / (a * log_term)
            cum_prob += prob_a
            if cum_prob >= u:
                samples.append(a)
                break
            a += 1
    return np.array(samples)

注意：当c接近1时，分布尾部会更长，迭代次数会增加，但对于常规的c值（比如c<0.9），这个方法的效率足够满足需求。

补充说明：为什么rv_discrete不适用？ #

scipy.stats.rv_discrete要求预先提供所有可能的取值和对应的概率元组，但对数级数分布的支持是无限的自然数，无法预先枚举所有a，因此确实不适合用这个类来实现。

内容的提问来源于stack exchange，提问作者serene