EnsekiTT Blog

EnsekiTTが書くブログです。

Python+NumpyでFIML(完全情報最尤推定法)試してたら色々勘違いをしていたっぽい話

スポンサーリンク

こんにちは、えんせきです。
タリーズブラジルコレクション ブラジル ファゼンダ プラナウトがうまくてハマりました。
いただきものの電動ミルがあるので、豆を買ってきて家で挽いて飲んでます。
*1

つまりなにしたの?

欠損値の補完方法を最近色々試していて、FIML(完全情報最尤推定法)という手法に行き着いた。
FIMLは欠損値を"補完"するわけじゃなくて、欠損値が含んでいても統計量を推定できるってことだったっぽい。

【追記: 2018/03/12】
正規化したら文献中MPlusで推定された結果と同じ値になりました。
この辺、今でもちょっと疑っているので詳しい文献があったら読もうと思います。また訂正いれるかもしれません。
f:id:ensekitt:20180313000033j:plain

前回はこちら
ensekitt.hatenablog.com


FIMLとは?

欠損データ分析(missing data analysis)-完全情報最尤推定法と多重代入法-
http://koumurayama.com/koujapanese/missing_data.pdf
欠損有りデータでも最尤推定で平均ベクトルと共分散行列が出せるってことだった。
*2

平均と共分散行列は作れるだろうけど結局何の値で欠損値を埋めたら良いんじゃろ?と思っていた。
いや、埋めなくても平均と共分散行列行列作れれば良いよね。見たいな話だった。

FIMLの実装

github.com
を写経しつつ、関数化されている部分を書き下してみた。

依存パッケージはこちら

import numpy as np
import scipy as sp
import scipy.optimize
  • 初期平均ベクトルと初期共分散行列を用意して最大化するために変更可能なパラメータベクトルにする
    size, dim = data.shape
    mean0 = np.zeros(dim)
    cov0 = np.eye(dim)

    params0 = np.empty(dim + dim * (dim + 1) // 2)
    params0[:dim] = mean0
    params0[dim:] = cov0[np.tril_indices(dim)]
  • 欠損のあるデータと無いデータのブロックに分ける
    obsmap = ~np.isnan(data)
    sortedidx = sorted(range(data.shape[0]), key=lambda i: list(obsmap[i]))
    blocks = [[sortedidx[0]]]
    for idx, prev in zip(sortedidx[1:], sortedidx[:-1]):
        if (obsmap[prev] == obsmap[idx]).all():
            blocks[-1].append(idx)
        else:
            blocks.append([idx])
    data_blocks = [(obsmap[b[0]], data[b][:, obsmap[b[0]]]) for b in blocks]
  • 対数尤度を平均ベクトルと共分散行列を更新しながら最大化(コード上は最小化)する
    result = sp.optimize.fmin_slsqp(
        _obj_func, params0, args=(dim, data_blocks), disp=True, iter=1000)

最小化したい関数

def _obj_func(params, dim, data_blocks):
    mean = params[0:dim]
    cov = np.empty((dim,dim))
    ii, jj = np.tril_indices(dim)
    cov[ii, jj] = params[dim:]
    cov[jj, ii] = params[dim:]
    
    if (np.linalg.eigvalsh(cov) < 0).any():
        # XXX Returning inf is not a good idea, because many solvers
        # cannot cope with it.
        return np.inf
    objval = 0.0
    for obs, obs_data in data_blocks:
        obs_mean = mean[obs]
        obs_cov = cov[obs][:, obs]
        objval += _log_likelihood_composed(obs_data, obs_mean, obs_cov)
    return -objval

対数尤度

def _log_likelihood_composed(x, mean, cov):
    xshift = x - mean
    size = x.shape[0]
    t1 = x.shape[-1] * np.log(2*np.pi)
    sign, logdet = np.linalg.slogdet(cov)
    t2 = logdet
    t3 = -0.5 * xshift.dot(np.linalg.inv(cov)) * xshift
    return (-0.5 * size * t1) + (-0.5 * size * t2) + t3.sum()

多変量正規分布を仮定した確率密度関数の対数尤度を最大化するので、
\log{L(u, \Sigma)} = -\frac{dN}{2}\log{2\pi} - \frac{N}{2}\log{|\Sigma|} + \sum_{i=0}^{N} \{-\frac{1}{2}(x_{i} - u)^T\Sigma^{-1}(x_{i} - u)\}
を最大化(コード上は最小化)した。
ただし、L(u, Σ)は尤度関数、|・|は行列式、dはデータの次元数、Nはデータ数、uは平均値、Σは共分散行列。

  • 平均ベクトルと共分散行列を返す
    mean = result[0:dim]
    cov = np.empty((dim, dim))
    ii, jj = np.tril_indices(dim)
    cov[ii, jj] = result[dim:]
    cov[jj, ii] = result[dim:]
    return mean, cov

データの例はこちら

data = np.array([
    [3, 83, np.nan],
    [4, 85, np.nan],
    [5, 95, np.nan],
    [2, 96, np.nan],
    [5, 103, 128],
    [3, 104, 102],
    [2, 109, 111],
    [6, 112, 113],
    [3, 115, 117],
    [3, 116, 133]], dtype=np.float)

http://koumurayama.com/koujapanese/missing_data.pdf
のTable1をそのまま作った。

【追記: 2018/03/12】
データセットの欠損値を含まない列を正規化する

mean_list = []
var_list = []
for i in range(data.shape[1]):
    if ~(~np.isnan(data[:,i]).any()):
        continue
    mean_list.append(data[:,i].mean())
    var_list.append(data[:,i].var())
    data[:,i] = (data[:,i] - mean_list[i]) / var_list[i]
mean, cov = fiml(data)
print(mean)
print(cov)
#===実行結果===#
mean:
array([-1.6090693e-06,  1.0780008e-06,  1.1091582e+02], dtype=float32)
cov:
array([[6.09783647e-01, 5.85975098e-04, 1.79081202e+00],
       [5.85975098e-04, 7.97747777e-03, 7.86245443e-01],
       [1.79081202e+00, 7.86245443e-01, 1.73048676e+02]])

mean[2]に110.9という値が出ていて、「欠損データ分析」のP7に書いてある「7Mplus で推定」した値と一致した。
最初、正規化しないでやってたら推定値が合わなくて泣きそうだった。
最尤推定で平均ベクトルと共分散行列を求められた。

つなげたコード

gist.github.com

クリエイティブ・コモンズ・ライセンス
この 作品 は クリエイティブ・コモンズ 表示 4.0 国際 ライセンスの下に提供されています。