EM算法(讲解+程序)

EM算法实验报告

一、算法简单介绍

EM 算法是Dempster，Laind，Rubin于1977年提出的求参数极大似然估计的一种方法，它可以从非完整数据集中对参数进行MLE估计，是一种非常简单实用的学习算法。这种方法可以广泛地应用于处理缺损数据、截尾数据以及带有噪声等所谓的不完全数据，可以具体来说，我们可以利用EM算法来填充样本中的缺失数据、发现隐藏变量的值、估计HMM中的参数、估计有限混合分布中的参数以及可以进行无监督聚类等等。

本文主要是着重介绍EM算法在混合密度分布中的应用，如何利用EM算法解决混合密度中参数的估计。

二、算法涉及的理论

我们假设X是观测的数据，并且是由某些高斯分布所生成的，X是包含的信息不完整（不清楚每个数据属于哪个高斯分布）。

，

此时，我们用k维二元随机变量Z（隐藏变量）来表示每一个高斯分布，将Z引入后，最终得到：，

，

然而Z的后验概率满足（利用条件概率计算）：

但是，Z nk为隐藏变量，实际问题中我们是不知道的，所以就用Z nk的期望值去估计它（利用全概率计算）。

然而我们最终是计算max：

最后，我们可以得到（利用最大似然估计可以计算）：

三、算法的具体描述

3.1 参数初始化

对需要估计的参数进行初始赋值，包括均值、方差、混合系数以及。

3.2 E-Step计算

利用上面公式计算后验概率，即期望。

3.3 M-step计算

重新估计参数，包括均值、方差、混合系数并且估计此参数下的期望值。

3.4 收敛性判断

将新的与旧的值进行比较，并与设置的阈值进行对比，判断迭代是否结束，若不符合条件，则返回到3.2，重新进行下面步骤，直到最后收敛才结束。

四、算法的流程图

五、实验结果

a_best=

0.8022 0.1978

mu_best=

2.7148

3.9307

4.9882 3.0102 cov_best=

(:,:,1) =

5.4082 -0.0693

-0.0693 0.2184

(:,:,2) =

0.0858 -0.0177

-0.0177 0.0769

f=

-1.6323

数据X的分布

每次迭代期望值

-50510

-50510利用EM估计的参量值与真实值比较（红色：真实值青绿色：估计值）

六、参考文献

1.M. Jordan. PatternRecognitionAndMachineLearning

2.Xiao Han. EM Algorithm

七、附录

close all;clear;clc;

% 参考书籍Pattern.Recognition.and.Machine.Learning.pdf

% https://www.wendangku.net/doc/b94653691.html,

% lwm@https://www.wendangku.net/doc/b94653691.html,

% 2009/10/15

%%

M=2; % number of Gaussian

N=200; % total number of data samples

th=0.000001; % convergent threshold

K=2; % demention of output signal

% 待生成数据的参数

a_real =[4/5;1/5];

mu_real=[3 4;

5 3];

cov_real(:,:,1)=[5 0;

0 0.2];

cov_real(:,:,2)=[0.1 0;

0 0.1];

% generate the data

x=[ mvnrnd( mu_real(:,1) , cov_real(:,:,1) , round(N*a_real(1)) )' , mvnrnd(mu_real(:,2),cov_real(:,:,2),N-round(N*a_real(1)))'];

% for i=1:round(N*a_real(1))

% while (~((x(1,i)>0)&&(x(2,i)>0)&&(x(1,i)<10)&&(x(2,i)<10)))

% x(:,i)=mvnrnd(mu_real(:,1),cov_real(:,:,1),1)';

% end

% end

%

% for i=round(N*a_real(1))+1:N

% while (~((x(1,i)>0)&&(x(2,i)>0)&&(x(1,i)<10)&&(x(2,i)<10)))

% x(:,i)=mvnrnd(mu_real(:,1),cov_real(:,:,1),1)';

% end

% end

figure(1),plot(x(1,:),x(2,:),'.')

%这里生成的数据全部符合标准

%% %%%%%%%%%%%%%%%% 参数初始化

a=[1/3,2/3];

mu=[1 2;2 1];%均值初始化完毕

cov(:,:,1)=[1 0;

0 1];

cov(:,:,2)=[1 0;

0 1];%协方差初始化

%% EM Algorothm

% loop

count=0;

figure(2),hold on

while 1

a_old = a;

mu_old = mu;

cov_old= cov;

rznk_p=zeros(M,N);

for cm=1:M

mu_cm=mu(:,cm);

cov_cm=cov(:,:,cm);

for cn=1:N

