首页 > 解决方案 > 使用 numpy 生成截断的对数正态分布

问题描述

我正在使用numpy.random.lognormal基于meanstd基础正态分布的分布来生成分布。

np.random.lognormal(mu, sigma, size=50)

我的问题是:我可以截断分布/样本以仅包含分布上下端的某些值吗?有没有办法指定分布的最小值和最大值?

标签: pythonnumpyrandom

解决方案


你的最后两句话对我来说是矛盾的。我在这里假设你想要一个区间内的值。如果您希望它们超出间隔,请给我写评论,我会更正我的代码。

为什么不只要求足够的对数正态分布值并获取区间内的结果呢?

import numpy as np
import math
import scipy
from scipy.stats import lognorm, binom
from itertools import count

def lognorm_in_interval(mu, sigma, k, loc=0, a=-np.inf, b=np.inf):
    s = sigma
    scale = math.exp(mu)
    dist = lognorm(s, loc, scale)
    p = dist.cdf(b)-dist.cdf(a)
    
    needed_tries = calc_needed_tries(p, k)
    x = dist.rvs(size=needed_tries)
    x = x[(a <= x) & (x <= b)]
    if len(x) >= k:
        return x[:k]
    else:
        np.array([*x, *lognorm_in_interval(mu, sigma, k-len(x), a, b)])
    
def calc_needed_tries(p, k):
    """calculates the amount of i.i.d. tries that are needed to have >= 95% 
    probability of an event with probability p to accur k times"""
    
    assert 0 < p, "this only works for events with positive probability"
    def prop(n): return 1-binom(n,p).cdf(k-1)-0.95
    m = next((m for m in count() if prop(10**m) >= 0))
    sol = scipy.optimize.root_scalar(prop, bracket=[10**(m-1),10**m])
    assert sol.converged
                                   
    return int(math.ceil(sol.root))
        
lognorm_in_interval(0,1,50,a=0.1,b=0.2)

推荐阅读