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]])