如何用matlab编最大似然法的程序

如何用matlab编最大似然法的程序,第1张

命令

利用mle函数进行参数估计

函数

mle

格式

phat=mle('dist′

,X)%返回用dist指定分布的最大似然估计值

[phat,

pci]=mle('dist′

,X)

%置信度为95

[phat,

pci]=mle('dist′

,X,alpha)%置信度由alpha确定

[phat,

pci]=mle('dist′

,X,apha,pl)%仅用于二项分布芦则兆,pl为试验次数。

说明

dist为分布函数名,如:beta(β分布)、bino(二项分布)等,盯纯X为数据样本,alpha为显著水平α,1-a为置信度。

phat为估计陪租值,pci为置信区间

希望能对你有所帮助!

基于对似然函数L(θ)形式(一般为连乘式且各因式>0)的考虑,求θ的最大似然估枯悉闹计的一般步骤如下:

(1)写出似然函数

总体X为离散型时:

总体X为连续型时:

(2)对似然函数两边取对数有

总体X为离散型时:

总体X为连续型时陆友:

(3)对

求导数并令之为0:

此方程为对数似然方程。解对数似然方程所得,即为未知参数 的最大似然估计值。

最大似然估计函数在采样样本总数趋于无穷的时候达到最小方差(其证明可见于Cramer-Rao lower bound)。当最大似然估计非偏时,等价的,在极限的情况下我们可以称其有没罩最小的均方差。对于独立的观察来说,最大似然估计函数经常趋于正态分布。

以下内容都是基于Christopher P. Randle教授在研究组上交流时的课件整理而来。 

我们都遇到过这样的问题:坐标图中给出一定数量的点,然后给出x值,预测y值,这是我们就需要根据已有的点建立适合模型,确定模型参数。

我们知道一般参数越多,模型就越适合(fit),预测就越准确:

但是同时,参数越多也会增加预测的错误,因为参数值也需要估计,而且很可能出错。

因此选择一个合适的模型,就是一个权衡利弊(trade off)的过程。

模型选择就是要选出能同时使得参数数量最小化,并且使系统发育树的似然值最大的模型。我们在前面说过模型都是嵌套的,在GTR的基础上通过增强假设条件,减少模型参数。比较这些模型最常用的方法就是likelihood ratio test,用来验证参数的增加是否能增加模型的适合度。

likelihood ratio的表汪搭达式如下:

max[p(data|model)]是给定的模型预测的观察数据对应的最大概率(似然值), 的参数数量要比 少。因为似然值数值通常很小,我们一般还是用自然对数表示:

对于嵌套模型来说,上图的数值应该类似卡方分布(χ2 distribution),自由度(k)等同于 和 的自由参数差异。显著性结果证实:模型越复杂,它的模拟结果越好。

最经常用到的,为最大似然法分析选择模型的程序是jModelTest2。在这个程序里,66个嵌套模型会拿来作比较,依据的是分层测试(hierarchical test),测试的项目包括likelihood ratio或者选择Akaike Information Criterion、Bayesian Information Criterion(比较的对象不需要必须是嵌套模型)。首先用likelihood ratio比较最简单的两个模型(JC69和F81),较好的一个会和一个更复杂的模型比较,直到可选择的模型都比较完。

对于系统学来说,最重要的结果就是系统发育树,即到底哪些分类群是单系群?分类群之间的系统发育关系如何?至于哪个模型与分子数据的符合度更好,模型参数的最优值是多少,这些并不是系统学家关心的。简单来说,这些反映的是系统发育的过程,而系统学只关心结果。

但是,我们能说这些参数没有用吗?它们对研究演化过程和机制的人就是很重要的信息。

比如,如茄陵唯果我们发现引入区分转换和颠换的参数能提高模型的适合度,我们就对演化的分子机制有了进一步认识。

再举几个例子:

如果我们想研究一个编码基因是否在净化选择压力下进化(净化选择会抑制碱基替代),那么我们可以比较同义替代速率(该位点碱基替代不改变编码的氨基酸)和非同义替代速率(该位点碱基替代改变编码的氨基酸)。如果这个选择压力是抑制表型的改变(即氨基酸的改变),那么非同义替代的速率应该低于替代随机发生在同义替代和非同义替代位点情况下的速率。

那怎么实现这个比较呢?我们选择两个模型,一个模型认为同义替代和非同颤培义替代位点的速率是相同的,另一个模型中,同义替代和非同义替代位点的速率是不同的。然后进行likelihood ratio test。看哪一个模型的拟合度更好,我们就能评估是发生了净化选择(purifying selection)还是正选择(positive selection,支持氨基酸的变化)。

目前我们讲到的模型都是GTR的嵌套模型,它们有一些通病:只适用于DNA系列数据;没有考虑到DNA序列中的indel(插入/缺失)情况;认为所有性状的在整个系统发育树上的演化速率都是同质的或者说均一的,或者认为用简单的参数I和

欢迎分享,转载请注明来源:内存溢出

原文地址: http://outofmemory.cn/yw/8260730.html

(0)
打赏 微信扫一扫 微信扫一扫 支付宝扫一扫 支付宝扫一扫
上一篇 2023-04-14
下一篇 2023-04-14

发表评论

登录后才能评论

评论列表(0条)

保存