mcmc算法
Ⅰ mcmc算法是让分布最终达到平稳分布吗
如果一个马尔科夫链满足 细致平稳条件，则其一定是收敛的，也就是会达到上述的平稳分布。注意：细致平稳条件只是马尔科夫链收敛的充分条件，不是必要条件。
细致平稳条件：，也就是从状态 转移到状态 的数量和从状态 转移到状态 的数量相一致，也就相互抵消，所以数量不发生改变。
Ⅱ 如何在matlab中使用metropolis-hasting算法
MH算法也是一种基于模拟的MCMC技术，一个很重要的应用是从给定的概率分布中抽样。主要原理是构造了一个精妙的Markov链，使得该链的稳态 是你给定的概率密度。它的好处，不用多说，自然是可以对付数学形式复杂的概率密度。有人说，单维的MH算法配上Gibbs Sampler几乎是“无敌”了。今天试验的过程中发现，MH算法想用好也还不简单，里面的转移参数设定就不是很好弄。即使用最简单的高斯漂移项，方差的确定也是个头疼的问题，需要不同问题不同对待，多试验几次。当然你也可以始终选择“理想”参数。还是拿上次的混合高斯分布来做模拟，模拟次数为500000次的时候，概率分布逼近的程度如下图。虽然几个明显的"峰"已经出来了，但是数值上还是 有很大差异的。估计是我的漂移方差没有选好。感觉还是inverse sampling好用，迭代次数不用很多，就可以达到相当的逼近程度。
试了一下MH算法，R Code：
p=function(x,u1,sig1,u2,sig2){
(1/3)*(1/(sqrt(2*pi)*15)*exp(-0.5*(x-70)^2/15^2)+1/(sqrt(2*pi)*11)*exp(-0.5*(x+80)^2/11^2)+1/(sqrt(2*pi)*sig1)*exp(-0.5*(x-u1)^2/sig1^2)+1/(sqrt(2*pi)*sig2)*exp(-0.5*(x-u2)^2/sig2^2))
}
MH=function(x0,n){
x=NULL
x[1] = x0
for (i in 1:n){
x_can= x[i]+rnorm(1,0,3.25)
d= p(x_can,10,30,-10,10)/p(x[i],10,30,-10,10)
alpha= min(1,d)
u=runif(1,0,1)
if (u<alpha){
x[i+1]=x_can}
else{
x[i+1]=x[i]
}
if (round(i/100)==i/100) print(i)
}
x
}
z=MH(10,99999)
z=z[-10000]
a=seq(-100,100,0.2)
plot(density(z),col=1,main='Estimated Density',ylim=c(0,0.02),lty=1)
points(a, p(a,10,30,-10,10),pch='.',col=2,lty=2)
legend(60,0.02,c("True","Sim (MH)"),col=c(1,2),lty=c(1,2))
Ⅲ 贝叶斯的mcmc算法怎么在matlab里实现呀 或者gibbs抽样 求代码
MATLAB中只有矩阵和cell，如果是想要在matlab里实现
要取出第r行、第c列的矩阵元素用：A(r,c)
要取出第r行、第c列的cell成员用：A{r,c}
Ⅳ 贝叶斯mcmc的matlab程序怎么写
%By Shelley from NCUT，April 14th 2011
%Email:[email protected]
%此程序利用贝叶斯分类算法，首先对两类样本进行训练，
%进而可在屏幕上任意取点，程序可输出属于第一类，还是第二类
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
close all
%读入两类训练样本数据
load data
%求两类训练样本的均值和方差
u1=mean(Sample1);
u2=mean(Sample2);
sigm1=cov(Sample1);
sigm2=cov(Sample2);
%计算两个样本的密度函数并显示
x=-20:0.5:40;
y= -20:0.5:20;
[X,Y] = meshgrid(x,y);
F1 = mvnpdf([X(:),Y(:)],u1,sigm1);
F2 = mvnpdf([X(:),Y(:)],u2,sigm2);
P1=reshape(F1,size(X));
P2=reshape(F2,size(X));
figure(2)
surf(X,Y,P1)
hold on
surf(X,Y,P2)
shading interp
colorbar
title('条件概率密度函数曲线');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%以下为测试部分
%利用ginput随机选取屏幕上的点（可连续取10个点）
%程序可根据点的位置自动地显示出属于那个类
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
pw1=0.4;pw2=0.6;
figure(1)
plot(Sample1(:,1),Sample1(:,2),'r.')
hold on
plot(Sample2(:,1),Sample2(:,2),'b.')
for i=1:10
[u,v]=ginput(1);
plot(u,v,'m*');
P1=pw1*mvnpdf([u,v],u1,sigm1);
P2=pw2*mvnpdf([u,v],u2,sigm2);
hold all
if(P1>P2)
disp('it belong to the first class');
else
disp('it belong to the second class');
end;
end
Ⅳ 如何理解贝叶斯估计
根据贝叶斯公式，进行统计推断，
在垃圾邮件分类方面应用很广，方法简单，具有很好的稳定性和健壮性
Ⅵ mcmc在r程序中用哪个软件包
%By Shelley from NCUT，April 14th 2011 %Email:[email protected] %此程序利用贝叶斯分类算法，首先对两类样本进行训练， %进而可在屏幕上任意取点，程序可输出属于第一类，还是第二类 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear; close all %读入两类训练样本数据 load data %求两类训练样本的均值和方差 u1=mean(Sample1); u2=mean(Sample2); sigm1=cov(Sample1); sigm2=cov(Sample2); %计算两个样本的密度函数并显示 x=-20:0.5:40; y= -20:0.5:20; [X,Y] = meshgrid(x,y); F1 = mvnpdf([X(:),Y(:)],u1,sigm1); F2 = mvnpdf([X(:),Y(:)],u2,sigm2); P1=reshape(F1,size(X)); P2=reshape(F2,size(X)); figure(2) surf(X,Y,P1) hold on surf(X,Y,P2) shading interp colorbar title('条件概率密度函数曲线'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %以下为测试部分 %利用ginput随机选取屏幕上的点（可连续取10个点） %程序可根据点的位置自动地显示出属于那个类 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% pw1=0.4;pw2=0.6; figure(1) plot(Sample1(:,1),Sample1(:,2),'r.') hold on plot(Sample2(:,1),Sample2(:,2),'b.') for i=1:10 [u,v]=ginput(1); plot(u,v,'m*'); P1=pw1*mvnpdf([u,v],u1,sigm1); P2=pw2*mvnpdf([u,v],u2,sigm2); hold all if(P1>P2) disp('it belong to the first class'); else disp('it belong to the second class'); end; end
Ⅶ 谁能详细说说Monte Carlo算法的历史
权权的《Monte Carlo方法》系列预告出来时我已注意到，但由于近来时间有限而直到今晚才有时间细读了第一篇《积分方法》。总体而言写得非常清楚，希望能够坚持继续下去。由于整个Bayesian统计学的根基就在Monte Carlo计算法，我对这个领域一向很感兴趣。但由于在研究工作中尚没有机会使用Bayesian统计学，因此有关的知识都还属间接经验，能够在本论坛探讨这个话题肯定会受益。
在点评之前先推荐几本参考资料，我相信下面这个书单是相当不错的，可惜本人尚无时间深入钻研：
* 对英文着作尚有心理障碍者可以参考一本出色的中文教科书：冯康先生所着《数
值计算方法》的第七章《蒙特卡洛方法》（国防工业出版社，1978）；
* 一本可读性极强的英文专着，美国哈佛大学教授Jun Liu所着"Monte Carlo
Strategies in Scientific Computing" (Springer 2002);
* 对Monte Carlo方法在Bayesian统计学中的广泛应用有兴趣者可以适当参考
Andrew Gelman等人所着的"Bayesian Data Analysis" (Second Edition,
2003)之第三部分。
权权将贴子发在物理论坛的目的显然是强调该方法在物理上的运用，我选择在数学论坛加以点评是更看重其统计学背景，着重点各有不同。
>> 蒙特卡洛（Monte Carlo）是摩纳哥公国一个城镇，位于地中海沿岸，以其赌场和豪华
>> 酒店而闻名，所以就有了以随机方法应用于数值计算的一类方法，被称为Monte Carlo
有关Monte Carlo方法历史背景的最精确描述来自Jun Liu的专着，他指出一批物理学家在二战期间为估算薛定谔方程的本征值而发明了一种基于统计抽样的数值计算法，其最初想法归功于Ulam。后来Ulam的同事Metropolis将该方法命名为Monte Carlo。1950年代Metropolis和几名统计物理学同事发表了一篇经典论文，提出了Markov Chain Monte Carlo(MCMC)算法。而MCMC法后来是Bayesian统计学能够不断前进的主要动力。
>> I = ∫ f(x)p(x)dx
这里可以强调一下x是个矢量。而这个积分是概率统计中数学期望的基本定义，可以写成E(f(x))。对于初学者而言，不要忘记概率密度函数p(x)的取值是可以大于1的，归一化条件是对累积密度函数而言。
>> 上述变换就是Monte Carlo积分的基本精神，因为需要用到随机抽样，必然伴随统计误差。
需要用到随机抽样，其动机是想用数值模拟实验中的频率来直接估计一个概率值，而这个概率值是计算许多复杂高维积分的关键。而数值模拟需要产生一个序列的随机数来保证抽样过程的随机性。
>> 因为x_i是按照概率密度p(x)分布的随机变量，f(x_i)也是随机变量
为了论述的清晰，应该说x_i是一个随机矢量，那么f(x_i)就是随机变量（标量）。
>> 而中心极限定理告诉我们，一组独立随机变量之和的概率分布是高斯，其方差等于每一
>> 项随机变量的方差之和
这里关于“中心极限定理”的表述不够精确，容易引起读者混淆，特将Kai-Lai Chung(钟开莱，我国着名数理统计大师许宝禄先生的弟子)概率论教科书中的定义按我的理解方式用英文转述一下：
[Central Limit Theorem] For mutually independent (or weakly correlated) random variables X_1, X_2, ..., X_n with mean mu and variance sigma^2,
√n ( Xbar - mu) / sigma --> N(0,1) in distribution,
where N(0,1) stands for standard Gaussian distribution. This means that the distribution shape of Xbar is more and more like a Gaussian random variable as n increases.
权权的中文表述中漏说了这组随机变量必须来自同一个总体(population)这个重要条件，而且“是高斯”必须改成“在n不断增大时趋向于高斯分布”。
>> 而计算高维积分时，Monte Carlo 方法是较优的选择。
权权只从收敛速度的视角来说明Monte Carlo方法在高维情形下的优越性是不够的，更关键的一点是---Monte Carlo模拟结果的精度和概型的维数D无关！结果的精度显然比收敛速度更为重要，因此Monte Carlo方法特别适合求解高维问题。
另外要指出Monte Carlo方法以O(1/√N)的速度收敛，这在理论上已经无法改善。关键要在实际应用中通过巧妙设计模拟概型和改进抽样方法来降低方差。降低方差的技巧是衡量各种Monte Carlo方法优劣的重要指标。
>> 不妨回顾一下布丰投针实验来结束本篇，在分布着等距平行木纹的地板上投针，要求
>> 针的长度小于木纹之间的距离，几何概形计算结果表明，针与一条木纹相交的概率可
>> 以用针的长度、木纹间距和圆周率π表示。而用几何概形计算概率实质上归结为面积的
>> 计算，也就是积分的计算，布丰投针实验可以说是用随机抽样计算积分的始祖。
Buffon投针实验看似简单，其中蕴含的几何概型思想值得细细品味。令针的长度为L，木纹间距为S，要求L < S。若针的中点到最近的一条平行线的距离为H，用a表示针与平行线的夹角。显然有约束条件0 <= H <= S/2 和 0 <= a <= π。为了使针与平行线相交，必须满足
H <= (L/2) sin(a)
这样针与平行线相交的概率就是两块面积的比值：
p = ∫_0^π (L/2) sin(a) da / (π S/2 ) = 2L / (π S)
这就是权权所说“而用几何概型计算概率实质上归结为面积的计算，也就是积分的计算”。倘若上式分子中的积分是一个复杂的高维积分，我们就可以用Monte Carlo方法模拟出的p值来估算它。当然假如我们感兴趣的是无理数π值的估算，那么由上式可推出：
π_hat = lim 2L / (S p_n)
极限中的n趋向于正无穷。
希望权权在接下来的系列文章中能谈到以下四种Monte Carlo抽样方法：
* Crude Sampling
* Acceptance-Rejection Sampling
* Stratified Sampling
* Importance Sampling
若能谈及MCMC类方法在统计物理学上的运用则更能引人入胜。
Ⅷ bnt matlab 怎么做mcmc有向无环贝叶斯网络结构学习
基于matlab的贝叶斯网络工具箱BNT是kevin p.murphy基于matlab语言开发的关于贝叶斯网络学习的开源包，提供了许多贝叶斯网络学习的底层基础函数库，支持多种类型的节点（概率分布）、精确推理和近似推理、参数学习及结构学习、静态模型和动态模型。
贝叶斯网络表示：BNT中使用矩阵方式表示贝叶斯网络，即若节点i到j有一条弧，则对应矩阵中（i，j）值为1，否则为0。
结构学习算法函数：BNT中提供了较为丰富的结构学习函数，都有：
学习树扩展贝叶斯网络结构的TANC算法learn_struct_tan().
2. 数据完整条件下学习一般贝叶斯网络结构的K2算法learn_struct_k2()、贪婪搜索GS（greedy search）算法learn_struct_gs()和爬山HC（hill climbing）算法learn_struct_hc()等。
3. 缺失数据条件下学习一般贝叶斯网络结构的最大期望EM（expectation maximization）算法learn_struct_EM（）和马尔科夫链蒙特卡罗MCMC（Markov Chain Monte Carlo）learn_struct_mcmc()算法等。
参数学习算法函数：BNT中也提供了丰富的参数学习函数，都有：
1． 完整数据时，学习参数的方法主要有两种：最大似然估计learn_params()和贝叶斯方法bayes_update_params();
2． 数据缺失时，如果已知网络拓扑结构，用EM算法来计算参数，倘若未知网络拓扑结构，使用结构最大期望SEM（structure EM）算法learn_struct_SEM()。
推理机制及推理引擎：为了提高运算速度，使各种推理算法能够有效应用，BNT工具箱采用了引擎机制，不同的引擎根据不同的算法来完成模型转换、细化和求解。这个推理过程如下：
BNT中提供了多种推理引擎，都有：
1． 联合树推理引擎jtree_inf_engine();
2． 全局联合树推理引擎global_joint_inf_engine();
3． 信念传播推理引擎 belprop_inf_engine();
4． 变量消元推理引擎 var_elim_inf_engine().
Ⅸ 贝叶斯后验概率算法有多少种
表达式有无数种
本质是一个
P(A|B)=P(AB)/P(B)
Ⅹ 偏最小二乘法与最小二乘法有什么区别
一、指代不同
1、偏最小二乘法：够在自变量存在严重多重相关性的条件下进行回归建模；允许在样本点个数少于变量个数的条件下进行回归建模。
2、最小二乘法：通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小。
二、特点不同
1、偏最小二乘法：在计算方差和协方差时，求和号前面的系数有两种取法：当样本点集合是随机抽取得到时，应该取1/(n-1)；如果不是随机抽取的，这个系数可取1/n。
2、最小二乘法：可用于曲线拟合。其他一些优化问题也可通过最小化能量或最大化熵用最小二乘法来表达。
三、用法不同
1、偏最小二乘法：在自变量的简单相关系数矩阵中，有某些自变量的相关系数值较大。回归系数的代数符号与专业知识或一般经验相反；或者，它同该自变量与y的简单相关系数符号相反。对重要自变量的回归系数进行t检验，其结果不显着。
2、最小二乘法：如果预测的变量是连续的，我们称其为回归。回归分析中，如果只包括一个自变量和一个因变量，且二者的关系可用一条直线近似表示，这种回归分析称为一元线性回归分析。