p_cm=exp(-0.5*(x(:,cn)-mu_cm)'/cov_cm*(x(:,cn)-mu_cm));

rznk_p(cm,cn)=p_cm;

end

rznk_p(cm,:)=rznk_p(cm,:)/sqrt(det(cov_cm));

end

rznk_p=rznk_p*(2*pi)^(-K/2);

%E step

%开始求rznk

rznk=zeros(M,N);%r(Z

pikn=zeros(1,M);%r(Z

pikn_sum=0;

for cn=1:N

for cm=1:M

pikn(1,cm)=a(cm)*rznk_p(cm,cn);

% pikn_sum=pikn_sum+pikn(1,cm);

end

for cm=1:M

rznk(cm,cn)=pikn(1,cm)/sum(pikn);

end

end

%求rank结束

% M step

nk=zeros(1,M);

for cm=1:M

for cn=1:N

nk(1,cm)=nk(1,cm)+rznk(cm,cn);

end

end

a=nk/N;

rznk_sum_mu=zeros(M,1);

% 求均值MU

for cm=1:M

rznk_sum_mu=0;%开始的时候就是错在这里，这里要置零。

for cn=1:N

rznk_sum_mu=rznk_sum_mu+rznk(cm,cn)*x(:,cn);

end

mu(:,cm)=rznk_sum_mu/nk(cm);

end

% 求协方差COV

for cm=1:M

rznk_sum_cov=zeros(K,M);

for cn=1:N

rznk_sum_cov=rznk_sum_cov+rznk(cm,cn)*(x(:,cn)-mu(:,cm))*(x(:,cn)-mu(:,cm))';

end

cov(:,:,cm)=rznk_sum_cov/nk(cm);

end

t=max([norm(a_old(:)-a(:))/norm(a_old(:));norm(mu_old(:)-mu(:))/norm(mu_old(:));norm(cov_ol d(:)-cov(:))/norm(cov_old(:))]);

temp_f=sum(log(sum(pikn)));

plot(count,temp_f,'r+')

count=count+1;

if t

break;

end

end %while 1

hold off

f=sum(log(sum(pikn)));

a_best=a;

mu_best=mu;

cov_best=cov;

f_best=f;

% 输出结果

disp('a_best=');

disp(a_best);

disp('mu_best=');

disp(mu_best);

disp('cov_best=');

disp(cov_best);

disp('f=');

disp(f);

figure(3),

hold on

plot(x(1,:),x(2,:),'.');

plot(mu_real(1,:),mu_real(2,:),'*r'); plot(mu_best(1,:),mu_best(2,:),'+c'); hold off

遗传算法MATLAB完整代码(不用工具箱)

遗传算法解决简单问题 %主程序：用遗传算法求解y=200*exp(-0.05*x).*sin(x)在区间[-2,2]上的最大值clc; clear all; close all; global BitLength global boundsbegin global boundsend bounds=[-2,2]; precision=0.0001; boundsbegin=bounds(:,1); boundsend=bounds(:,2); %计算如果满足求解精度至少需要多长的染色体 BitLength=ceil(log2((boundsend-boundsbegin)'./precision)); popsize=50; %初始种群大小 Generationmax=12; %最大代数 pcrossover=0.90; %交配概率 pmutation=0.09; %变异概率 %产生初始种群 population=round(rand(popsize,BitLength)); %计算适应度，返回适应度Fitvalue和累计概率cumsump [Fitvalue,cumsump]=fitnessfun(population); Generation=1; while Generation

聚类分析Matlab程序实现

2. Matlab程序 2.1 一次聚类法 X=[11978 12.5 93.5 31908;…;57500 67.6 238.0 15900]; T=clusterdata(X,0.9) 2.2 分步聚类 Step1 寻找变量之间的相似性 用pdist函数计算相似矩阵，有多种方法可以计算距离，进行计算之前最好先将数据用zscore 函数进行标准化。 X2=zscore(X); %标准化数据 Y2=pdist(X2); %计算距离 Step2 定义变量之间的连接 Z2=linkage(Y2); Step3 评价聚类信息 C2=cophenet(Z2,Y2); //0.94698 Step4 创建聚类，并作出谱系图 T=cluster(Z2,6); H=dendrogram(Z2); Matlab提供了两种方法进行聚类分析。 一种是利用 clusterdata函数对样本数据进行一次聚类，其缺点为可供用户选择的面较窄，不能更改距离的计算方法； 另一种是分步聚类：（1）找到数据集合中变量两两之间的相似性和非相似性，用pdist函数计算变量之间的距离；（2）用 linkage函数定义变量之间的连接；（3）用 cophenetic函数评价聚类信息；（4）用cluster函数创建聚类。 1．Matlab中相关函数介绍 1.1 pdist函数 调用格式：Y=pdist(X,’metric’) 说明：用‘metric’指定的方法计算 X 数据矩阵中对象之间的距离。’ X：一个m×n的矩阵，它是由m个对象组成的数据集，每个对象的大小为n。 metric’取值如下： ‘euclidean’：欧氏距离（默认）；‘seuclidean’：标准化欧氏距离； ‘mahalanobis’：马氏距离；‘cityblock’：布洛克距离； ‘minkowski’：明可夫斯基距离；‘cosine’： ‘correlation’：‘hamming’： ‘jaccard’：‘chebychev’：Chebychev距离。 1.2 squareform函数 调用格式：Z=squareform(Y,..) 说明：强制将距离矩阵从上三角形式转化为方阵形式，或从方阵形式转化为上三角形式。 1.3 linkage函数 调用格式：Z=linkage(Y,’method’) 说明：用‘method’参数指定的算法计算系统聚类树。 Y：pdist函数返回的距离向量；

PID算法Matlab仿真程序和C程序

增量式PID控制算法Matlab仿真程序设一被控对象G（s）=50/(0.125s^2+7s),用增量式PID控制算法编写仿真程序（输入分别为单位阶跃、正弦信号，采样时间为1ms，控制器输出限幅：[-5,5],仿真曲线包括系统输出及误差曲线，并加上注释、图例）。程序如下clear all; close all; ts=0.001; sys=tf(50,[0.125,7, 0]); dsys=c2d(sys,ts,'z'); [num,den]=tfdata(dsys,'v'); u_1=0.0;u_2=0.0; y_1=0.0;y_2=0.0; x=[0,0,0]'; error_1=0; error_2=0; for k=1:1:1000 time(k)=k*ts; S=2; if S==1 kp=10;ki=0.1;kd=15; rin(k)=1; % Step Signal elseif S==2 kp=10;ki=0.1;kd=15; %Sin e Signal rin(k)=0.5*sin(2*pi*k*ts); end du(k)=kp*x(1)+kd*x(2)+ki*x(3); % PID Controller u(k)=u_1+du(k); %Restricting the output of controller if u(k)>=5 u(k)=5; end if u(k)<=-5 u(k)=-5; end %Linear model yout(k)=-den(2)*y_1-den(3)*y_2+nu m(2)*u_1+num(3)*u_2; error(k)=rin(k)-yout(k); %Return of parameters u_2=u_1;u_1=u(k); y_2=y_1;y_1=yout(k); x(1)=error(k)-error_1; %C alculating P x(2)=error(k)-2*error_1+error_2; %Calculating D x(3)=error(k); %Calculating I error_2=error_1; error_1=error(k); end figure(1); plot(time,rin,'b',time,yout,'r'); xlabel('time(s)'),ylabel('rin,yout'); figure(2); plot(time,error,'r') xlabel('time(s)');ylabel('error'); 微分先行PID算法Matlab仿真程序%PID Controler with differential in advance clear all; close all; ts=20; sys=tf([1],[60,1],'inputdelay',80); dsys=c2d(sys,ts,'zoh'); [num,den]=tfdata(dsys,'v'); u_1=0;u_2=0;u_3=0;u_4=0;u_5=0;

MATLAB课程遗传算法实验报告及源代码

硕士生考查课程考试试卷 考试科目： 考生姓名：考生学号： 学院：专业： 考生成绩： 任课老师(签名) 考试日期：年月日午时至时

《MATLAB 教程》试题： A 、利用MATLA B 设计遗传算法程序，寻找下图11个端点最短路径，其中没有连接端点表示没有路径。要求设计遗传算法对该问题求解。 a e h k B 、设计遗传算法求解f (x)极小值，具体表达式如下： 321231(,,)5.12 5.12,1,2,3i i i f x x x x x i =?=???-≤≤=? ∑ 要求必须使用m 函数方式设计程序。 C 、利用MATLAB 编程实现：三名商人各带一个随从乘船渡河，一只小船只能容纳二人，由他们自己划行，随从们密约，在河的任一岸，一旦随从的人数比商人多，就杀人越货，但是如何乘船渡河的大权掌握在商人手中，商人们怎样才能安全渡河？ D 、结合自己的研究方向选择合适的问题，利用MATLAB 进行实验。 以上四题任选一题进行实验，并写出实验报告。

选择题目： B 、设计遗传算法求解f (x)极小值，具体表达式如下： 321231(,,)5.12 5.12,1,2,3i i i f x x x x x i =?=???-≤≤=? ∑ 要求必须使用m 函数方式设计程序。 一、问题分析（10分） 这是一个简单的三元函数求最小值的函数优化问题，可以利用遗传算法来指导性搜索最小值。实验要求必须以matlab 为工具，利用遗传算法对问题进行求解。 在本实验中，要求我们用M 函数自行设计遗传算法，通过遗传算法基本原理，选择、交叉、变异等操作进行指导性邻域搜索，得到最优解。 二、实验原理与数学模型（20分） （1）试验原理： 用遗传算法求解函数优化问题，遗传算法是模拟生物在自然环境下的遗传和进化过程而形成的一种自适应全局优化概率搜索方法。其采纳了自然进化模型，从代表问题可能潜在解集的一个种群开始，种群由经过基因编码的一定数目的个体组成。每个个体实际上是染色体带有特征的实体；初始种群产生后，按照适者生存和优胜劣汰的原理，逐代演化产生出越来越好的解：在每一代，概据问题域中个体的适应度大小挑选个体；并借助遗传算子进行组合交叉和主客观变异，产生出代表新的解集的种群。这一过程循环执行，直到满足优化准则为止。最后，末代个体经解码，生成近似最优解。基于种群进化机制的遗传算法如同自然界进化一样，后生代种群比前生代更加适应于环境，通过逐代进化，逼近最优解。 遗传算法是一种现代智能算法，实际上它的功能十分强大，能够用于求解一些难以用常规数学手段进行求解的问题，尤其适用于求解多目标、多约束，且目标函数形式非常复杂的优化问题。但是遗传算法也有一些缺点，最为关键的一点，即没有任何理论能够证明遗传算法一定能够找到最优解，算法主要是根据概率论的思想来寻找最优解。因此，遗传算法所得到的解只是一个近似解，而不一定是最优解。 （2）数学模型 对于求解该问题遗传算法的构造过程： （1）确定决策变量和约束条件；

MATLAB实现FCM 聚类算法

本文在阐述聚类分析方法的基础上重点研究FCM 聚类算法。FCM 算法是一种基于划分的聚类算法，它的思想是使得被划分到同一簇的对象之间相似度最大，而不同簇之间的相似度最小。最后基于MATLAB实现了对图像信息的聚类。 第 1 章概述 聚类分析是数据挖掘的一项重要功能，而聚类算法是目前研究的核心，聚类分析就是使用聚类算法来发现有意义的聚类，即“物以类聚” 。虽然聚类也可起到分类的作用，但和大多数分类或预测不同。大多数分类方法都是演绎的，即人们事先确定某种事物分类的准则或各类别的标准，分类的过程就是比较分类的要素与各类别标准，然后将各要素划归于各类别中。确定事物的分类准则或各类别的标准或多或少带有主观色彩。 为获得基于划分聚类分析的全局最优结果，则需要穷举所有可能的对象划分，为此大多数应用采用的常用启发方法包括：k-均值算法，算法中的每一个聚类均用相应聚类中对象的均值来表示；k-medoid 算法，算法中的每一个聚类均用相应聚类中离聚类中心最近的对象来表示。这些启发聚类方法在分析中小规模数据集以发现圆形或球状聚类时工作得很好，但当分析处理大规模数据集或复杂数据类型时效果较差，需要对其进行扩展。 而模糊C均值(Fuzzy C-means, FCM)聚类方法，属于基于目标函数的模糊聚类算法的范畴。模糊C均值聚类方法是基于目标函数的模糊聚类算法理论中最为完善、应用最为广泛的一种算法。模糊c均值算法最早从硬聚类目标函数的优化中导出的。为了借助目标函数法求解聚类问题，人们利用均方逼近理论构造了带约束的非线性规划函数，以此来求解聚类问题，从此类内平方误差和WGSS(Within-Groups Sum of Squared Error)成为聚类目标函数的普遍形式。随着模糊划分概念的提出，Dunn [10] 首先将其推广到加权WGSS 函数，后来由Bezdek 扩展到加权WGSS 的无限族，形成了FCM 聚类算法的通用聚类准则。从此这类模糊聚类蓬勃发展起来，目前已经形成庞大的体系。 第 2 章聚类分析方法 2-1 聚类分析 聚类分析就是根据对象的相似性将其分群，聚类是一种无监督学习方法，它不需要先验的分类知识就能发现数据下的隐藏结构。它的目标是要对一个给定的数据集进行划分，这种划分应满足以下两个特性：①类内相似性：属于同一类的数据应尽可能相似。②类间相异性：属于不同类的数据应尽可能相异。图2.1是一个简单聚类分析的例子。

最短路径算法_matlab程序[1]

算法描述： 输入图G，源点v0，输出源点到各点的最短距离D 中间变量v0保存当前已经处理到的顶点集合，v1保存剩余的集合 1.初始化v1,D 2.计算v0到v1各点的最短距离，保存到D for each i in v0;D(j)=min[D(j),G(v0(1),i)+G(i,j)] ,where j in v1 3.将D中最小的那一项加入到v0，并且从v1删除这一项。 4.转到2，直到v0包含所有顶点。 %dijsk最短路径算法 clear,clc G=[ inf inf 10 inf 30 100; inf inf 5 inf inf inf; inf 5 inf 50 inf inf; inf inf inf inf inf 10; inf inf inf 20 inf 60; inf inf inf inf inf inf; ]; %邻接矩阵 N=size(G,1); %顶点数 v0=1; %源点 v1=ones(1,N); %除去原点后的集合 v1(v0)=0; %计算和源点最近的点 D=G(v0,:); while 1 D2=D; for i=1:N if v1(i)==0 D2(i)=inf; end end D2 [Dmin id]=min(D2); if isinf(Dmin),error,end v0=[v0 id] %将最近的点加入v0集合，并从v1集合中删除 v1(id)=0; if size(v0,2)==N,break;end %计算v0（1）到v1各点的最近距离 fprintf('计算v0（1）到v1各点的最近距离\n');v0,v1 id=0; for j=1:N %计算到j的最近距离 if v1(j)

遗传算法Matlab程序

% f(x)=11*sin(6x)+7*cos(5x),0<=x<=2*pi; %%初始化参数 L=16;%编码为16位二进制数 N=32;%初始种群规模 M=48;%M个中间体，运用算子选择出M/2对母体，进行交叉；对M个中间体进行变异 T=100;%进化代数 Pc=0.8;%交叉概率 Pm=0.03;%%变异概率 %%将十进制编码成16位的二进制，再将16位的二进制转成格雷码 for i=1:1:N x1(1,i)= rand()*2*pi; x2(1,i)= uint16(x1(1,i)/(2*pi)*65535); grayCode(i,:)=num2gray(x2(1,i),L); end %% 开始遗传算子操作 for t=1:1:T y1=11*sin(6*x1)+7*cos(5*x1); for i=1:1:M/2 [a,b]=min(y1);%找到y1中的最小值a，及其对应的编号b grayCodeNew(i,:)=grayCode(b,:);%将找到的最小数放到grayCodeNew中grayCodeNew(i+M/2,:)=grayCode(b,:);%与上面相同就可以有M/2对格雷码可以作为母体y1(1,b)=inf;%用来排除已找到的最小值 end for i=1:1:M/2 p=unidrnd(L);%生成一个大于零小于L的数，用于下面进行交叉的位置if rand()

matlab实现Kmeans聚类算法

题目：matlab实现Kmeans聚类算法 姓名吴隆煌 学号41158007

背景知识 1.简介： Kmeans算法是一种经典的聚类算法，在模式识别中得到了广泛的应用，基于Kmeans的变种算法也有很多，模糊Kmeans、分层Kmeans 等。 Kmeans和应用于混合高斯模型的受限EM算法是一致的。高斯混合模型广泛用于数据挖掘、模式识别、机器学习、统计分析。Kmeans 的迭代步骤可以看成E步和M步，E：固定参数类别中心向量重新标记样本，M：固定标记样本调整类别中心向量。K均值只考虑（估计）了均值，而没有估计类别的方差，所以聚类的结构比较适合于特征协方差相等的类别。 Kmeans在某种程度也可以看成Meanshitf的特殊版本，Meanshift 是一种概率密度梯度估计方法（优点：无需求解出具体的概率密度，直接求解概率密度梯度。），所以Meanshift可以用于寻找数据的多个模态（类别），利用的是梯度上升法。在06年的一篇CVPR文章上，证明了Meanshift方法是牛顿拉夫逊算法的变种。Kmeans 和EM算法相似是指混合密度的形式已知（参数形式已知）情况下，利用迭代方法，在参数空间中搜索解。而Kmeans和Meanshift相似是指都是一种概率密度梯度估计的方法，不过是Kmean选用的是特殊的核函数（uniform kernel），而与混合概率密度形式是否已知无关，是一种梯度求解方式。 k-means是一种聚类算法，这种算法是依赖于点的邻域来决定哪些

点应该分在一个组中。当一堆点都靠的比较近，那这堆点应该是分到同一组。使用k-means，可以找到每一组的中心点。 当然，聚类算法并不局限于2维的点，也可以对高维的空间（3维，4维，等等）的点进行聚类，任意高维的空间都可以。 上图中的彩色部分是一些二维空间点。上图中已经把这些点分组了，并使用了不同的颜色对各组进行了标记。这就是聚类算法要做的事情。 这个算法的输入是： 1：点的数据（这里并不一定指的是坐标，其实可以说是向量） 2：K，聚类中心的个数（即要把这一堆数据分成几组） 所以，在处理之前，你先要决定将要把这一堆数据分成几组，即聚成几类。但并不是在所有情况下，你都事先就能知道需要把数据聚成几类的。但这也并不意味着使用k-means就不能处理这种情况，下文中会有讲解。 把相应的输入数据，传入k-means算法后，当k-means算法运行完后，该算法的输出是： 1：标签（每一个点都有一个标签，因为最终任何一个点，总会被分到某个类，类的id号就是标签） 2：每个类的中心点。 标签，是表示某个点是被分到哪个类了。例如，在上图中，实际上

遗传算法经典MATLAB代码资料讲解

遗传算法经典学习Matlab代码 遗传算法实例: 也是自己找来的，原代码有少许错误，本人都已更正了，调试运行都通过了的。 对于初学者，尤其是还没有编程经验的非常有用的一个文件 遗传算法实例 % 下面举例说明遗传算法% % 求下列函数的最大值% % f(x)=10*sin(5x)+7*cos(4x) x∈[0,10]% % 将x 的值用一个10位的二值形式表示为二值问题，一个10位的二值数提供的分辨率是每为(10-0)/(2^10-1)≈0.01。% % 将变量域[0,10] 离散化为二值域[0,1023], x=0+10*b/1023, 其 中 b 是[0,1023] 中的一个二值数。% % % %--------------------------------------------------------------------------------------------------------------% %--------------------------------------------------------------------------------------------------------------% % 编程 %----------------------------------------------- % 2.1初始化(编码) % initpop.m函数的功能是实现群体的初始化，popsize表示群体的大小，chromlength表示染色体的长度(二值数的长度)，

% 长度大小取决于变量的二进制编码的长度(在本例中取10位)。 %遗传算法子程序 %Name: initpop.m %初始化 function pop=initpop(popsize,chromlength) pop=round(rand(popsize,chromlength)); % rand随机产生每个单元 为{0,1} 行数为popsize，列数为chromlength的矩阵， % roud对矩阵的每个单元进行圆整。这样产生的初始种群。 % 2.2 计算目标函数值 % 2.2.1 将二进制数转化为十进制数(1) %遗传算法子程序 %Name: decodebinary.m %产生[2^n 2^(n-1) ... 1] 的行向量，然后求和，将二进制转化为十进制 function pop2=decodebinary(pop) [px,py]=size(pop); %求pop行和列数 for i=1:py pop1(:,i)=2.^(py-i).*pop(:,i); end pop2=sum(pop1,2); %求pop1的每行之和 % 2.2.2 将二进制编码转化为十进制数(2) % decodechrom.m函数的功能是将染色体(或二进制编码)转换为十进制，参数spoint表示待解码的二进制串的起始位置

最短距离聚类的matlab实现-1(含聚类图-含距离计算)

最短距离聚类的matlab实现-1 【2013-5-21更新】 说明：正文中命令部分可以直接在Matlab中运行， 作者(Yangfd09)于2013-5-21 19:15:50在MATLAB R2009a(7.8.0.347)中运行通过 %最短距离聚类（含距离计算，含聚类图） %说明：此程序的优点在于每一步都是自己编写的，很少用matlab现成的指令， %所以更适合于初学者，有助于理解各种标准化方法和距离计算方法。 %程序包含了极差标准化（两种方法）、中心化、标准差标准化、总和标准化和极大值标准化等标准化方法， %以及绝对值距离、欧氏距离、明科夫斯基距离和切比雪夫距离等距离计算方法。 %==========================>>导入数据<<============================== %变量名为test（新建一个以test变量，双击进入Variable Editor界面，将数据复制进去即可）%数据要求：m行n列，m为要素个数，n为区域个数（待聚类变量）。 % 具体参见末页测试数据。 testdata=test; %============================>>标准化<<=============================== %变量初始化，m用来寻找每行的最大值，n找最小值，s记录每行数据的和 [M,N]=size(testdata);m=zeros(1,M);n=9999*ones(1,M);s=zeros(1,M);eq=zeros(1,M); %为m、n和s赋值 for i=1:M for j=1:N if testdata(i,j)>=m(i) m(i)=testdata(i,j); end if testdata(i,j)<=n(i) n(i)=testdata(i,j); end s(i)=s(i)+testdata(i,j); end eq(i)=s(i)/N; end %sigma0是离差平方和，sigma是标准差 sigma0=zeros(M); for i=1:M for j=1:N sigma0(i)=sigma0(i)+(testdata(i,j)-eq(i))^2; end end sigma=sqrt(sigma0/N);

数学实验05聚类分析---用matlab做聚类分析

用matlab做聚类分析 Matlab提供了两种方法进行聚类分析。 一种是利用clusterdata函数对样本数据进行一次聚类，其缺点为可供用户选择的面较窄，不能更改距离的计算方法； 另一种是分步聚类：（1）找到数据集合中变量两两之间的相似性和非相似性，用pdist函数计算变量之间的距离；（2）用linkage函数定义变量之间的连接；（3）用cophenetic函数评价聚类信息；（4）用cluster函数创建聚类。1．Matlab中相关函数介绍 1.1pdist函数 调用格式：Y=pdist(X,’metric’) 说明：用‘metric’指定的方法计算X数据矩阵中对象之间的距离。’X：一个m×n的矩阵，它是由m个对象组成的数据集，每个对象的大小为n。 metric’取值如下： ‘euclidean’：欧氏距离（默认）；‘seuclidean’：标准化欧氏距离； ‘mahalanobis’：马氏距离；‘cityblock’：布洛克距离； ‘minkowski’：明可夫斯基距离；‘cosine’： ‘correlation’：‘hamming’： ‘jaccard’：‘chebychev’：Chebychev距离。 1.2squareform函数 调用格式：Z=squareform(Y,..)

说明：强制将距离矩阵从上三角形式转化为方阵形式，或从方阵形式转化为上三角形式。 1.3linkage函数 调用格式：Z=linkage(Y,’method’) 说明：用‘method’参数指定的算法计算系统聚类树。 Y：pdist函数返回的距离向量； method：可取值如下： ‘single’：最短距离法（默认）；‘complete’：最长距离法； ‘average’：未加权平均距离法；‘weighted’：加权平均法； ‘centroid’：质心距离法；‘median’：加权质心距离法； ‘ward’：内平方距离法（最小方差算法） 返回：Z为一个包含聚类树信息的（m-1）×3的矩阵。 1.4dendrogram函数 调用格式：[H，T，…]=dendrogram(Z,p，…) 说明：生成只有顶部p个节点的冰柱图（谱系图）。 1.5cophenet函数 调用格式：c=cophenetic(Z,Y) 说明：利用pdist函数生成的Y和linkage函数生成的Z计算cophenet相关系数。 1.6cluster函数 调用格式：T=cluster(Z,…) 说明：根据linkage函数的输出Z创建分类。

最优化算法-Matlab程序

CG程序代码 function [x,y] = cg(A,b,x0) %%%%%%%%%%%%%%%%%CG算法%%%%%%%%%%%% r0 = A*x0-b; p0 = -r0; k = 0; r = r0; p = p0; x = x0; while r~=0 alpha = -r'*p/(p'*A*p); x = x+alpha*p; rold = r; r = rold+alpha*A*p; beta = r'*r/(rold'*rold); p = -r+beta*p; plot(k,norm(p),'.--'); hold on k = k+1; end y.funcount = k; y.fval = x'*A*x/2-b'*x;

function [x,y] = cg_FR(fun,dfun,x0) %%%%%%%%%%%%%%%CG_FR算法%%%%%%%%%%%%%%% error = 10^-5; f0 = feval(fun,x0); df0 = feval(dfun,x0); p0 = -df0; f = f0; df = df0; p = p0; x = x0; k = 0; while ((norm(df)>error)&&(k<1000)) f = feval(fun,x); [alpha,funcNk,exitflag] = lines(fun,0.01,0.15,0.85,6,f,df'*p,x,p);%%用线搜索找下降距离%% if exitflag == -1 disp('Break!!!'); break; end x = x+alpha*p; dfold = df; df = feval(dfun,x); beta = df'*df/(dfold'*dfold); p = -df+beta*p; plot(k,norm(df),'.--'); hold on k = k+1; end y.funcount = k; y.fval = feval(fun,x); y.error = norm(df);

MVDR算法matlab程序

clc clear all close all %% 常量定义 Freqs=1.6e9; %工作频率 c=3e8; %光速 lamda=c/Freqs; %波长 d=0.5*lamda; %单元间距 M=16; %天线阵元数 fs=2e6; %采样频率 pd=10; %快拍数 %% 模型建立 %--------------第一个干扰模型-------------------- thetaJ1=20*pi/180; %干扰方向 FreqJ1=5e5; %第一个干扰的频率 J1NR=sqrt(10^(60/10)); J1one=J1NR*exp(j*(2*pi*FreqJ1*(1:1:pd)/fs)); %1*pd %--------------第二个干扰模型-------------------- ThetaJ2=60*pi/180; %干扰方向 FreqJ2=6e5; %第二个干扰的频率 J2NR=sqrt(10^(60/10)); J2one=J2NR*exp(j*(2*pi*FreqJ2*(1:1:pd)/fs)); %1*pd %--------------信号模型-------------------- ThetaS=30*pi/180; FreqS=7e5; SNR=sqrt(10^(40/10)); Sone=SNR*exp(j*(2*pi*FreqS*(1:1:pd)/fs)); %1*pd %--------------空域处理------------------------- I1=zeros(M,1); I2=zeros(M,1); IS=zeros(M,1); for n=1:M I1(n)=exp(j*2*pi*(n-1)*d*sin(thetaJ1)/lamda); I2(n)=exp(j*2*pi*(n-1)*d*sin(ThetaJ2)/lamda); IS(n)=exp(j*2*pi*(n-1)*d*sin(ThetaS)/lamda); end J1=I1*J1one; J1=J1.'; J2=I2*J2one; J2=J2.'; %------------产生噪声------------------------- noise=sqrt(1/2)*randn(pd,M)+j*sqrt(1/2)*randn(pd,M);

基于遗传算法的matlab源代码

function youhuafun D=code; N=50;%Tunable maxgen=50;%Tunable crossrate=0.5;%Tunable muterate=0.08;%Tunable generation=1; num=length(D); fatherrand=randint(num,N,3); score=zeros(maxgen,N); while generation<=maxgen ind=randperm(N-2)+2;%随机配对交叉 A=fatherrand(:,ind(1:(N-2)/2)); B=fatherrand(:,ind((N-2)/2+1:end)); %多点交叉 rnd=rand(num,(N-2)/2); ind=rnd tmp=A(ind); A(ind)=B(ind); B(ind)=tmp; %%两点交叉 %for kk=1:(N-2)/2 %rndtmp=randint(1,1,num)+1; %tmp=A(1:rndtmp,kk); %A(1:rndtmp,kk)=B(1:rndtmp,kk); %B(1:rndtmp,kk)=tmp; %end fatherrand=[fatherrand(:,1:2),A,B]; %变异 rnd=rand(num,N); ind=rnd[m,n]=size(ind); tmp=randint(m,n,2)+1; tmp(:,1:2)=0; fatherrand=tmp+fatherrand; fatherrand=mod(fatherrand,3); %fatherrand(ind)=tmp; %评价、选择 scoreN=scorefun(fatherrand,D);%求得N个个体的评价函数 score(generation,:)=scoreN; [scoreSort,scoreind]=sort(scoreN); sumscore=cumsum(scoreSort); sumscore=sumscore./sumscore(end); childind(1:2)=scoreind(end-1:end); for k=3:N tmprnd=rand; tmpind=tmprnd difind=[0,diff(t mpind)]; if~any(difind) difind(1)=1; end childind(k)=scoreind(logical(difind)); end fatherrand=fatherrand(:,childind); generation=generation+1; end %score maxV=max(score,[],2); minV=11*300-maxV; plot(minV,'*');title('各代的目标函数值'); F4=D(:,4); FF4=F4-fatherrand(:,1); FF4=max(FF4,1); D(:,5)=FF4; save DData D function D=code load youhua.mat %properties F2and F3 F1=A(:,1); F2=A(:,2); F3=A(:,3); if(max(F2)>1450)||(min(F2)<=900) error('DATA property F2exceed it''s range (900,1450]') end %get group property F1of data,according to F2value F4=zeros(size(F1)); for ite=11:-1:1 index=find(F2<=900+ite*50); F4(index)=ite; end D=[F1,F2,F3,F4]; function ScoreN=scorefun(fatherrand,D) F3=D(:,3); F4=D(:,4); N=size(fatherrand,2); FF4=F4*ones(1,N); FF4rnd=FF4-fatherrand; FF4rnd=max(FF4rnd,1); ScoreN=ones(1,N)*300*11; %这里有待优化

图论算法及matlab程序的三个案例

图论实验三个案例 单源最短路径问题 Dijkstra 算法 Dijkstra 算法是解单源最短路径问题的一个贪心算法。其基本思想是，设置一个顶点集合S 并不断地作贪心选择来扩充这个集合。一个顶点属于集合S 当且仅当从源到该顶点的最短路径长度已知。设v 是图中的一个顶点，记()l v 为顶点 v 到源点v 1的最短距离， ,i j v v V ?∈，若 (,)i j v v E ?，记i v 到j v 的权ij w =∞。 Dijkstra 算法： ① 1{}S v =,1()0l v =；1{}v V v ??-,()l v =∞,1i =,1{}S V v =-； ② S φ=，停止，否则转③； ③ ()min{(),(,)} j l v l v d v v =， j v S ∈，v S ?∈； ④ 存在 1 i v +，使 1()min{()} i l v l v +=，v S ∈； ⑤ 1{} i S S v +=， 1{} i S S v +=-，1i i =+，转②； 实际上，Dijkstra 算法也是最优化原理的应用：如果12 1n n v v v v -是从1v 到 n v 的最短路径，则 12 1 n v v v -也必然是从1v 到 1 n v -的最优路径。 在下面的MATLAB 实现代码中，我们用到了距离矩阵，矩阵第i 行第j 行元 素表示顶点i v 到j v 的权ij w ，若i v 到j v 无边，则realmax ij w =，其中realmax 是 MATLAB 常量，表示最大的实数+308)。 function re=Dijkstra(ma)

遗传算法的MATLAB程序实例

遗传算法的程序实例 如求下列函数的最大值 f(x)=10*sin(5x)+7*cos(4x) x∈[0,10] 一、初始化(编码) initpop.m函数的功能是实现群体的初始化，popsize表示群体的大小，chromlength表示染色体的长度(二值数的长度)， 长度大小取决于变量的二进制编码的长度(在本例中取10位)。 代码: %Name: initpop.m %初始化 function pop=initpop(popsize,chromlength) pop=round(rand(popsize,chromlength)); % rand随机产生每个单元为 {0,1} 行数为popsize，列数为chromlength的矩阵， % roud对矩阵的每个单元进行圆整。这样产生的初始种群。 二、计算目标函数值 1、将二进制数转化为十进制数(1) 代码: %Name: decodebinary.m %产生 [2^n 2^(n-1) ... 1] 的行向量，然后求和，将二进制转化为十进制 function pop2=decodebinary(pop) [px,py]=size(pop); %求pop行和例数 for i=1:py pop1(:,i)=2.^(py-1).*pop(:,i); py=py-1; end pop2=sum(pop1,2); %求pop1的每行之和 2、将二进制编码转化为十进制数(2) decodechrom.m函数的功能是将染色体(或二进制编码)转换为十进制，参数spoint表示待解码的二进制串的起始位置。(对于多个变量而言，如有两个变量，采用20为表示，每个变量10为，则第一个变量从1开始，另一个变量从11开始。本例为1)，参数1ength表示所截取的长度（本例为10）。 代码: %Name: decodechrom.m %将二进制编码转换成十进制 function pop2=decodechrom(pop,spoint,length) pop1=pop(:,spoint:spoint+length-1); pop2=decodebinary(pop1); 3、计算目标函数值 calobjvalue.m函数的功能是实现目标函数的计算，其公式采用本文示例仿真，可根据不同优化问题予以修改。

matlab实现Kmeans聚类算法

matlab实现Kmeans聚类算法 1.简介： Kmeans和应用于混合高斯模型的受限EM算法是一致的。高斯混合模型广泛用于数据挖掘、模式识别、机器学习、统计分析。Kmeans 的迭代步骤可以看成E步和M步，E：固定参数类别中心向量重新标记样本，M：固定均值只考虑（估计）了均值，而没有估计类别的方差，所以聚类的结构比较适合于特征协方差相等的类别。 Kmeans在某种程度也可以看成Meanshitf的特殊版本，Meanshift 是所以Meanshift可以用于寻找数据的多个模态（类别），利用的是梯度上升法。在06年的一篇CVPR文章上，证明了Meanshift方法是牛顿拉夫逊算法的变种。Kmeans和EM算法相似是指混合密度的形式已知（参数形式已知）情况下，利用迭代方法，在参数空间中搜索解。而Kmeans和Meanshift相似是指都是一种概率密度梯度估计的方法，不过是Kmean选用的是特殊的核函数（uniform kernel），而与混合概率密度形式是否已知无关，是一种梯度求解方式。 k-means是一种聚类算法，这种算法是依赖于点的邻域来决定哪些点应该分在点，也可以对高维的空间（3维，4维，等等）的点进行聚类，任意高维的空间都可以。 上图中的彩色部分是一些二维空间点。上图中已经把这些点分组了，并使用了不同的颜色对各组进行了标记。这就是聚类算法要做的事情。 这个算法的输入是： 1：点的数据（这里并不一定指的是坐标，其实可以说是向量）

2：K，聚类中心的个数（即要把这一堆数据分成几组） 所以，在处理之前，你先要决定将要把这一堆数据分成几组，即聚成几类。但并不是在所有情况下，你都事先就能知道需要把数据聚成几类的。意味着使用k-means就不能处理这种情况，下文中会有讲解。 把相应的输入数据，传入k-means算法后，当k-means算法运行完后，该算法的输出是： 1：标签（每一个点都有一个标签，因为最终任何一个点，总会被分到某个类，类的id号就是标签） 2：每个类的中心点。 标签，是表示某个点是被分到哪个类了。例如，在上图中，实际上有4中“标签”，每个“标签”使用不同的颜色来表示。所有黄色点我们可以用标签以看出，有3个类离的比较远，有两个类离得比较近，几乎要混合在一起了。 当然，数据集不一定是坐标，假如你要对彩色图像进行聚类，那么你的向量就可以是(b,g,r)，如果使用的是hsv颜色空间，那还可以使用(h,s,v),当然肯定可以有不同的组合例如(b*b,g*r,r*b) ，(h*b,s*g,v*v)等等。 在本文中，初始的类的中心点是随机产生的。如上图的红色点所示，是本文随机产生的初始点。注意观察那两个离得比较近的类，它们几乎要混合在一起，看看算法是如何将它们分开的。 类的初始中心点是随机产生的。算法会不断迭代来矫正这些中心点，并最终得到比较靠5个中心点的距离,选出一个距离最小的(例如该点与第2个中心点的距离是5个距离中最小的),那么该点就归属于该类.上图是点的归类结果示意图. 经过步骤3后,每一个中心center(i)点都有它的”管辖范围”,由于这个中心点不一定是这个管辖范围的真正中心点,所以要重新计算中心点,计算的方法有很多种,最简单的一种是,直接计算该管辖范围内所有点的均值,做为心的中心点new_center(i). 如果重新计算的中心点new_center(i)与原来的中心点center(i)的距离大于一定的阈值（该阈值可以设定），那么认为算法尚未收敛，使用new_center(i)代替center(i)（如图，中心点从红色点