統計とイモリとバーベルと

勉強の備忘録メイン

ブラックゴーストを飼いたい欲に駆られている.あとスネークヘッド.

メモ: scipy.optimize.minimizeによる最尤法の実装

Missing Data関連の準備として. MCARとかMARとかMNARとかの話をするときにどうしても最尤法は必要になってくるので今のうちに実装.

# -*- coding: utf-8 -*-

import numpy as np
import pandas as pd

from scipy.stats import norm, poisson, uniform
from scipy.optimize import minimize

missing_ratio = 0.3
sample_size = 10000
seed = 123

def generate_sample(sample_size, random_state):
    np.random.seed(random_state)
    df = pd.DataFrame({
        "y1": np.random.normal(2, 5, sample_size),
        "y2": np.random.randn(sample_size),
        "r": np.ones(sample_size, dtype="int")
    })
    df["r"].where(df.y2 >= norm.ppf(missing_ratio), other=0, inplace=True)
    df["y1_obs"] = df["y1"].where(df.r == 1, other=np.nan)
    return df

def _loglike_norm(x, *args):
    return -np.sum(np.log(norm.pdf(args[0], x[0], x[1]))) 

def mle_norm(data, x0=None):
    if x0 == None:
        x0 = np.random.randn(2)

    estimates_ = minimize(_loglike_norm,
                          x0,
                          args=(data),
                          method='Nelder-Mead',
                          options={'maxiter': 100,
                                   'disp': True})
    return estimates_

def main():
    df = generate_sample(sample_size, seed)
    print(mle_norm(data=df.y1, x0=[1, 10]))

if __name__ == "__main__":
    main()

初期値によってはうまく収束してくれないので注意が必要. このあたり,最初にグリッドサーチをしていい塩梅の値を見つけるみたいなやり方ができそう. 最適化アルゴリズムにはNelder-Mead法を選択しているのですが,このあたりは不勉強につき... アウトプットはこんな感じ.欠測のないy1はN(2, 5)を想定したので,ほぼ期待通りの数字が出てきています.

Optimization terminated successfully.
         Current function value: 30264.885192
         Iterations: 43
         Function evaluations: 85
 final_simplex: (array([[ 2.04856667,  4.99058462],
       [ 2.04856811,  4.99053511],
       [ 2.0486284 ,  4.9905318 ]]), array([ 30264.88519227,  30264.88519265,  30264.88519369]))
           fun: 30264.885192273156
       message: 'Optimization terminated successfully.'
          nfev: 85
           nit: 43
        status: 0
       success: True
             x: array([ 2.04856667,  4.99058462])

おまけ

optimizeの仕様をよく理解していなかったがために,作ってしまったけど使わなかった関数(泣). 確率密度関数を引数にとれたらいいなと思って組んだのですが,optimizeの引数となる関数はf(params, args)の形をとる必要があるので,これはちょっと難しい.argsに確率密度関数を含められるような作りにしてもいいけど,args = (pdf, data, ...) みたいな感じになってわかりやすい作りではなくなりますよね...
クラスを定義すればいけるけど,まぁ今回はいいかなと.

def loglike(pdf, data, **kwargs):
    return np.sum(np.log(pdf(data, **kwargs)))

    df = generate_sample(sample_size, seed)
    params = {"loc": 0,
              "scale": 1}
    print(loglike(data=df.y1, pdf=norm.pdf, **params))

    d = np.random.random_integers(0, 10, 3)
    params = {"loc": 0,
              "scale": 10}
    print(loglike(data=d, pdf=uniform.pdf, **params))