正規分布の母分散をTensorFlowを使って最尤推定で数値計算してみる

正規分布に従う標本から母集団の分散を最尤推定で計算してみました。最尤推定での計算はTensorFlowで勾配降下法を使いました。

最尤推定で計算した分散は不偏分散ではなく標本の分散と一致することが確認できました。

手順

  1. 正規分布に従う乱数を100個作成
  2. その平均と分散と不偏分散を計算
  3. 母集団が正規分布に従うと仮定したときの母集団の平均と分散をパラメータにとる尤度関数を設定
  4. 平均と分散を勾配降下法で変化させて、尤度関数が最大となる平均と分散を探す
  5. 尤度関数が最大となる平均と分散が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()

CustomOptimizerTensorFlowの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回ループ後の最尤推定の平均、分散

横軸はループ回数、縦軸は損失関数の値

f:id:suzuki-navi:20200926144719p:plain

OptimizerがAdamの場合

Adamですと収束が不十分でよくわかりません。

(0.1080843015362518, 1.0207865528862736, 1.031097528167953) # 平均、分散、不偏分散
(0.19999417662620544, 1.1999610662460327)
(0.15666580200195312, 1.0132553577423096)
(0.12916100025177002, 1.0091229677200317) # 30回ループ後の最尤推定の平均、分散

横軸はループ回数、縦軸は損失関数の値

f:id:suzuki-navi:20200926144734p:plain

わかること

わかることは2つ

リンク

関連する私の記事