こんにちは、えんせきです。
タリーズブラジルコレクション ブラジル ファゼンダ プラナウトがうまくてハマりました。
いただきものの電動ミルがあるので、豆を買ってきて家で挽いて飲んでます。
*1
つまりなにしたの?
欠損値の補完方法を最近色々試していて、FIML(完全情報最尤推定法)という手法に行き着いた。
FIMLは欠損値を"補完"するわけじゃなくて、欠損値が含んでいても統計量を推定できるってことだったっぽい。
【追記: 2018/03/12】
正規化したら文献中MPlusで推定された結果と同じ値になりました。
この辺、今でもちょっと疑っているので詳しい文献があったら読もうと思います。また訂正いれるかもしれません。
前回はこちら
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()
多変量正規分布を仮定した確率密度関数の対数尤度を最大化するので、
を最大化(コード上は最小化)した。
ただし、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 で推定」した値と一致した。
最初、正規化しないでやってたら推定値が合わなくて泣きそうだった。
最尤推定で平均ベクトルと共分散行列を求められた。
つなげたコード
*1:タリーズブラジルコレクション ブラジル ファゼンダ プラナウト|バライエタル(ストレート)|コーヒー豆|商品情報|TULLY'S COFFEE - タリーズコーヒー
*2:Table1の平均値間違っている気がするけど