出力が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の両方を推定する必要はありません。2 つのうちの 1 つだけで十分です)。推測によると、bだけを推定するのが最も簡単な方法のようです ( a << bの場合、 a を推定するために使用するものはほとんどなく、bを推定するのに十分なものがあり、a = bの場合、両方ともまだ十分に推定することができます見積もりb )。理想的には Python で実行したいのですが、任意のフリー ソフトウェアを使用できます。
私の最初のアプローチはsklearn.optimize
、尤度関数を最適化するために使用することでした。ここで、データセットの各数値について、レートaの指数関数に対して P(X=x) を計算し、レートbの指数関数に対して同じ計算を行い、単純に大きい方を選択します。 2の:
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))
残念ながら、これはあまりうまくいきませんでした。パラメータのいくつかの選択では、それは一桁以内ですが、他のものでは、ばかげて外れています. 私の問題 (制約付き) と、2 つの指数関数の大きい方のパラメーターを推定するという私の目標 (小さい方のパラメーターやどちらからのポイントの数も気にせずに) を考えると、アイデアはありますか?