正規分布に従う標本から母集団の分散を最尤推定で計算してみました。最尤推定での計算はTensorFlowで勾配降下法を使いました。
最尤推定で計算した分散は不偏分散ではなく標本の分散と一致することが確認できました。
手順
- 正規分布に従う乱数を100個作成
- その平均と分散と不偏分散を計算
- 母集団が正規分布に従うと仮定したときの母集団の平均と分散をパラメータにとる尤度関数を設定
- 平均と分散を勾配降下法で変化させて、尤度関数が最大となる平均と分散を探す
- 尤度関数が最大となる平均と分散が2で計算した値と一致するかどうかを見る
Pythonコード
Google Colaboratoryで実行しました。
import numpy as np import scipy.stats as stats import matplotlib.pyplot as plt import math import tensorflow as tf # 標本を作成 n = 100 x = range(n) y = np.random.randn(n) # 標本の平均と分散を計算 desc = stats.describe(y) mean = desc.mean uvariance = desc.variance variance = uvariance * (n-1) / n # 尤度を計算する関数 logpi2 = 0.5 * math.log(2.0 * math.pi) def prob(): d = y - l_mean d2 = -0.5 * d * d / l_variance return tf.math.reduce_sum(d2) / n - 0.5 * tf.math.log(l_variance) - logpi2 # 目的となる損失関数 def loss(): return -prob() maxLoopCount = 30 # 適当な初期値 mean_ini = 0.0 variance_ini = 1.0 # 複数のOptimizerを試す opts = [ tf.optimizers.Adam(learning_rate=0.2), CustomOptimizer(learning_rate=0.1), ] for i in range(len(opts)): # 最初に正解を表示 print((mean, variance, uvariance)) # グラフにするための配列 lossHistory = [] l_mean = tf.Variable(mean_ini) l_variance = tf.Variable(variance_ini) for loopCount in range(maxLoopCount): l = float(loss()) # グラフにするために記録 lossHistory.append(l) # 最適化 opts[i].minimize(loss, var_list = [l_mean, l_variance]) # 学習の推移を数値で表示 if loopCount % 10 == 0: print((float(l_mean), float(l_variance))) # グラフ化 plt.rcParams['figure.figsize'] = (6.4, 4.8) plt.plot(range(maxLoopCount), lossHistory) plt.show()
CustomOptimizer
はTensorFlowのOptimizerを自作するを参照。
結果
実行結果例
OptimizerがCustomOptimizerの場合
コードの実行結果です。1行目は100個の乱数の平均・分散・不偏分散です。4行目は最尤推定の平均、分散です。最尤推定の分散は乱数の不偏分散ではなく通常の分散と一致していることがわかります。
(0.1080843015362518, 1.0207865528862736, 1.031097528167953) # 平均、分散、不偏分散 (0.010808431543409824, 1.0016233921051025) (0.10808192193508148, 1.0207873582839966) (0.1080842986702919, 1.0207866430282593) # 30回ループ後の最尤推定の平均、分散
横軸はループ回数、縦軸は損失関数の値
OptimizerがAdamの場合
Adamですと収束が不十分でよくわかりません。
(0.1080843015362518, 1.0207865528862736, 1.031097528167953) # 平均、分散、不偏分散 (0.19999417662620544, 1.1999610662460327) (0.15666580200195312, 1.0132553577423096) (0.12916100025177002, 1.0091229677200317) # 30回ループ後の最尤推定の平均、分散
横軸はループ回数、縦軸は損失関数の値
わかること
わかることは2つ
リンク
関連する私の記事