目录
3、采用贝叶斯估计方法，求男女生肺活量的分布参数（方差已知，注明自己选定的参数情况）；
4、基于身高和体重，采用最小错误率贝叶斯决策，画出类别判定的决策面。并判断某样本的身高体重分别为(165,50)时应该属于男生还是女生？为(175,55)时呢？
模式识别第一次作业
一、背景：
基于提供的数据进行作业的分析。并假定男生、女生的身高、体重、鞋码、50m成绩、肺活量都服从正态分布。
二、作业内容：
1、以肺活量为例，画出男女生肺活量的直方图并做对比；
2、采用最大似然估计方法，求男女生肺活量的分布参数；
3、采用贝叶斯估计方法，求男女生肺活量的分布参数（方差已知，注明自己选定的参数情况）；
4、基于身高和体重，采用最小错误率贝叶斯决策，画出类别判定的决策面。并判断某样本的身高体重分别为(165,50)时应该属于男生还是女生？为(175,55)时呢？
三、参数估计公式推导
1、最大似然估计方法估计分布参数
前提：男女生肺活量服从正态分布，且待估参数是一个确定性的量。此题使用最大似然估计单变量肺活量的均值和方差两个参数。
思路：任何一个子集Xk 样本是从总体中独立抽取的，独立同分布。根据样本信息独立处理，利用类概率密度形式来估计参数，使似然函数最大对应的参数是待估参数。
样本满足独立且同分布，所以类概率密度为：
则似然函数为：
对均值进行最大似然估计(对μ 求导并令导函数为0)：
得到估计均值：
对方差进行最大似然估计（对σ2 求导并令导函数为0）：
得到估计方差为：
2、用贝叶斯方法估计分布参数
前提：男女生的肺活量服从正态分布（均值、方差已知）即先验概率分布已知（μ0,σ02 ），类概率密度呈正态分布(μ,σ2) ，后验概率服从正态分布（均值未知、方差都已知）。样本信息独立同分布，均值为X 。
思路：引入估计风险使得贝叶斯估计中满足最小风险原则，使得求出的估计参数满足最小风险。通过定义损失函数（平方误差损失函数）来进行贝叶斯估计，通常最小方差贝叶斯估计是贝叶斯最优估计。
定义损失函数（平方误差损失函数）：
经计算贝叶斯统计中总平均风险：
定义样本下的条件风险：
通过求导可得最小风险的贝叶斯估计量为：
已知先验概率满足正态分布，均值、方差已知：
后验概率由贝叶斯公式（分母恒定设为C1 ）得：
其中，
解的上式得：
将后验概率带入∎ 式中的最小风险贝叶斯估计量求出均值估计量：
所以对均值估计：
3、决策面方程推导：
前提：男女生身高体重都服从正态分布，是两个特征向量。
思路：基于最小错误率贝叶斯一个分类器，得到判别函数下，令相邻判别函数相等从而求出决策面。
判别函数初始形式：
似然函数是正态分布：
对初始判别函数采取对数形式：
相邻判别函数相等得决策方程:
四、题目解答过程：
1、直方图绘制过程：
2、最大似然估计和贝叶斯法估计分布参数解答过程：
3、决策判断解答过程：
五、结果展示并分析：
1、直方图绘制：
结果为：
分析：男生女生的肺活量符合正态分布，男生均值在4300左右，女生在3300左右。由图可看出男生比女生更符合正态分布图形，因为男生样本数更多一点，这也符合我们的尝试，当样本数据越大时，分布图会越接近正态分布图。
在做此题时，对于数据的清理很重要，排除不完整或错误信息。直到第四题做完之后我才反应过来可以将这个性别错误信息通过决策判断出最有可能的性别。但是这个题中仅仅是将性别错误更改，对于没有填写肺活量的同学，并没有根据别的特征来估算他的肺活量，也是这个题完成不够好的一面。
2、最大似然估计与贝叶斯估计结果展示并分析：
结果图为：
从结果中可以看出最大似然估计和贝叶斯估计分布参数差不多，我们将结果画图展现并分析他们的特点。最大似然估计计算出方差太大，说明肺活量每个样本起伏很大，说明每个人差别还是有点大，这也激励着我们要勤于锻炼身体。
两个参数估计分布结果比较图为：
上面两个图是把男女生的肺活量样本真实数据直方图（柱状图）与最大似然估计分布图（点横曲线）和贝叶斯估计分布图（实曲线）作比较，图中三者分布图相差不大。为了更明显看出两个估计方法造成的误差，图中选中一部分进行放大来更好看出二者差别。为了比较估计出来的分布图与样本数据图进行比较，分别对估计出来的两个分布概率系数进行一定程度的放大。其中，最大似然估计就是利用样本数据得到估计，没有先验概率，所以最大似然估计得到的结果和原数据基本一致；而贝叶斯估计分布图将先验知识看的比较重要，设的方差比较小，使得均值估计量更受先验知识的影响，样本数据受的影响比较小。
根据男女生数量来讲，男生样本数量更多，使得估计的参数分布与实际分布十分相似，女生估计分布相对于实际就会有所偏差，这说明当样本数量越大时，估计的越会接近真实值。
3、判断决策结果以及特征点分布图为：
分析：上图为决策面以及特征点分布情况，并且将题目中两个点以及性别错误点标注在图中，并在图中指出结果。在这个决策判断的过程中，对于算法推导了几次，明白他的算法的推导，所以相应的书写代码也会很快写出来。所以，对于算法的熟练便于代码的书写。
六、代码
import xlrd
import matplotlib.pyplot as plt
import numpy as np
import random
import scipy.stats as st
import math
import sympy
from matplotlib.patches import ConnectionPatch
# 打开所需要的excel
workbook = xlrd.open_workbook(filename=r'G:\研究生一年级\模式识别\模式识别作业\模式识别第一次作业\作业数据.xls')
# 获取工作簿的sheet名称
names = workbook.sheet_names()
print('工作簿的sheet名称:', names)
# 通过sheet名称来获取表格数据
table = workbook.sheet_by_name(sheet_name='Sheet1')
# 获取第一行表头内容
row1 = table.row_values(0)
print('表头内容：', row1)
# 获取sheet中有效行数
row = table.nrows
# 获取表头所在列号
for i in range(len(row1)):
if row1[i] == '性别 男1女0':
oneindex1 = i
elif row1[i] == '肺活量':
oneindex2 = i
elif row1[i] == '身高(cm)':
oneindex3 = i
elif row1[i] == '体重(kg)':
oneindex4 = i
# 获取对应列中的数据
sex = table.col_values(colx=oneindex1, start_rowx=1, end_rowx=None)
lung = table.col_values(colx=oneindex2, start_rowx=1, end_rowx=None)
height = table.col_values(colx=oneindex3, start_rowx=1, end_rowx=None)
weight = table.col_values(colx=oneindex4, start_rowx=1, end_rowx=None)
# 女生的肺活量
girl = list()
# 男生的肺活量
boy = list()
list1 = list()
# 错误数据的数量
a = 0
error = list()
# 进行数据处理，将性别，肺活量不符合的去除
for j in range(row-1):
if sex[j] == 0: # 性别是女性
if type(lung[j]) is float: # 判断肺活量是否是数值类型
girl.append(lung[j])
else:
h = height[j]
bw = weight[j]
list1.append(sex[j])
list1.append(ah)
list1.append(bw)
a += 1
error.append(a) # 收集错误信息数量
elif sex[j] == 1: # 性别是男性
if type(lung[j]) is float: # 判断肺活量是否是数值类型
boy.append(lung[j])
else:
ah = height[j]
bw = weight[j]
list1.append(sex[j])
list1.append(ah)
list1.append(bw)
a += 1
error.append(a) # 收集错误信息数量
else:
if type(sex[j]) is float:
if type(lung[j]) is float:
print('该信息可以修改弥补，优第四题中程序判断')
print('第四题中判断决策该信息是男性的')
boy.append(lung[j])
a += 1
error.append(a) # 收集错误信息数量
print('肺活量信息错误，提取的性别身高体重为：', list1)
# 显示错误信息数量
print('错误信息数：', error)
# 设置字体样式
plt.rcParams['font.family'] = 'Fangsong'
# 设置字体大小
plt.rcParams['font.size'] = '10'
# 绘制直方图
plt.figure()
plt.hist(boy, bins=10, alpha = 0.4, color= 'blue', histtype = "bar")
plt.hist(girl, bins=10, color='red', alpha = 0.5, histtype = "bar")
plt.legend(["男生肺活量直方图", "女生肺活量直方图"])
plt.title("男女生肺活量直方图比较")
plt.xlabel("肺活量")
plt.ylabel("频数")
plt.grid()
plt.show()
'''第二题解答代码'''
# 男生,女生肺活量服从正态分布N(μ,σ^2)
# 最大似然估计:对于正态分布来说，实际上最后的均值的估计值就是等于样本的均值，方差的估计值等于样本二阶中心矩
# 男生肺活量的均值：
boy_u = np.mean(boy)
# 显示3位小数
boy_u_1 = float("%.3f"%boy_u)
# 女生肺活量的均值：
girl_u = np.mean(girl)
# 显示3位小数
girl_u_1 = float("%.3f"%girl_u)
# 计算方差:
# 男生方差
boy_num = 0
for i in range(len(boy)-1):
boy_num += (boy[i]-float(boy_u_1))*(boy[i]-float(boy_u_1))
boy_S = boy_num/len(boy)
boy_S = float("%.3f"%boy_S)
# 女生方差
girl_num = 0
for i in range(len(girl)-1):
girl_num += (girl[i]-float(girl_u_1))*(girl[i]-float(girl_u_1))
girl_S = girl_num/len(girl)
girl_S = float("%.3f"%girl_S)
print('\n')
print(' 第二题解答如下:')
print('---------------------最大似然估计结果展示-------------------------')
print('最大似然估计男生肺活量均值为:', boy_u_1, '最大似然估计男生肺活量方差为:', boy_S)
print('最大似然估计女生肺活量均值为:', girl_u_1, '最大似然估计女生肺活量方差为:', girl_S)
'''第三题解答代码'''
# 贝叶斯估计方法去估计参数，对于正态分布来说，实际上最后的均值的估计值就是等于样本的均值和我们的先验均值的一个加权求和
# 先验概率分布参数选为最大似然估计求取的值：男生p(u)~N（4500,500）女生p(u)~N(3500,500)
#设置均值u的先验概率的方差均为500,总体方差大约是600000
sigma =600000.0
# 下面采用贝叶斯估计方法求男女肺活量中后验概率分布的u=uN
# 由于p(u/X)~N(uN,sigma)
# 最小错误率贝叶斯估计所估计出的均值为
uNb = (sum(boy)*500000+sigma*4300)/(500000*len(boy)+sigma)
uNg = (sum(girl)*500000+sigma*3000)/(500000*len(girl)+sigma)
# 取三位小数
uNb = float("%.3f"%uNb)
uNg = float("%.3f"%uNg)
print('\n')
print(' 第三题解答如下:')
print('------------------贝叶斯估计后验概率均值结果展示---------------------------')
print('贝叶斯估计男生的肺活量均值为：{}\n(男生先验概率均值取{}，先验方差取{}，假设数据分布方差为{})'.format(uNb,4500,500,sigma))
print('贝叶斯估计女生的肺活量均值为：{}\n(女生先验概率均值取{}，先验方差取{}，假设数据分布方差为{})'.format(uNg,3500,500,sigma))
'''估计分布与真实估计之间的比较'''
# 设置字体样式
plt.rcParams['font.family'] = 'Fangsong'
# 设置字体大小
plt.rcParams['font.size'] = '10'
# 正常显示负号
plt.rcParams['axes.unicode_minus'] = False
#plt.figure()
fig, ax = plt.subplots(1, 1)
# 绘制直方图
plt.hist(boy, bins=10, alpha = 0.4, color= 'blue', histtype = "bar")
# plt.hist(girl, bins=10, color='red', alpha = 0.5, histtype = "bar")
# 绘制曲线图
# 定义正态分布函数
def normal_distribution(x, mu, sigma):
return np.exp( -1 * ( (x-mu) ** 2) / ( 2 * (sigma ** 2)) )*250000 / (math.sqrt( 2 * np.pi ) * sigma)
# 初始化mu, sigma
mub1, sigmab1 = boy_u_1, boy_S/1000
mub2, sigmab2 = uNb, 600
xb1 = np.linspace(mub1 - 6 * sigmab1, mub1 + 6 * sigmab1, 100)
yb1 = normal_distribution(xb1, mub1, sigmab1)
plt.plot(xb1, yb1, color= 'red', alpha = 0.4, linestyle="-.")
xb2 = np.linspace(mub2 - 6 * sigmab2, mub2 + 6 * sigmab2, 100)
yb2 = normal_distribution(xb2, mub2, sigmab2)
plt.plot(xb2, yb2, color= 'green', alpha = 0.4)
plt.legend(["-.男生肺活量最大似然估计", "-男生肺活量贝叶斯估计", "男生肺活量样本直方图"])
plt.title("男生肺活量估计分布与真实分布比较")
# 嵌入局部放大图的坐标系
axins = ax.inset_axes((0.7, 0.3, 0.4, 0.3))
axins.plot(xb1, yb1, color='r', alpha=0.8, linestyle="-.")
axins.plot(xb2, yb2, color='g', alpha=0.8)
# 设置放大区间
# X轴的显示范围
xlim0 = 4000
xlim1 = 4700
# Y轴的显示范围
ylim0 = 150
ylim1 = 170
# 调整子坐标系的显示范围
axins.set_xlim(xlim0, xlim1)
axins.set_ylim(ylim0, ylim1)
# 原图中画方框
tx0 = xlim0
tx1 = xlim1
ty0 = ylim0
ty1 = ylim1
sx = [tx0,tx1,tx1,tx0,tx0]
sy = [ty0,ty0,ty1,ty1,ty0]
ax.plot(sx,sy, color="gray", linestyle="--")
# 画两条线
xy = (xlim0,ylim0)
xy2 = (xlim0,ylim1)
con = ConnectionPatch(xyA=xy2, xyB=xy, coordsA="data", coordsB="data",
axesA=axins, axesB=ax, linestyle="--")
axins.add_artist(con)
xy = (xlim1,ylim0)
xy2 = (xlim1,ylim1)
con = ConnectionPatch(xyA=xy2, xyB=xy, coordsA="data", coordsB="data",
axesA=axins, axesB=ax, linestyle="--")
axins.add_artist(con)
plt.xlabel("肺活量")
plt.ylabel("频数")
plt.grid()
axins.grid()
plt.show()
#plt.figure()
fig, ax = plt.subplots(1, 1)
# 绘制直方图
plt.hist(girl, bins=10, alpha = 0.4, color= 'mediumpurple', histtype = "bar")
# 绘制曲线图
# 定义正态分布函数
def normal_distribution(x, mu, sigma):
return np.exp( -1 * ( (x-mu) ** 2) / ( 2 * (sigma ** 2)) )*40000 / (math.sqrt( 2 * np.pi ) * sigma)
# 初始化mu, sigma
mub1, sigmab1 = girl_u_1, girl_S/1000
mub2, sigmab2 = uNg, 600
xb1 = np.linspace(mub1 - 6 * sigmab1, mub1 + 6 * sigmab1, 100)
yb1 = normal_distribution(xb1, mub1, sigmab1)
plt.plot(xb1, yb1, color= 'red', alpha = 0.4, linestyle="-.")
xb2 = np.linspace(mub2 - 6 * sigmab2, mub2 + 6 * sigmab2, 100)
yb2 = normal_distribution(xb2, mub2, sigmab2)
plt.plot(xb2, yb2, color= 'green', alpha = 0.4)
plt.legend(["-.女生肺活量最大似然估计", "-女生肺活量贝叶斯估计", "女生肺活量样本直方图"])
plt.title("女生肺活量估计分布与真实分布比较")
axins = ax.inset_axes((0.7, 0.2, 0.4, 0.3))
axins.plot(xb1, yb1, color='r', alpha=0.8, linestyle="-.")
axins.plot(xb2, yb2, color='g', alpha=0.8)
# 设置放大区间
# X轴的显示范围
xlim0 = 3000
xlim1 = 3500
# Y轴的显示范围
ylim0 = 25
ylim1 = 30
# 调整子坐标系的显示范围
axins.set_xlim(xlim0, xlim1)
axins.set_ylim(ylim0, ylim1)
# 原图中画方框
tx0 = xlim0
tx1 = xlim1
ty0 = ylim0
ty1 = ylim1
sx = [tx0,tx1,tx1,tx0,tx0]
sy = [ty0,ty0,ty1,ty1,ty0]
ax.plot(sx,sy, color="gray", linestyle="--")
# 画两条线
xy = (xlim0,ylim0)
xy2 = (xlim0,ylim1)
con = ConnectionPatch(xyA=xy2,xyB=xy,coordsA="data",coordsB="data",
axesA=axins,axesB=ax,linestyle="--")
axins.add_artist(con)
xy = (xlim1,ylim0)
xy2 = (xlim1,ylim1)
con = ConnectionPatch(xyA=xy2,xyB=xy,coordsA="data",coordsB="data",
axesA=axins,axesB=ax,linestyle="--")
axins.add_artist(con)
plt.xlabel("肺活量")
plt.ylabel("频数")
plt.grid()
axins.grid()
plt.show()
'''第四题代码'''
# 获取表头所在列号
for i in range(len(row1)):
if row1[i] == '身高(cm)':
oneindex1 = i
elif row1[i] == '体重(kg)':
oneindex2 = i
elif row1[i] == '性别 男1女0':
oneindex3 = i
# 获取对应列中的数据
height = table.col_values(colx=oneindex1, start_rowx=1, end_rowx=None)
weight = table.col_values(colx=oneindex2, start_rowx=1, end_rowx=None)
sex = table.col_values(colx=oneindex3, start_rowx=1, end_rowx=None)
j = 0
# 创建男生女生身高体重列表进行对数据的存储
heightb = list()
heightg = list()
weightb = list()
weightg = list()
# 数据处理，将男女身高体重都满足有值的情况下分别放于相应的列表中
for i in range(len(sex)):
if sex[i] == 0: # 性别是女性
if type(height[i]) is float: # 判断身高是否是数值类型
if type(weight[i]) is float: # 判断体重是否是数值类型
weightg.append(weight[i]) # 将女生体重添加到weightg
heightg.append(height[i]) # 将女生身高数据添加到heightg
elif sex[i] == 1: # 性别是男性
if type(height[i]) is float: # 判断身高是否是数值类型
if type(weight[i]) is float: # 判断体重是否是数值类型
weightb.append(weight[i]) # 将男生体重添加到weightb
heightb.append(height[i]) # 将男生身高添加到heightb
else:
c3 = i
a3 = height[i]
b3 = weight[i]
# 计算均值
heightb_u = np.mean(heightb)
heightb_u = float("%.3f"%heightb_u)
weightb_u =np.mean(weightb)
weightb_u = float("%.3f"%weightb_u)
heightg_u = np.mean(heightg)
heightg_u = float("%.3f"%heightg_u)
weightg_u = np.mean(weightg)
weightg_u = float("%.3f"%weightg_u)
# 男女生身高体重矩阵为：
boyx = np.array([[heightb], [weightb]])
girlx = np.array([[heightg], [weightg]])
# 计算均值矩阵
boymean = np.array([[heightb_u], [weightb_u]])
girlmean = np.array([[heightg_u], [weightg_u]])
# 协方差矩阵计算
boycov = np.cov(heightb, weightb)
girlcov = np.cov(heightg, weightg)
BoyCov = sympy.Matrix([[float(boycov[[0], [0]]), float(boycov[[0], [1]])], \
[float(boycov[[1], [0]]), float(boycov[[1], [1]])]])
GirCov = sympy.Matrix([[float(girlcov[[0], [0]]), float(girlcov[[0], [1]])], \
[float(girlcov[[1], [0]]), float(girlcov[[1], [1]])]])
# 计算ln（|Σboy| / |Σgirl|）
ln_Value = math.log(np.linalg.det(boycov) /np.linalg.det(girlcov), math.e)
# 计算求ln( P(w1) / P(w2) )
ln_Pw = math.log((len(heightb) / len(sex)) / (len(heightg) / len(sex)), math.e)
# 求协方差矩阵的逆
boycov_N = np.linalg.inv(boycov)
girlcov_N = np.linalg.inv(girlcov)
# 判别函数求解
h, w = sympy.symbols('h w')
Func = (0.5 * (sympy.Matrix([[h], [w]]) - boymean).T * (BoyCov ** (-1)) * (
sympy.Matrix([[h], [w]]) - boymean)).det() \
- (0.5 * (sympy.Matrix([[h], [w]]) - girlmean).T * (GirCov ** (-1)) * (
sympy.Matrix([[h], [w]]) - girlmean)).det() \
+ 0.5 * ln_Value - ln_Pw
SolveFunc = sympy.solve(Func, [h, w])
# 决策面显示
plt.figure()
Y_Weight = np.linspace(40, 70, 200)
X_Hight = 0.25 * Y_Weight + 80 * np.sqrt(-0.00093 * Y_Weight ** 2 + 0.084 * Y_Weight - 1) + 83
plt.plot(Y_Weight, X_Hight)
# 设置字体样式
plt.rcParams['font.family'] = 'Fangsong'
# 设置字体大小
plt.rcParams['font.size'] = '10'
# 正常显示负号
plt.rcParams['axes.unicode_minus'] = False
# 男女特征点显示
plt.scatter(weightb, heightb, s=16, color='green', alpha=0.7, marker='*')
plt.scatter(weightg, heightg, s=16,color='pink', alpha=0.8, marker='+')
plt.xlabel("身高")
plt.ylabel("体重")
plt.legend(["分界面", "男生身高体重分布点", "女生身高体重分布点"], loc="center right")
plt.title('身高体重分布')
plt.plot(55, 175, 'o')
plt.plot(50, 165, 'o')
plt.plot(60, 169, 'o')
plt.annotate(text='我身高175体重55在这里，是男性哦', xy=(55, 175), xytext=(125, 190),weight='bold', color='purple',\
arrowprops=dict(arrowstyle='-|>',connectionstyle='arc3',color='purple', alpha=0.8),\
bbox=dict(boxstyle='round,pad=0.5', fc='black', ec='k',lw=0.5 ,alpha=0.4))
plt.annotate(text='我身高165体重50在这里，是女性哦', xy=(50, 165), xytext=(25, 140),weight='bold', color='purple',\
arrowprops=dict(arrowstyle='-|>',connectionstyle='arc3',color='purple', alpha=0.8),\
bbox=dict(boxstyle='round,pad=0.5', fc='black', ec='k',lw=0.5 ,alpha=0.4))
plt.annotate(text='我身高169体重60在这里，是男性哦', xy=(60, 169), xytext=(125, 140),weight='bold', color='purple',\
arrowprops=dict(arrowstyle='-|>',connectionstyle='arc3',color='purple', alpha=0.8),\
bbox=dict(boxstyle='round,pad=0.5', fc='black', ec='k',lw=0.5 ,alpha=0.4))
plt.show()
# 判别固定的点来判断类型是男是女
x1 = np.array([[175], [55]])
x2 = np.array([[165], [50]])
x3 = np.array([[a3], [b3]])
Func1 = (0.5 * (x1 - boymean).T * (BoyCov ** (-1)) * (
x1 - boymean)).det() \
- (0.5 * (x1 - girlmean).T * (GirCov ** (-1)) * (
x1 - girlmean)).det() \
+ 0.5 * ln_Value - ln_Pw
Func2 = (0.5 * (x2 - boymean).T * (BoyCov ** (-1)) * (
x2 - boymean)).det() \
- (0.5 * (x2 - girlmean).T * (GirCov ** (-1)) * (
x2 - girlmean)).det() \
+ 0.5 * ln_Value - ln_Pw
Func3 = (0.5 * (x3 - boymean).T * (BoyCov ** (-1)) * (
x3 - boymean)).det() \
- (0.5 * (x3 - girlmean).T * (GirCov ** (-1)) * (
x3 - girlmean)).det() \
+ 0.5 * ln_Value - ln_Pw
print('\n')
print(' 第四题解答如下:')
print('性别出错人信息是第', c3, '列')
print('这个人的身高体重为:', a3, b3)
print('----判断结果类型如下:------')
if Func1 > 0:
print(' （175，55）的类型是女性')
else:
print(' （175，55）的类型是男性')
if Func2 > 0:
print(' （165，50）的类型是女性')
else:
print(' （165，50）的类型是男性')
if Func3 > 0:
print('性别错误身高体重（169，60）的类型是女性，供第一题使用')
else:
print('性别错误身高体重（169，60）的类型是男性，供第一题使用')
扫一扫
9480
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
