1 处理数据 import numpy as np # 预处理数据 def loadData(filename): dataSet = [] fr = open(filename) for line in fr.readlines(): curLine = line.strip('\n').split('\t') fltLine = list(map(float, curLine)) dataSet.append(fltLine) fr.close() return dataSet # 计算高斯函数 def Gaussian(data,mean,cov): dim = np.shape(cov)[0] # 计算维度 covdet = np.linalg.det(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) else: covinv = np.linalg.inv(cov) # 计算cov的逆 #print(data,mean) 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) # 返回概率密度值

2 获取初始聚类中心 # 获取最初的聚类中心 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

3 GMM主程序 def GMM(data,K,ITER): N = data.shape[0] dim = data.shape[1] means= Kmeans(data,K) #means=GetInitialMeans(data,K,0.03) convs=[0]*K # 初始方差等于整体data的方差 for i in range(K): convs[i]=np.cov(data.T) #convs=np.full((K,dim,dim),np.diag(np.full(dim,0.1))) phi = [1.0/K] * K omega = [np.zeros(K) for i in range(N)] loglikelyhood = 0 oldloglikelyhood = 1 while np.abs(loglikelyhood - oldloglikelyhood) > 0.00001: #print(np.abs(loglikelyhood - oldloglikelyhood)) #while ITER: oldloglikelyhood = loglikelyhood # E步 for i in range(N): res = [phi[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个混合高斯的概率 omega[i][k] = res[k] / sumres # M步 for k in range(K): Nk = np.sum([omega[n][k] for n in range(N)]) # N[k] 表示N个样本中有多少属于第k个高斯 phi[k] = 1.0 * Nk/N means[k] = (1.0/Nk)*np.sum([omega[n][k] * data[n] for n in range(N)],axis=0) xdiffs = data - means[k] convs[k] = (1.0/ Nk)*np.sum([omega[n][k]* xdiffs[n].reshape(dim,1) * xdiffs[n] for n in range(N)],axis=0) # 计算最大似然函数 loglikelyhood = np.sum( [np.log(np.sum([phi[k] * Gaussian(data[n], means[k], convs[k]) for k in range(K)])) for n in range(N)]) ITER-=1 #print(oldloglikelyhood,loglikelyhood) return phi,means,convs

在GMM中用到的Kmeans算法如下: # K均值算法,估计大约几个样本属于一个GMM import copy def Kmeans(data,K): N = data.shape[0] # 样本数量 dim = data.shape[1] # 样本维度 means = GetInitialMeans(data,K,0.03) 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.0001: means_old = copy.deepcopy(means) numlog = [1] * 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 computeOmega(X,mu,sigma,phi,multiGaussian): n_samples=X.shape[0] n_clusters=len(phi) gamma=np.zeros((n_samples,n_clusters)) p=np.zeros(n_clusters) g=np.zeros(n_clusters) for i in range(n_samples): for j in range(n_clusters): p[j]=multiGaussian(X[i],mu[j],sigma[j]) g[j]=phi[j]*p[j] for k in range(n_clusters): gamma[i,k]=g[k]/np.sum(g) return gamma def predict(data,m,c,p): pred=computeOmega(np.array(data),m,c,p,Gaussian) cluster_results=np.argmax(pred,axis=1) return cluster_results

4 测试身高体重数据集 d=[] with open('d:/dataset1.txt','r') as f: for line in f.readlines(): d.append(line.strip('\n').split('\t')) d1=np.array(d) d2=[list(map(float,i))for i in d1[:,:-1]] data=np.array(d2) t=d1[:,-1] p2,m2,c2=GMM(data,2,50) pt2=predict(data,m2,c2,p2) pt2

结果如下: array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0], dtype=int64)

p2,m2,c2的值分别为: p2 [0.7822324209405439, 0.21776757905945612] m2 [array([170.5300851 , 59.69283016]), array([175.14653505, 73.79246884])] c2 [array([[63.72099898, 54.66242696], [54.66242696, 73.46041058]]), array([[ 21.00207568, 14.73278643], [ 14.73278643, 140.29202475]])]

误差和错误率: t=d1[:,-1] t[t=='f']=1 t[t=='m']=0 t=list(map(int,t)) t=np.array(t) c=0 for i in t==pt2: if i==False: c+=1 print('错误数为:',c) print('错误率为:',round(c/len(t),3))

结果为:

错误数为: 176

错误率为: 0.389

