【问题标题】:Is there an inverse probability mass function for Poisson statistics in Python?Python中的泊松统计是否有逆概率质量函数?
【发布时间】:2020-05-16 13:36:06
【问题描述】:

我正在寻找泊松统计的逆 pmf。 inverse 我的意思是一个函数inv_pmf(p, k),它返回分布参数 lambda。为清楚起见,参数使用如下:p = lambda^k / k! * exp(-lambda)。 谢谢

【问题讨论】:

    标签: python statistics distribution poisson


    【解决方案1】:

    所以你有概率方程 p(k,λ) = λk e/k!。你知道pk 但想知道λ。那么,从 lhs 和 rhs 中获取 log 并得到简单的方程。

    log(p) = k*log(λ) - λ - log(k!)

    λ = k*log(λ) - log(p) - log(G(k+1)),其中G() 是 Gamma 函数,在 Python lib 中可用。您可以绘制 RHS 和 LHS 之间的差异,并看到它可能有多种解决方案。他们使用 Python fsolve 函数,你可以从那个非线性方程中求根。

    代码(Python 3.7、Anaconda、Windows 10 x64)

    #%%    
    import math
    import numpy as np
    import matplotlib.pyplot as plt
    
    def logrhs(p, k, λ):
        return k*np.log(λ) - math.log(p) - math.lgamma(k+1)
    
    def poissonPMF(k, λ):
        lp = k*np.log(λ) - λ - math.lgamma(k+1)
        return np.exp(lp)
    
    p = 0.2
    k = 3.0
    
    λλ = np.linspace(0.001, 10.0, 101)
    
    #%%    
    rhs = logrhs(p, k, λλ)
    lhs = np.copy(λλ)
    
    pmf = poissonPMF(k, λλ)
    
    plt.plot(λλ, lhs - rhs, 'r')
    plt.plot(λλ, pmf, 'g')
    plt.show()
    
    # %%
    from scipy.optimize import fsolve
    
    def f(x):
        return x - logrhs(p, k, x)
    
    starting_guess = 4.0
    λ = fsolve(f, starting_guess)
    print((λ, poissonPMF(k, λ)))
    
    starting_guess = 1.9
    λ = fsolve(f, starting_guess)
    print((λ, poissonPMF(k, λ)))
    

    例如,我使用 k=3 和 p=0.2 进行测试。对我来说,它打印两个根

    (array([3.90263215]), array([0.2]))
    (array([2.24859448]), array([0.2]))
    

    并验证确实概率等于 0.2。图表显示红线在两个地方穿过 0

    【讨论】:

    • 谢谢,工作正常。但我假设在 scipy 或 numpy 中没有实现这样的功能?在那种情况下,我确实需要像你这样的功能。
    • @maximumphysics 我不知道。请注意,您以某种方式拥有所有概率 p_k,您可以取回 λ 作为平均值,λ = Sum(0...Infinity) k*p_k。对于像我这样的单值函数应该可以工作,但我演示的问题很简单——你可以得到两个答案,所以你必须使用一些额外的(特定领域的)知识来选择一个。
    • @maximumphysics 如果您认为答案可以接受,您可以按照 stackoverflow 规则接受它stackoverflow.com/help/someone-answers
    • @Severin_Pappadeux 谢谢,所以我将不得不使用像你这样的功能。
    猜你喜欢
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 2017-10-17
    • 1970-01-01
    • 1970-01-01
    相关资源
    最近更新 更多