0

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

我的目标是估计a或b（不关心k，我不需要同时估计a和b：两者之一就可以了）。从推测来看，似乎只估计b可能是最简单的路径（当a << b时，几乎没有什么可用于估计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))

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

4

1 回答 1

0

我在 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个分位数是 ( i *100)-指数分布（第i个分位数 = −ln(1 − i ) / λ，因此 λ = −ln(1 − i ) / ( i 第 i个分位数））。结果是一个列表，其中每个第i个元素对应于使用第 ( i +1) 个百分位数的b估计值。

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

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

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

于 2017-07-15T07:57:07.570 回答