李航《統計學習方法》第2版 第9章課後習題答案

習題9.1

習題9.2

習題9.3

2個分模型，5個參數μ0，σ0，μ1，σ1，（a0,a1）

代碼實作：

import numpy as np
import itertools
import matplotlib.pyplot as plt
import time


class MyEM():
	def __init__(self, y, args, tol=0.001, K=2, max_iter=50):
		self.y = y
		self.r = None #分模型k對觀測資料yi的響應度

		#寫成清單的形式是為了讓後面疊代時友善把更新的參數append到後面
		self.al = [np.array(args[0],dtype="float16").reshape(K,1)] # 分模型權重
		self.u = [np.array(args[1], dtype="float16").reshape(K, 1)]  # 分模型均值
		self.si = [np.array(args[2], dtype="float16").reshape(K, 1)]  # 分模型方差的平方

		self.score = 0  # 局部最優解評分,即似然函數值
		self.tol = tol  # 疊代中止門檻值
		self.K = K  # 高斯混合模型分量個數
		self.max_iter = max_iter  # 最大疊代次數

	def gaussian(self):
		"""計算高斯分布機率密度"""
		return 1/np.sqrt(2*np.pi*self.si[-1])*(np.exp(-(self.y-self.u[-1])**2/(2*self.si[-1])))

	def updata_r(self):
		"""更新分模型k對觀測資料yi的響應度"""
		r_jk = self.al[-1]*self.gaussian()
		self.r = r_jk/r_jk.sum(axis=0)

	def updata_u_al_si(self):
		"""更新u al si 每個分模型k的均值、權重、方差的平方"""
		u = self.u[-1]
		self.u.append(((self.r * self.y).sum(axis=1) / self.r.sum(axis=1)).reshape(self.K, 1))
		self.si.append((((self.r * (self.y - u) ** 2).sum(axis=1) / self.r.sum(axis=1))).reshape(self.K, 1))
		self.al.append((self.r.sum(axis=1) / self.y.size).reshape(self.K, 1))

	def judge_stop(self):
		"""中止條件判斷"""
		#就是判斷參數前後整體變化不大
		a = np.linalg.norm(self.u[-1] - self.u[-2])
		b = np.linalg.norm(self.si[-1] - self.si[-2])
		c = np.linalg.norm(self.al[-1] - self.al[-2])
		return True if np.sqrt(a ** 2 + b ** 2 + c ** 2) < self.tol else False

	def _score(self):
		"""計算該局部最優解的'評分'，也就是似然函數值"""
		self.score = (self.al[-1]*self.gaussian()).sum()

	def fit(self):
		"""疊代訓練獲得預估參數"""
		self.updata_r() # 更新r_jk 分模型k對觀測資料yi的響應度
		self.updata_u_al_si()  # 更新u al si 每個分模型k的均值、權重、方差的平方

		for i in range(self.max_iter):
			if not self.judge_stop():  # 若未達到中止條件，重複疊代
				self.updata_r()
				self.updata_u_al_si()
			else:  # 若達到中止條件，計算該局部最優解的“評分”，疊代中止
				self._score()
				break


def main():
	start = time.time()
	#觀測資料
	y = np.array([-67, -48, 6, 8, 14, 16, 23, 24, 28, 29, 41, 49, 56, 60, 75]).reshape(1, 15)

	y_mean = y.mean()//1
	y_std = (y.std()**2)//1


	#選取多個不同初值進行疊代，最後對得到的各個估計值進行比較，從中選擇最好的
	al = [[i, 1 - i] for i in np.linspace(0.1, 0.9, 9)]
	u = [[y_mean + i, y_mean + j] for i in range(-10, 10, 5) for j in range(-10, 10, 5)]
	si=[[y_std+i,y_std+j] for i in range(-1000,1000,500) for j in range(-1000,1000,500)]

	results=[]  #用來儲存不同初值選擇對應的最後的似然函數值，用來比較

	#産生這三個參數的所有可能組合i
	for i in itertools.product(al,u,si):
		clf = MyEM(y, i)
		clf.fit()
		results.append([clf.al[-1],clf.u[-1],clf.si[-1],clf.score]) # 收集不同初值收斂的局部最優解
		best_value = max(results, key=lambda x: x[3])  # 根據score即似然函數值，從所有局部最優解找到相對最優解,

	print("alpha:{}".format(best_value[0]))
	print("mean:{}".format(best_value[1]))
	print("std:{}".format(best_value[2]))


	# 可視化看看,似然函數值的分布（看樣子主要分布在0.15）
	re = [res[3] for res in results]
	plt.hist(re)
	plt.show()

	end = time.time()
	print("用時{:.2f}s".format(end-start))


if __name__ == '__main__':
    main()
           

習題9.4

參考這裡：習題9.4