1硬币问题
先看一个抛硬币问题,如果我们有A和B两个不均匀硬币,选择任意一个硬币抛10次(这里我们知道选择是的哪一个硬币),共计选择5次。正面记为H,背面记为T。记录实验结果,求A和B再抛正面向上的概率?
使用极大似然估计(Maximum likelihood)来算:
- 统计出每次实验,正反面的次数
- 多次实验结果相加
- 相除得到结果,P(A)=0.8,P(B)=0.45
但是在实际过程中,很有可能我们只知道有两个硬币,不知道每次选择的哪一个硬币,问是否能求出每个硬币抛出正面的概率?
是不是第一感觉这个问题很无解,我都不知道选择是哪个硬币,我怎么求解啊?
事实上是可以求出的。这里使用的时EM算法。
使用期望最大值(Expectation maximization,EM)算法来算:
- 假设
$\hat{\theta}_A^{(0)}=0.6,\hat{\theta}_B^{(0)}=0.5$
- 统计每次的实验结果,记录正反面
- 通过贝叶斯公式,估计每次实验选择的A硬币或是B硬币的概率
- 依据计算出的选择硬币概率得到该概率下的正反面结果
- 5次平均得到
$\hat{\theta}_A^{(1)}\approx 0.71,\hat{\theta}_B^{(1)}\approx 0.58$
- 重复上面的过程,例如迭代10次后,得到
$\hat{\theta}_A^{(2)}\approx 0.8,\hat{\theta}_B^{(2)}\approx 0.52$
- $ `就是使用EM算法计算出的概率值
这里比较难理解的是:如何利用贝叶斯公式计算每次实验选择的A硬币或是B硬币的概率?
那么先看下面这个例子:
假如现在有射击运动员甲和乙,甲乙射击靶心的概率为$\hat{\theta}_{\text{甲}}=0.9$
,$
`,如果现在有一组实验结果为
**直观上的来看,非常大的概率是甲射击的,但是也有可能是乙走狗屎运了。那么该如何从概率的角度计算是谁射击的? **
首先我们知道选择甲和乙的概率为p(甲)=p(乙)=0.5(先验概率),本次实验记为E。通过贝叶斯公式
$ `
故本次实验有98%可能是A射击的,2%的可能是B射击的。
有了上面的案例,我们再回到抛硬币的问题上,由贝叶斯公式:
p(A\vert E)=\frac{p(E|A)P(A)}{p(E)}
A 为选用硬币A,E为本次实验。而选择两个硬币的概率是相同的:$p(A)=p(B)=0.5$
,且$p(E|A)=0.6^5\times 0.4^5\approx 0.0008,P(E|B)=0.5^5\times 0.5^5\approx 0.001$
\begin{array}{lll}
p(A|E)&=&\frac{p(E|A)p(A)}{p(E)}\\
&=&\frac{p(E|A)p(A)}{p(E|A)p(A)+p(E|B)p(B)}\\
&\approx & 0.45
\end{array}
有$0.45\times 10\times 0.5\approx 2.2$
,故第一次实验结果平均下来,有2.22.2个A硬币正面的可能。
最后将5次结果平均得到新的A,B抛硬币正面估计值$\hat{\theta}_A^{(1)}\approx 0.71,\hat{\theta}_B^{(1)}\approx 0.58$
,,这是我们第一次迭代的值(这就是一次学习过程),照着这个流程迭代多次,得到最后的估测值。
上面我们计算出每次实验中是抛A或抛B的概率值就是隐变量.这个过程就是EM算法的简单案例。
2形式化EM算法
2.1 算法理论
` $
` $
$ `
这是一个非凸问题,只能求局部极值,一般采用EM算法。
步骤如下:
(1).随机化` $
(2).E步:计算
上式表示:第i个样本落在第j个高斯的概率.
(3). M步:
(4).回到(2)直至收敛.
2.2 python实现
以西瓜数据集4.0为例:
编号 | 密度 | 含糖率 |
---|---|---|
1 | 0.697 | 0.46 |
2 | 0.774 | 0.376 |
3 | 0.634 | 0.264 |
4 | 0.608 | 0.318 |
5 | 0.556 | 0.215 |
6 | 0.403 | 0.237 |
7 | 0.481 | 0.149 |
8 | 0.437 | 0.211 |
9 | 0.666 | 0.091 |
10 | 0.243 | 0.267 |
11 | 0.245 | 0.057 |
12 | 0.343 | 0.099 |
13 | 0.639 | 0.161 |
14 | 0.657 | 0.198 |
15 | 0.36 | 0.37 |
16 | 0.593 | 0.042 |
17 | 0.719 | 0.103 |
18 | 0.359 | 0.188 |
19 | 0.339 | 0.241 |
20 | 0.282 | 0.257 |
21 | 0.748 | 0.232 |
22 | 0.714 | 0.346 |
23 | 0.483 | 0.312 |
24 | 0.478 | 0.437 |
25 | 0.525 | 0.369 |
26 | 0.751 | 0.489 |
27 | 0.532 | 0.472 |
28 | 0.473 | 0.376 |
29 | 0.725 | 0.445 |
30 | 0.446 | 0.459 |
import numpy as np
# 预处理数据
def loadData(filename):
dataSet = []
fr = open(filename)
for line in fr.readlines():
curLine = line.strip().split(' ')
fltLine = list(map(float, curLine))
dataSet.append(fltLine)
return dataSet
# 计算高斯函数
def Gaussian(data,mean,cov):
dim = np.shape(cov)[0] # 计算维度
covdet = np.linalg.det(cov) # 计算|cov|
covinv = np.linalg.inv(cov) # 计算cov的逆
if covdet==0: # 以防行列式为0
covdet = np.linalg.det(cov+np.eye(dim)*0.01)
covinv = np.linalg.inv(cov+np.eye(dim)*0.01)
m = data - mean
z = -0.5 * np.dot(np.dot(m, covinv),m) # 计算exp()里的值
return 1.0/(np.power(np.power(2*np.pi,dim)*abs(covdet),0.5))*np.exp(z) # 返回概率密度值
# 获取最初的聚类中心
def GetInitialMeans(data,K,criterion):
dim = data.shape[1] # 数据的维度
means = [[] for k in range(K)] # 存储均值
minmax=[]
for i in range(dim):
minmax.append(np.array([min(data[:,i]),max(data[:,i])])) # 存储每一维的最大最小值
minmax=np.array(minmax)
while True:
for i in range(K):
means[i]=[]
for j in range(dim):
means[i].append(np.random.random()*(minmax[j][1]-minmax[j][0])+minmax[j][0] ) #随机产生means
means[i]=np.array(means[i])
if isdistance(means,criterion):
break
return means
# 用于判断初始聚类簇中的means是否距离离得比较近
def isdistance(means,criterion=0.03):
K=len(means)
for i in range(K):
for j in range(i+1,K):
if criterion>np.linalg.norm(means[i]-means[j]):
return False
return True
dataSet=loadData('d:/watermelon4.txt')
means=GetInitialMeans(np.array(dataSet),3,0.03)
# K均值算法,估计大约几个样本属于一个GMM
def Kmeans(data,K):
N = data.shape[0] # 样本数量
dim = data.shape[1] # 样本维度
means = GetInitialMeans(data,K,15)
means_old = [np.zeros(dim) for k in range(K)]
# 收敛条件
while np.sum([np.linalg.norm(means_old[k] - means[k]) for k in range(K)]) > 0.01:
means_old = cp.deepcopy(means)
numlog = [0] * K # 存储属于某类的个数
sumlog = [np.zeros(dim) for k in range(K)]
# E步
for i in range(N):
dislog = [np.linalg.norm(data[i]-means[k]) for k in range(K)]
tok = dislog.index(np.min(dislog))
numlog[tok]+=1 # 属于该类的样本数量加1
sumlog[tok]+=data[i] # 存储属于该类的样本取值
# M步
for k in range(K):
means[k]=1.0 / numlog[k] * sumlog[k]
return means
def GMM(data,K):
N = data.shape[0]
dim = data.shape[1]
means= Kmeans(data,K)
convs=[0]*K
# 初始方差等于整体data的方差
for i in range(K):
convs[i]=np.cov(data.T)
pis = [1.0/K] * K
gammas = [np.zeros(K) for i in range(N)]
loglikelyhood = 0
oldloglikelyhood = 1
while np.abs(loglikelyhood - oldloglikelyhood) > 0.0001:
oldloglikelyhood = loglikelyhood
# E步
for i in range(N):
res = [pis[k] * Gaussian(data[i],means[k],convs[k]) for k in range(K)]
sumres = np.sum(res)
for k in range(K): # gamma表示第n个样本属于第k个混合高斯的概率
gammas[i][k] = res[k] / sumres
# M步
for k in range(K):
Nk = np.sum([gammas[n][k] for n in range(N)]) # N[k] 表示N个样本中有多少属于第k个高斯
pis[k] = 1.0 * Nk/N
means[k] = (1.0/Nk)*np.sum([gammas[n][k] * data[n] for n in range(N)],axis=0)
xdiffs = data - means[k]
convs[k] = (1.0/ Nk)*np.sum([gammas[n][k]* xdiffs[n].reshape(dim,1) * xdiffs[n] for n in range(N)],axis=0)
# 计算最大似然函数
loglikelyhood = np.sum(
[np.log(np.sum([pis[k] * Gaussian(data[n], means[k], convs[k]) for k in range(K)])) for n in range(N)])
print(means)
print(loglikelyhood)
GMM(np.array(dataSet),3)