2017-10-29 19 views
0

私はscipy.minimizeを使って簡単な最小化(シミュレーションされた最尤の基本的な例)を実行しようとしています。何らかの理由で初期値を返すだけです。私は間違って何をしていますか?ここでscipy minimの戻り値の初期値

は私のコードです:

import numpy as np 
from scipy.optimize import minimize 

# Simulated likelihood function 
# Arguments: 
# theta: vector representing probabilities 
# sims: vector representing uniform simulated data, e.g. [0.43, 0.11, 0.02, 0.97, 0.77] 
# dataCounts: vector representing counts of actual data, e.g. [4, 10, 7] 
def simLogLikelihood(theta, sims, dataCounts): 

     # Categorise sims using theta 
     simCounts = np.bincount(theta.cumsum().searchsorted(sims)) 

     # Calculate probabilities using simulated data 
     simProbs = simCounts/simCounts.sum() 
     # Calculate likelihood using simulated probabilities and actual data 
     logLikelihood = (dataCounts*np.log(simProbs)).sum() 

     return -logLikelihood 


# Set seed 
np.random.seed(121) 

# Generate 'true' data 
trueTheta = np.array([0.1, 0.4, 0.5]) 
dataCounts = np.bincount(np.random.choice([0, 1, 2], 1000, p=trueTheta)) 

# Generate simulated data (random draws from [0, 1)) 
sims = np.random.random(1000) 

# Choose theta to maximise likelihood 
thetaStart = np.array([0.33, 0.33, 0.34]) 
bnds = ((0, 1), (0, 1), (0, 1)) 
cons = ({'type': 'eq', 'fun': lambda x: x.sum() - 1.0}) 

result = minimize(simLogLikelihood, x0=thetaStart, args=(sims, dataCounts), method='SLSQP', bounds=bnds, constraints=cons) 

bndsで境界が確率は0と1の間である必要があるという事実を反映consでの制約は、確率の合計が1にしなければならないということである。。)

私はこのコードを実行した場合

resultが含まれています

 fun: 1094.7593617864004 
    jac: array([ 0., 0., 0.]) 
message: 'Optimization terminated successfully.' 
    nfev: 5 
    nit: 1 
    njev: 1 
    status: 0 
success: True 
     x: array([ 0.33, 0.33, 0.34]) 

をだから、ただ1回の繰り返しを行い、返します私が始めた確率のベクトル。しかし、より低い目的を有する確率の別のベクトルを見つけることは容易である。 [0.1,0.4,0.5]。何がうまくいかないのですか?

+1

scipyではなくnumpyへの誤った参照に気付いたsaschaさんに感謝します。これらは修正されました。 – jms202

答えて

0

あなたの最適化問題は、非平滑きれいに見える(おそらくnp.bincount()のが、私はそれに深くつもりはありません)実際にそこにほとんどのオプティマイザのために悪いことです。また、制約があるため、2つのオプティマイザ(SLSQP、COBYLA)のみが残されています。どちらも滑らかさを前提としています。以下のようなプリントの追加

:最後にsimLogLikelihood

print(theta, -logLikelihood) 

は、数値微分(グラデーションを提供しなかったとして)の間に、scipyのダウンロードは、いくつかの小さな摂動をしようとしているが、客観的であること、を示します全く変化しない(滑らかでない)!

[ 0.33 0.33 0.34] 1094.75936179 
[ 0.33 0.33 0.34] 1094.75936179 
[ 0.33000001 0.33  0.34  ] 1094.75936179 
[ 0.33  0.33000001 0.34  ] 1094.75936179 
[ 0.33  0.33  0.34000001] 1094.75936179 
    fun: 1094.7593617864004 
    jac: array([ 0., 0., 0.]) 
message: 'Optimization terminated successfully.' 
    nfev: 5 
    nit: 1 
    njev: 1 
    status: 0 
success: True 
     x: array([ 0.33, 0.33, 0.34]) 

num-diffは大きなステップを取るように調整することはできますが、私はあなたの問題がここにうってつけだとは思わないでしょう!

クイックデモ(推奨されません):

result = minimize(simLogLikelihood, x0=thetaStart, args=(sims, dataCounts), 
       method='SLSQP', bounds=bnds, constraints=cons, options={'eps': 1e-2}) 
                  # much bigger num-diff steps 

出力:あなたが見

[ 0. 0. 1.] inf 
[ 0.21587719 0.2695045 0.51461833] 1013.80776084 
[ 0.23010601 0.28726799 0.48262602] 1012.05516321 
[ 0.23627513 0.29496961 0.46875527] 1010.48916647 
[ 0.2386537 0.29793905 0.46340726] 1010.13774627 
[ 0.23957593 0.29909039 0.46133369] 1009.0850268 
[ 0.2397671 0.29932904 0.46090387] 1008.96044271 
[ 0.23981532 0.29938924 0.46079545] 1008.96044271 
[ 0.23983943 0.29941934 0.46074124] 1008.96044271 
[ 0.23985149 0.29943439 0.46071414] 1008.96044271 
[ 0.23985751 0.29944192 0.46070058] 1008.96044271 
    fun: 1008.960442706361 
    jac: array([ 947.81880269, -52.71300484, 0.  ]) 
message: 'Optimization terminated successfully.' 
    nfev: 44 
    nit: 6 
    njev: 5 
    status: 0 
success: True 
     x: array([ 0.23985751, 0.29944192, 0.46070058]) 

、いくつかのケースでは、あなたの関数によって返された非有限値があること。何か非常に悪い!

だから、私は非常に滑らかなものを公式化しようとすることをお勧めします。