次元数が大きい多変量正規分布の計算

Question

以下のような式で多変量正規分布の計算をしています
また，例として，1500次元の変数，平均，分散を与えます

import numpy as np
x = np.random.rand(1500)
mean = np.random.rand(1500)
cov = np.random.rand(1500,1500)

def gaussian(x,mean,cov):
    A = 1. / ((2. * np.pi) ** (x.size/2.0))
    B = 1. / (np.linalg.det(cov) ** 0.5)
    C = -0.5 * np.dot(np.dot(x - mean, np.linalg.inv(cov)), x - mean)
    return A * B * np.exp(C)

gaussian(x,mean,cov)

このときにxの数が大きすぎて，Aの計算の段階で，オーバーフローしてしまいます
また，Bの計算ではnp.linalg.det(cov)が-infとなってしまいます

このような場合，応急措置として，それっぽい値を代入しているのですが
適切な方法があれば教えていただけるとありがたいです

Answer 1

logをとるんじゃないでしょうか
log は単調増加関数
L(θ1)<L(θ2) であれば logL(θ1)<logL(θ2)が成り立つ

Comments

[kohekoh]

式までのせていただきありがとうございます

Answer 2

オーバーフロー対策について
対数(log)変換が王道かと思います。これで乗算が加算に置き換わるのでかなり桁数増が抑制できるようになるはずです。ただし、変換は多変量正規分布の確率密度関数全体に対して実施しないとつじつまが合わないので注意願います。

値Bの-infについて
計算式というよりcovに問題があるように思います。ご存知の通りcovは共分散行列なので値が完全なランダムということはありえません。対角成分を軸にして表裏の位置にある成分は同値であるはずです。

なお1500次元ほど大きい多変量正規分布だと平均ベクトルでも得られる期待値はほぼゼロになりますが、この点は問題にならないのでしょうか。

Comments

[kohekoh]

なるほど
たしかに共分散行列なので対称行列だから
ということですね
ある変数が正規分布にしたがっていると仮定して、分散、平均を学習しています
更新のコードが間違っていて、
対称になっていない可能性があるかもです

[kohekoh]

共分散行列の更新式の間違いを見つけました
ありがとうございました

[kohekoh]

「なお1500次元ほど大きい多変量正規分布だと平均ベクトルでも得られる期待値はほぼゼロになりますが、この点は問題にならないのでしょうか。」

期待値がほぼ0の問題というのはどういうのがあるのでしょうか
そこについて理解ができていませんでした

[R.Shigemori]

期待値が常にほぼゼロしか出力されない確率密度関数だと、推計に使えないように思いました。ここでいうほぼゼロとは10の-10乗分の1といったスケール感ですが、これだと「ある変量の組み合わせがサンプル10000件に出現する件数」は常にゼロになります。

[kohekoh]

推計に使えないというのはどういうことでしょうか
やっても値が正しいとはいえないという感じですか？

[R.Shigemori]

値は正しいのですが、現実と整合しないという意味です。1500次元の確率密度関数が常に10の-10乗程度の期待値を出力しても値そのものは正当です。また、これを使ってある1500次元のデータがサンプル数10000件に含まれる件数が0件になることも正当です。
しかしながら、現実のデータをサンプリングすると推計上0件とされたデータが10000種類観測され、理論値と整合しない状況になります。

[kohekoh]

なるほど、そういうことですね
例えばなのですが、以下のような場合、文書数分の観測データを推計することになると思うのですが
この場合も同じように期待値が非常に小さくなってしまうようなきがするのですが
どうでしょうか
https://papers.nips.cc/paper/3328-supervised-topic-models.pdf

[R.Shigemori]

紹介いただいた論文は初見なので詳しいことはわかりませんが、Abstractを読む限り、トピックモデルの核となっている技術であるLDAを拡張したものかと思います。そうだとすると、ご質問の「文書数分の観測データを推計することになる」というものが何を意味しているのか私には分かりませんでした。
おそらく先の論文を十分に理解していないためのような気がするので、詳しく解説していただかないとコメントしようがないような気がしています。

[kohekoh]

この論文だと，単語の生成モデルと共に，教師データ（例えば，評価など）の生成モデルも含んでいます．
その教師データが正規分布を仮定していて，その平均，分散を学習していくというふうになっています．
論文中だと，5000や4000の文書を使っていて，その文書分の評価値があるので，教師データ（自分の質問でいうところのx）の数も5000ほどになって，それこそ期待値が非常に小さくなると思うのですが，という疑問です

[kohekoh]

勘違いしていたかもしれません
5000次元のxに対して正規分布を考えているのではなくて
スカラーとして，xを一つ一つ正規分布の確率密度を計算しているみたいです