前提・実現したいこと
現在、Python3.x系にて最尤推定を行おうとしている者です。
最小化問題を解くツールとしてはscipyを使用しております。
一応コードは回るのですがどうも結果の数字が大きすぎる、コードに問題があるのではないかと思い質問させていただきます。
目的関数(尤度関数)は以下URLの第3式になっています。
リンク内容
第4式、第5式は第3式の中で再帰的に計算される部分なのでこちらも目的関数の中身になります。
"ts2"という名前の配列は関数外で定義されているもので、リンク先の式では添え字付きのtで表されています。
条件式は①全てのパラメータ(theta=[alpha,beta,mu])が正、②beta/alphaが1未満、といったものになります。
該当のソースコード
Python
1#objective function 2def func(theta): #theta={alpha, beta, mu} 3 n=len(ts2) 4 R = np.zeros(n) 5 secondterm = 0 6 thirdterm = 0 7 for i in range(n-1): 8 if i != 0 : 9 R[i]= np.exp(-theta[1]*(ts2[i]-ts2[i-1]))*(1+R[i-1]) 10 secondterm += np.exp(-theta[1]*(ts2[n-1]-ts2[i]))-1 11 thirdterm += math.log(theta[2]+theta[0]*R[i]) 12 loglikelihood = -ts2[n-1]*theta[2]+(theta[0]/theta[1])*secondterm+thirdterm 13 14 return -loglikelihood 15 16#constraints 17def cons (theta): 18 return 1.0001-theta[0]/theta[1] 19 20cons = ( 21 {'type': 'ineq', 'fun': cons}, 22 {'type': 'ineq', 'fun': lambda theta:theta[0]-0.0001}, 23 {'type': 'ineq', 'fun': lambda theta:theta[1]-0.0001}, 24 {'type': 'ineq', 'fun': lambda theta:theta[2]-0.0001} 25) 26 27#initial values of theta 28theta = [0.03, 0.05, 0.06] 29 30# 31result = minimize(fun=func, x0=theta, constraints=cons, method='SLSQP', callback=cbf) 32print(result)
補足情報(FW/ツールのバージョンなど)
ここにより詳細な情報を記載してください。
あなたの回答
tips
プレビュー