admin管理员组

文章数量:1122861

ML

---

视频讲解地址

=copy_web

假如有一个事件，其有四种可能，发生的概率分别为： ( 1 2 − θ 4 ) , ( 1 4 − θ 4 ) , ( 1 4 + θ 4 ) , ( θ 4 ) (\frac{1}{2} - \frac{\theta}{4}),(\frac{1}{4}-\frac{\theta}{4}),(\frac{1}{4}+\frac{\theta}{4}),(\frac{\theta}{4}) (21​−4θ​),(41​−4θ​),(41​+4θ​),(4θ​)

概率	1 2 − θ 4 \frac{1}{2} - \frac{\theta}{4} 21​−4θ​	1 4 − θ 4 \frac{1}{4}-\frac{\theta}{4} 41​−4θ​	1 4 + θ 4 \frac{1}{4}+\frac{\theta}{4} 41​+4θ​	θ 4 \frac{\theta}{4} 4θ​
观测到的次数	y 1 y_1 y1​	y 2 y_2 y2​	y 3 y_3 y3​	y 4 y_4 y4​

试估计 θ \theta θ

极大似然估计

L ( θ ) = ( 1 2 − θ 4 ) y 1 ( 1 4 − θ 4 ) y 2 ( 1 4 + θ 4 ) y 3 ( θ 4 ) y 4 L(\theta) = (\frac{1}{2}-\frac{\theta}{4})^{y_1}(\frac{1}{4}-\frac{\theta}{4})^{y_2} (\frac{1}{4}+\frac{\theta}{4})^{y_3}(\frac{\theta}{4})^{y_4} L(θ)=(21​−4θ​)y1​(41​−4θ​)y2​(41​+4θ​)y3​(4θ​)y4​

l o g L ( θ ) = y 1 l o g ( 2 − θ 4 ) + y 2 l o g ( 1 − θ 4 ) + y 3 l o g ( 1 + θ 4 ) + y 4 l o g ( θ 4 ) logL(\theta) = y_1log(\frac{2-\theta}{4})+y_2log(\frac{1-\theta}{4}) + y_3log(\frac{1+\theta}{4}) + y_4log(\frac{\theta}{4}) logL(θ)=y1​log(42−θ​)+y2​log(41−θ​)+y3​log(41+θ​)+y4​log(4θ​)

对 l o g L ( θ ) logL(\theta) logL(θ)求导并令其为0，可得极值点。

欲估计 θ \theta θ，需求 l o g L ( θ ) logL(\theta) logL(θ)的极值点处对应的 θ \theta θ。遗憾的是，这里即使靠机器也无法算出，哈哈。

在以下代码中， t t t代表 θ \theta θ

# %%
import numpy as np
import pandas as pdimport sympy as sp
import matplotlib.pyplot as plt
import seaborn as sns# %%
# 含参数的概率模型
# 一个事件的四种可能的概率t = sp.symbols("t")
probs = [1.0 / 2 - t / 4.0, 1.0 / 4 - t / 4.0, 1.0 / 4 + t / 4.0, t / 4.0]# 观测到每个事件发生的次数
y1, y2, y3, y4 = sp.symbols("y1 y2 y3 y4")
nums = [y1, y2, y3, y4]# %%
# 传统的最大似然估计MLE
import mathLt = (probs[0] ** nums[0]) * (probs[1] ** nums[1]) * (probs[2] ** nums[2]) * (probs[3] ** nums[3])logLt = sp.log(Lt, t)logLt_diff = sp.diff(logLt, t)# sp.solve([logLt_diff],[t]) 很遗憾,解不出来

优化后的极大似然估计

为什么解不出来，因为 θ \theta θ最后包含了至少3次项，不容易搞啊。那么，如何优化这个问题呢？将各个事件分解，使得其概率分布尽量与其他事件有相同点。

将事件重新划分。

概率	1 4 − θ 4 \frac{1}{4} - \frac{\theta}{4} 41​−4θ​	1 4 \frac{1}{4} 41​	1 4 − θ 4 \frac{1}{4}-\frac{\theta}{4} 41​−4θ​	θ 4 \frac{\theta}{4} 4θ​	1 4 \frac{1}{4} 41​	θ 4 \frac{\theta}{4} 4θ​
观测到次数	z 1 z_1 z1​	y 1 − z 1 y_1-z_1 y1​−z1​	y 2 y_2 y2​	z 1 z_1 z1​	y 3 − z 1 y_3-z_1 y3​−z1​	y 4 y_4 y4​

此时，有：

L ( θ ) = ( 1 4 − θ 4 ) z 1 + y 2 ( 1 4 ) y 1 − z 1 ( θ 4 ) z 2 + y 4 ( 1 4 ) y 3 − z 2 L(\theta)=\left(\frac{1}{4}-\frac{\theta}{4}\right)^{z_{1}+y_{2}}\left(\frac{1}{4}\right)^{y_{1}-z_{1}}\left(\frac{\theta}{4}\right)^{z_{2}+y_4}\left(\frac{1}{4}\right)^{y_{3}-z_{2}} L(θ)=(41​−4θ​)z1​+y2​(41​)y1​−z1​(4θ​)z2​+y4​(41​)y3​−z2​

ln ⁡ L ( θ ) = ( z 1 + y 2 ) ln ⁡ 1 − θ 4 + ( z 2 + y 4 ) ln ⁡ θ 4 + ( y 1 − z 1 + y 3 − z 2 ) ln ⁡ 1 4 \ln L(\theta)=\left(z_{1}+y_{2}\right) \ln \frac{1-\theta}{4}+\left(z_{2}+y_{4}\right) \ln \frac{\theta}{4}+\left(y_{1}-z_{1}+y_{3}-z_{2}\right) \ln \frac{1}{4} lnL(θ)=(z1​+y2​)ln41−θ​+(z2​+y4​)ln4θ​+(y1​−z1​+y3​−z2​)ln41​