python身高体重程序代码_python EM算法4(身高体重数据集)相关推荐

  1. python入门小程序代码_Python入门小程序(二)

    学习了Python编程从入门到实践的第九章,本次的内容是类的创建.对象的实例化以及继承等内容. 1. 创建一个名为Restaurant的类,其方法__init__()设置两个属性:restaurant ...

  2. python抖音表白程序代码_python教程之抖音同款表白神器——Python让你的七夕更完美!-Go语言中文社区...

    又到了一年一度的七夕!大家都准备送什么给自己心仪的对象呢?今天给大家带来python教程之抖音同款表白神器--Python让你的七夕更完美! 先上效果: python表白小程序 代码: from tk ...

  3. python程序代码_python基础二

    Python基础-注释的引入 注释的分类: <1>单行注释:以#开头,#右边的所有文字当作说明,而不是真正要执行的程序,起辅助说明作用 多行注释用三个单引号 ''' 或者三个双引号 &qu ...

  4. python抢票软件代码_Python抢票程序优化,可以选择车次和座次

    通过程序自动化去刷新并点击抢票,就有了这个 Python 抢票程序.这个程序是 Python 模拟手工去操作浏览器的,所以会因为各种网络或者其他因素导致程序终止.本文主要讲解增加车次选择功能和座次选择 ...

  5. python求平方根的代码_Python求解平方根的方法

    本文实例讲述了Python求解平方根的方法.分享给大家供大家参考.具体如下: 主要通过SICP的内容改写而来.基于newton method求解平方根.代码如下: #!/usr/bin/python ...

  6. python的简单程序代码_小白学编程?从一个简单的程序开始学习Python编程

    笔者思虑再三还是决定选择图文(因为百家的视频发布画质真不怎么样[囧]). 笔者学习编程的时间也挺长的,因为业余,因为时间不多,各种原因,自学编程的路特别难走.然后笔者发现,自己能为小白贡献一些力量,然 ...

  7. python画正方形的代码_Python编程练习:使用 turtle 库完成正方形的绘制

    Vuforia点击屏幕自动对焦,过滤UGUI的按钮 //点击屏幕自对对焦 #if UNITY_EDITOR )) #elif UNITY_ANDROID || UNITY_IPHONE &&a ...

  8. python多线程执行同样代码_Python 多线程、多进程 (一)之 源码执行流程、GIL

    一.python程序的运行原理 许多时候,在执行一个python文件的时候,会发现在同一目录下会出现一个__pyc__文件夹(python3)或者.pyc后缀(python2)的文件 Python在执 ...

  9. python读文件完整代码_Python读写文件的代码示例

    本篇文章给大家带来的内容是关于Python读写文件的代码示例,有一定的参考价值,有需要的朋友可以参考一下,希望对你有所帮助 一.读取文件 读取文件步骤: 1.找到文件 2.打开文件 3.读取文件内容 ...

  10. python中空格的代码_python 空格

    初学python,不明白代码之间时空格的用处 比如: print "Hens", 25 + 30 / 6 print"hens",25+30/6 一个有空格一个 ...

最新文章

  1. 【CSS练习】IT修真院--练习4-移动端界面
  2. Spark详解(十一):Spark运行架构原理分析
  3. Scala中=gt;的用法
  4. 计算机学office有必要吗,计算机二级office要学多久
  5. [2018.09.08 T1] 炉石
  6. 《关于动态社交网络建模和分析的教程》的读书笔记
  7. cad脚本合适_CAD脚本学习
  8. 抽奖随机滚动_如何在party上用来宾的照片抽奖
  9. [转载]三、二、一 …… Geronimo!,第 3 部分: 状态问题
  10. 家用路由器与企业路由器有什么区别
  11. Extracting Multiple-Relations in One-Pass with Pre-Trained Transformers [论文研读]
  12. java 数组 eqlue_Java源码浅析,Character(3)
  13. 零点漂移、零点补偿问题
  14. 一款手机电脑都能用的进销存财务软件
  15. Png图片换色的方法
  16. java+如何画一个扇形_实现一个扇形的几种方法
  17. 河北计算机软件职业技术学院,河北软件职业技术学院2021年排名
  18. TI DM36X 名词
  19. robotframework中文乱码---robotframework日志输出时出现中文以unicode编码方式
  20. Unity URP 2020 设置DOTS

热门文章

  1. python爬取公众号之 创建个人微信公众号
  2. python彩色螺旋线_python绘制彩色螺旋线
  3. 鸡啄米c语言入门,鸡啄米编程课堂-最适合程序员在线学习和参考的教程站
  4. 吴式太极拳的特点-和基础要求
  5. C++编程-牛客网-雀魂启动
  6. 51单片机、STM32中生成QRCode二维码
  7. 服务器抓不到mrcp信息,MRCP学习笔记-语音识别资源的事件和headers详解
  8. 电驴服务器软性文件,电驴服务器.doc
  9. Oracle审计与数据库防火墙(AVDF)介绍
  10. 水溶性量子点CdSe/ZnS