지수 분포 및 최대 우도 추정 (MLE) 파이썬 코드 예제

여기서는 지수 분포에 관한 기본적인 문제와 그와 관련된 최대 우도 추정량의 문제를 파이썬으로 풀어보고 싶습니다.

문제 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층 미분도 계산해, 항상 마이너스가 되는 것도 확인하지 않으면 안됩니다. 그러나 이 계산은 인간이 하면 간단하지만, 파이썬은 중간에 식변형을 하고 계산을 편하게 하는 일을 하지 않기 때문에 잘 안 된다. (´·ω·`) 쇼본.

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