d ln ⁡ L ( θ ) d θ = − z 1 + y 2 1 − θ + z 2 + y 4 θ = 0 \frac{d \ln L(\theta)}{d \theta}=-\frac{z_{1}+y_{2}}{1-\theta}+\frac{z_{2}+y_{4}}{\theta}=0 dθdlnL(θ)​=−1−θz1​+y2​​+θz2​+y4​​=0

θ ^ = z 2 + y 4 z 1 + z 2 + y 2 + y 4 \hat{\theta}=\frac{z_{2}+y_{4}}{z_1 + z_{2}+y_{2}+y_{4}} θ^=z1​+z2​+y2​+y4​z2​+y4​​

没错，很容易的就解出来了，但是真的解出来了吗？ 毕竟我们不知道 z 1 ， z 2 z_1，z_2 z1​，z2​啊。所以这个问题被我们引入了另外的参数 z 1 z_1 z1​和 z 2 z_2 z2​。 这个时候，就需要我们的EM算法了。

# 优化后的MLE
t = sp.symbols("t")
probs = [1.0 / 4 - t / 4.0, 1.0 / 4, 1.0 / 4 - t / 4.0, t / 4.0, 1.0 / 4, t / 4.0]# 观测到每个事件发生的次数
y1, y2, y3, y4 = sp.symbols("y1 y2 y3 y4")
z1, z2 = sp.symbols("z1 z2")
nums = [z1, y1 - z1, y2, z2, y3 - z2, y4]Lt = ((1.0 / 4 - t / 4.0) ** (z1 + y1)) * ((1.0 / 4) ** (y1 - z1)) * ((t / 4.0) ** (z2 + y4)) * ((1.0 / 4) ** (y3 - z2))logLt = sp.log(Lt)logLt_diff = sp.diff(logLt, t)result = sp.solve([logLt_diff], [t])

EM

EM算法适用于多参数估计、数据缺失（假设为参数）的情况。就如上例中引入的新参数。

E步 消去参数

如上例中， Z 1 ∼ B ( y 1 , 1 − θ 2 − θ ) . Z 2 ∼ B ( y 3 , θ 1 + θ ) Z_{1} \sim B\left(y_1 , \frac{1-\theta}{2-\theta}\right) . \quad Z_{2} \sim B\left(y_{3}, \frac{\theta}{1+\theta}\right) Z1​∼B(y1​,2−θ1−θ​).Z2​∼B(y3​,1+θθ​)，那么就用它们的期望分别替代 z 1 z_1 z1​和 z 2 z_2 z2​。

E ( z 1 ) = 1 − θ 2 − θ y 1 E ( z 2 ) = θ 1 + θ y 3 E\left(z_{1}\right)=\frac{1-\theta}{2-\theta} y_{1} \quad E\left(z_{2}\right)=\frac{\theta}{1+\theta} y_{3} E(z1​)=2−θ1−θ​y1​E(z2​)=1+θθ​y3​

则代入之后可以得到：

θ ^ = θ 1 + θ y 3 + y 4 1 − θ 2 − θ y 1 + θ 1 + θ y 3 + y 2 + y 4 \hat{\theta}=\frac{\frac{\theta}{1+\theta} y_{3}+y_{4}}{\frac{1-\theta}{2-\theta} y_{1}+\frac{\theta}{1+\theta} y_{3}+y_{2}+y_{4}} θ^=2−θ1−θ​y1​+1+θθ​y3​+y2​+y4​1+θθ​y3​+y4​​

M步 不断迭代，求得最值

任意取 θ \theta θ，不断迭代，直到收敛。（已经有人证明一定会收敛）

θ ( i + 1 ) = θ ( i ) 1 + θ ( i ) y 3 + y 4 1 − θ ( i ) 2 − θ ( i ) y 1 + θ ( i ) 1 + θ ( i ) y 3 + y 2 + y 4 \theta^{(i+1)}=\frac{\frac{\theta^{(i)}}{1+\theta^{(i)}} y_{3}+y_{4}}{\frac{1-\theta^{(i)}}{2-\theta^{(i)}} y_{1}+\frac{\theta^{(i)}}{1+\theta^{(i)}} y_{3}+y_{2}+y_{4}} θ(i+1)=2−θ(i)1−θ(i)​y1​+1+θ(i)θ(i)​y3​+y2​+y4​1+θ(i)θ(i)​y3​+y4​​

# EM 算法
threshold = 1e-6
iter_func = result[0][0]t_tmp = 0.3t_old = 0.3
t_new = 0
idx = 0
max_iter = 1000
while (idx < max_iter):t_old = t_newz1 = ((1.0 - t) / (2.0 - t)) * y1z2 = (t / (1.0 + t)) * y3z1_val = z1.evalf(subs={t: t_old, y1: 100.0, y2: 100.0, y3: 100.0, y4: 100.0})z2_val = z2.evalf(subs={t: t_old, y1: 100.0, y2: 100.0, y3: 100.0, y4: 100.0})z1_temp = float(z1_val)z2_temp = float(z2_val)t_new = iter_func.evalf(subs={t: t_old, y1: 100.0, y2: 100.0, y3: 100.0, y4: 100.0})t_new = t_new.evalf(subs={'z1': z1_temp, 'z2': z2_temp})idx += 1print(f"t || {t_old}")

最后迭代到收敛的 θ \theta θ值，就是我们对 θ \theta θ的估计值了。

本文标签： ml