Как оценить параметры смеси 2 экспоненциальных случайных величин (в идеале в Python)

Question

Представьте себе имитационный эксперимент, в котором на выходе получается 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, то и то, и другое достаточно для оценка б). В идеале я хочу сделать это на 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))

К сожалению, это не сработало. Для некоторых выборов параметров это в пределах порядка величины, но для других, это абсурдно. Учитывая мою проблему (с ее ограничениями) и мою цель оценить больший параметр двух экспонент (не заботясь о меньшем параметре или количестве точек, которые были получены от любого из них), есть идеи?

Answer 1

Я разместил вопрос в более общих статистических терминах на статистике 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) - th процентиль) для экспоненциального распределения (i-й квантиль = −ln(1- i) / λ, поэтому λ = −ln(1- i) / (i-й квантиль)). Результатом является список, в котором каждый i-й элемент соответствует оценке b с использованием (i+1) -го процентиля.

Затем я выполняю пиковые вызовы в этом списке с помощью реализации Python функции пиковых вызовов Matlab. Затем я беру список полученных пиков и возвращаю минимум. Кажется, это работает довольно хорошо.

Я также реализую решение EM в посте Stack Exchange и посмотрю, какой из них работает лучше.

РЕДАКТИРОВАТЬ: Я внедрил решение EM, и, кажется, он работает прилично хорошо в моих моделированиях (n = 1000, различные a и b).