指数分布と最尤推定量(MLE)のpythonのコードの例

ここでは指数分布に関する基本的な問題と、それに関連した最尤推定量の問題をpythonで解いてみたいと思います。

問題1

N氏のホームページのアクセスの間隔のT(単位：時間)は、指数分布\\
f(t) = 2 e^{-2t}~~~~ (t>0)\\

に従うという。このときP(T\leq 3)を求めよ。

まずは諸々の準備。

import sympy as sym
from sympy.plotting import plot
sym.init_printing(use_unicode=True)
%matplotlib inline
oo = sym.oo # 無限大
(x,t) = sym.symbols('x t')

次に確率密度関数のグラフを描いてみます。

expr =2* sym.exp((-2*t))
# 得られた関数のグラフを0から3まで図示。
plot(expr, (t, 0, 3))
#念のため確率密度関数になっているかも確認。
sym.integrate(expr, (t, 0, oo))

最後に与えられた問題を解きます。

sym.integrate(expr, (t, 0, 3))

問題2

今度は最尤推定量の問題を解いてみましょう。

N氏のホームページのアクセスの間隔のT(単位：時間)は、指数分布\\
f(t) = x e^{-xt}~~~~ (t>0)\\

に従うとする。Tの独立な観測値の組を任意に入力し、その値の組に対応するxの最尤推定量を求めよ。

尤度関数を定義。変数がtではなくxに変わります。(普通はxではなくλを使います。）

#まずは観測値の組を入力。
A = list(map(float, input().split()))
#次に尤度関数を定義。
expr = 1
for t in A:
    expr = expr*(x * sym.exp((-x)*t))

print(expr)
# xを変数とする尤度関数のグラフを0から3まで図示。ここでは3という数字は特に意味なし。何でもよい。
plot(expr, (x, 0, 3))

次に対数尤度関数を定義。

L = sym.log(expr)
print(L)
# 得られた関数のグラフを図示。
plot(L, (x, 0, 3))

微分して0になる場所を計算します。

L_1 = L.diff(x,1)
plot(L_1, (x, 1, 3))
sym.solve(L_1, x)

本来なら以下の様に2階微分も計算して、常にマイナスとなることも確かめなければいけません。しかしこの計算は人間が行えば簡単ですが、Pythonは途中で式変形をして計算を楽にするということをしないため、上手くいきません。(´·ω·`)ｼｮﾎﾞｰﾝ。

L_2 = L_1.diff(x,1)
plot(L_2, (x, 1, 2))
print(L_2)