如何估计 2 个指数随机变量的混合参数（最好在 Python 中）答案

【问题标题】：how to estimate parameters of mixture of 2 exponential random variables (ideally in Python)如何估计 2 个指数随机变量的混合参数（最好在 Python 中）
【发布时间】：2017-07-15 16:52:30
【问题描述】：

想象一个模拟实验，其中输出是 n 个总数，其中 k 个是从具有 a率的指数随机变量中采样的> 和 n-k 是从具有 b 率的指数随机变量中采样的。约束条件是 0 a ≤ b 和 0 ≤ k ≤ n，但是 a, b 和 k 都是未知的。此外，由于模拟实验的细节，当a ，k ≈ 0，当a = b， k ≈ n/2.

我的目标是估计 a 或 b （不关心 k，我不需要同时估计两者a 和 b：两者之一即可）。从推测来看，似乎只估计 b 可能是最简单的路径（当 a 时，几乎没有什么可用于估计 a 和大量可估计 b，当 a = b 时，两者仍有大量可估计 b)。理想情况下，我想用 Python 来做，但我对任何免费软件都持开放态度。

我的第一种方法是使用sklearn.optimize 来优化似然函数，其中，对于我的数据集中的每个数字，我计算 P(X=x) 的指数为 a，计算与比率 b 的指数相同，只需选择两者中的较大者：

from sys import stdin
from math import exp,log
from scipy.optimize import fmin
DATA = None

def pdf(x,l): # compute P(X=x) for an exponential rv X with rate l
    return l*exp(-1*l*x)

def logML(X,la,lb): # compute the log-ML of data points X given two exponentials with rates la and lb where la < lb
    ml = 0.0
    for x in X:
       ml += log(max(pdf(x,la),pdf(x,lb)))
    return ml

def f(x): # objective function to minimize
    assert DATA is not None, "DATA cannot be None"
    la,lb = x
    if la > lb: # force la <= lb
        return float('inf')
    elif la <= 0 or lb <= 0:
        return float('inf') # force la and lb > 0
    return -1*logML(DATA,la,lb)

if __name__ == "__main__":
    DATA = [float(x) for x in stdin.read().split()] # read input data
    Xbar = sum(DATA)/len(DATA) # compute mean
    x0 = [1/Xbar,1/Xbar] # start with la = lb = 1/mean
    result = fmin(f,x0,disp=DISP)
    print("ML Rates: la = %f and lb = %f" % tuple(result))

不幸的是，这并没有很好地工作。对于某些参数选择，它在一个数量级之内，但对于其他参数，它是荒谬的。鉴于我的问题（有其约束）和我估计两个指数的较大参数的目标（不关心较小的参数，也不关心来自任何一个的点数），有什么想法吗？

【问题讨论】：

标签： python scikit-learn statistics estimation exponential-distribution

【解决方案1】：

我在 stats Stack Exchange 上以更一般的统计术语发布了这个问题，它得到了答案：

https://stats.stackexchange.com/questions/291642/how-to-estimate-parameters-of-mixture-of-2-exponential-random-variables-ideally

另外，我尝试了以下方法，效果很好：

首先，对于每一个整数百分位数（第 1 个百分位数、第 2 个百分位数、...、第 99 个百分位数），我使用分位数封闭式方程计算 b 的估计值（其中 i-th 分位数是指数分布的 (i *100)-th 个百分位数 (i-th 分位数 = -ln(1 - i) / λ，所以 λ = −ln(1 - i) / (i-th quantile))。结果是一个列表，其中每个第 i 个元素对应于使用第 (i+1) 个百分位数的 b 估计。

然后，我使用 Matlab 峰值调用函数的 Python 实现对该列表执行峰值调用。然后，我获取结果峰值列表并返回最小值。它似乎工作得很好。

我还将在 Stack Exchange 帖子中实施 EM 解决方案，看看哪个效果更好。

编辑：我实现了 EM 解决方案，它似乎在我的模拟中运行良好（n = 1000，各种 a 和 b)。

【讨论】：