2014-06-24
From [1], soft-decision is better than hard-decision, requiring nearly half the repetition bits for error correction with FRR 10^(-7). Generate 128-bit key with 3900 SRAM cells, only a single measurement. If multiple Maes et al only needs 1536.
[1] Vincent van der Leest, Bart Preneel etc, "Soft Decision Error Correction for Compact Memory-based PUFs using Single Enrollment", 2012 CHES
2014年6月24日星期二
2014年6月11日星期三
Bio-Informatics
1, Motif finding sampling
找n个基因序列的共同motif, 这个motif可能是某段RNA复制的起始点。需要确定的模型参数有,
(1) motif中的概率矩阵
(2) background的概率矩阵
另外,找motif的起始位置也是很困难,涉及到二维采样(n*基因长度)。
Gibbs sampling可以减少采样难度。
首先固定除第1段基因外的其他基因段的motif位置,在第一段基因中找到打分最高的位置,固定此位置;固定除第二段外的其他motif位置,在第二段中中找打分最高的;依次类推。
2,
找n个基因序列的共同motif, 这个motif可能是某段RNA复制的起始点。需要确定的模型参数有,
(1) motif中的概率矩阵
(2) background的概率矩阵
另外,找motif的起始位置也是很困难,涉及到二维采样(n*基因长度)。
Gibbs sampling可以减少采样难度。
首先固定除第1段基因外的其他基因段的motif位置,在第一段基因中找到打分最高的位置,固定此位置;固定除第二段外的其他motif位置,在第二段中中找打分最高的;依次类推。
2,
MCMC算法简介
Other interesting/ ML MAP Bayesian: https://engineering.purdue.edu/kak/Tutorials/
Gibbs Motif finding: http://www.cs.cmu.edu/~ckingsf/bioinfo-lectures/gibbs.pdf
原文链接:百度贴吧 http://tieba.baidu.com/p/1681867134
前天一个应用数学专业的朋友问我一个问题。这种类型的问题在工程上叫做求解反问题(inverse problem)。这种问题在他们看来是非常麻烦的。但是作为一个学统计的人而言,这种问题是再常见不过的了。也许是不同学科对问题的理解不同吧。我这里就引用这个例子介绍一下统计学的基本思想以及MCMC算法。
Gibbs Motif finding: http://www.cs.cmu.edu/~ckingsf/bioinfo-lectures/gibbs.pdf
原文链接:百度贴吧 http://tieba.baidu.com/p/1681867134
前天一个应用数学专业的朋友问我一个问题。这种类型的问题在工程上叫做求解反问题(inverse problem)。这种问题在他们看来是非常麻烦的。但是作为一个学统计的人而言,这种问题是再常见不过的了。也许是不同学科对问题的理解不同吧。我这里就引用这个例子介绍一下统计学的基本思想以及MCMC算法。
原始问题是这样的,有一个正问题,背景似乎是研究岩石断层什么的:
幻灯片中的等式成立(由相关物理知识决定),称为forward problem。其中theta和phi是可以控制的量,称为输入。左侧的量称为输出,因此是输入量的函数。问题是:有颜色的参数是未知的。
幻灯片中的等式成立(由相关物理知识决定),称为forward problem。其中theta和phi是可以控制的量,称为输入。左侧的量称为输出,因此是输入量的函数。问题是:有颜色的参数是未知的。
我们需要通过一组输入和输出的数据来估计参数,这样的问题称为反问题。从统计学角度上看,这个问题属于参数估计问题。
----------------------------------------
1、统计学/贝叶斯统计
我们首先要用统计语言重新描述这个问题:
假设我们有n次试验,对应于theta和phi的不同组合。我们将第i次试验表述为
其中epsilon_i是第i次实验的随机误差。
假定epsilon_i服从独立同分布的正态分布。不妨假设正态分布的均值为0,因为如果不为0均值也可以被截距项A吸收掉。设正态分布方差为sigma^2,未知。
统计学中的一个核心概念称为似然函数(likelihood function),它是观测到的量的概率分布。这里我们观测到R,所以就要写出R的联合概率分布。注意到R_i也是独立的正态分布,因此联合概率密度为
由于概率密度要满足积分为1的正则性条件。因此统计上常常省略密度前面的常数项,然后用正比于表示所缺常数由正则性条件确定。
其中epsilon_i是第i次实验的随机误差。
假定epsilon_i服从独立同分布的正态分布。不妨假设正态分布的均值为0,因为如果不为0均值也可以被截距项A吸收掉。设正态分布方差为sigma^2,未知。
统计学中的一个核心概念称为似然函数(likelihood function),它是观测到的量的概率分布。这里我们观测到R,所以就要写出R的联合概率分布。注意到R_i也是独立的正态分布,因此联合概率密度为
由于概率密度要满足积分为1的正则性条件。因此统计上常常省略密度前面的常数项,然后用正比于表示所缺常数由正则性条件确定。
然后我们这里应用贝叶斯统计的观点,认为参数也是有一定分布的。这是因为在样本量有限的情况下参数没办法准确确定,因此用一个分布来描述是合理的。假设我们给所有参数赋予无信息先验,根据贝叶斯公式,就得到后验密度
贝叶斯统计认为参数的一切信息都包含在后验分布当中。于是我们只要研究这个后验密度函数就可以了。
然而不幸的是,这个后验密度有5维,而且还不知道前面的归一化常数,因此没办法直观把握。于是贝叶斯统计中有了一个重要的边缘化的思想:考虑每个参数的边缘分布,即将其他参数积分掉。尽管贝叶斯统计历史悠久,这个边缘化思想却是近30年来随着Markov Chain Monte Carlo (MCMC)的兴起才开始出现的。
简单地说,就是通过计算形如:
的式子求出每个参数的边缘分布,然后统计推断就结束了。
贝叶斯统计认为参数的一切信息都包含在后验分布当中。于是我们只要研究这个后验密度函数就可以了。
然而不幸的是,这个后验密度有5维,而且还不知道前面的归一化常数,因此没办法直观把握。于是贝叶斯统计中有了一个重要的边缘化的思想:考虑每个参数的边缘分布,即将其他参数积分掉。尽管贝叶斯统计历史悠久,这个边缘化思想却是近30年来随着Markov Chain Monte Carlo (MCMC)的兴起才开始出现的。
简单地说,就是通过计算形如:
的式子求出每个参数的边缘分布,然后统计推断就结束了。
________________________________
*边缘分布(对除theta外其他变量的积分,从而消除其他变量影响)
在介绍MCMC算法之前。我们先讲一下计算这个边缘分布的难点。第一个难点是归一化常数不知道,第二个则是高维数值积分本身的困难。
因为这两个困难,目前认为最有效的算法是Monte Carlo算法。算法的基本思想是,产生服从后验分布的样本,设想如果我们能够产生服从P(A,B_1,B_2,\phi,\sigma^2|R)的样本,那么它的第一维就服从P(A|R),依此类推。如果我能产生N个样本,N足够大,我们就可以用样本的边缘分布近似总体的边缘分布,也就是说,取出样本的第一维,就可以得到A的近似后验分布。因此问题转化为了已知密度函数(缺归一化常数)的抽样问题。
_______________________________
2、MCMC算法
MCMC不是一个算法,而是一类算法的总称。从数学上讲,算法的思想是产生一个Markov Chain,以目标分布为平稳分布。根据Markov Chain的理论,一个Markov Chain从任意初值出发,都会收敛到其平稳分布(如果存在)。注意这个收敛是依分布的。换句话说,如果我模拟了一条这样的马尔科夫链,去除掉前面一部分样本之后,就可以认为后面的样本来自于平稳分布。这样上面所述的抽样问题就解决了。
MCMC与传统的独立抽样的关键区别在于:抽样的结果是相关的。理论上这对最终结果没有影响(数学上称为遍历理论)。然而实际上仍然有一系列的问题。其中最主要的是收敛诊断问题。和传统算法不同,MCMC的结果无法直接看出是否收敛。因此有的人认为使用MCMC需要慎重。但是根据我的经验,MCMC的表现基本还是可以相信的。
MCMC与传统的独立抽样的关键区别在于:抽样的结果是相关的。理论上这对最终结果没有影响(数学上称为遍历理论)。然而实际上仍然有一系列的问题。其中最主要的是收敛诊断问题。和传统算法不同,MCMC的结果无法直接看出是否收敛。因此有的人认为使用MCMC需要慎重。但是根据我的经验,MCMC的表现基本还是可以相信的。
2.1、Metropolis算法
史上第一个MCMC算法诞生于1942年,现在称为Metropolis算法。
这个算法首先要求一个状态空间上的对称转移函数。即设x,y是状态空间上任意两点,则从x转移到y的概率和从y转移到x的概率相等。
Metropolis算法描述如下:
设x是上一次的状态,通过对称转移函数得到新状态x',然后抽取随机变量U服从(0,1)上的均与分布,如果
1)U<=q(x')/q(x),则markov chain转移到x';
2)U>q(x')/q(x),则markov chain保持不变,仍为x,
其中q是抽样的目标的分布。
容易证明这样构造的markov chain以q为平稳分布。
这个算法有一个推广,称为Hasting's Generalization。如果转移函数不是对称的,只需要将q(x')/q(x)改为{q(x')p(x'->x)}/{q(x)p(x->x')},其中p(x'->x)是从x'转移到x的概率。
2.2、Gibbs抽样
Metropolis是一个包打天下的算法。但是和其他一切适用范围广的算法一样,这个算法收敛是很慢的。统计中最常见的MCMC算法称为Gibbs抽样,它可以看做是一种特殊的Metropolis算法。
Metropolis算法的主要问题,是算法效率严重依赖于转移函数的选取。如果选得不好,算法表现会很差。Gibbs抽样则不存在这样的问题,但是它要求被抽样的函数具有一定的形式。Gibbs抽样主要针对高维抽样问题。这种问题如果直接用Metropolis会变得很难。Gibbs抽样的想法是:既然n维抽样很难,我们就把他化成n个一维的抽样。
Metropolis算法的主要问题,是算法效率严重依赖于转移函数的选取。如果选得不好,算法表现会很差。Gibbs抽样则不存在这样的问题,但是它要求被抽样的函数具有一定的形式。Gibbs抽样主要针对高维抽样问题。这种问题如果直接用Metropolis会变得很难。Gibbs抽样的想法是:既然n维抽样很难,我们就把他化成n个一维的抽样。
这里为了简单起见,我们举一个二维的例子,更高的维度可以直接推广:
设目标分布为q(x,y),我们取初值(x_0,y_0),Gibbs抽样一次迭代为
1)抽取x_1服从q(x|y_0),
2)抽取y_1服从q(y|x_1),
这样就得到了(x_1,y_1)。一直重复下去就得到需要的链。
很显然这个算法要求两个条件分布q(x|y)和q(y|x)具有简单的形式,否则这两个一维抽样也是不可以实现的。当然如果一维抽样不能直接实现,还可以通过嵌入Metropolis算法,用一个接受/拒绝的过程代替直接从条件分布中抽样。
和Metropolis算法相比Gibbs抽样少了拒绝的步骤,这是一般情况下Gibbs抽样效果更好的原因。
——————————————
在介绍了两种基本算法之后,我们回到原来的问题。由于这个问题是多维的,我们首先考虑Gibbs抽样。
I、让我们来看看每个参数的条件分布是什么:
1、首先,A在给定B_1,B_2,phi,sigma^2下的分布是什么呢?
答案是正态分布,因为分布的形式为exp{-“A的二次函数”},这是正态分布的核(回顾一下我们看分布的时候是不关心前面的常数系数的,所以看到这样的函数形式,就可以确定它是正态分布了)。事实上,(A,B_1,B_2)这个三维向量在给定phi和sigma^2的条件下服从联合正态分布(核的形式为exp{-“A,B_1,B_2的二次型”})。经过配方(由矩阵运算可以迅速得到),可以算出这个正态分布的均值和协方差矩阵。
2、然后,sigma^2在给定A,B_1,B_2,phi下是什么分布呢?
我们注意到如果把sigma^2换成1/sigma^2就具有gamma分布的形式。所以这个分布被称为inverse gamma分布。
3、最后我们来看phi在给定A,B_1,B_2,sigma^2下的分布。
这里补充一下:原始物理问题有个约束条件,phi取值于0°到45°之间。非常不幸的是,phi的分布不属于任何已知的分布类型。因此对于phi的抽样我们只能借助metropolis算法来实现了。
当然这里要怎么选取转移函数也是一个可以讨论的问题。但是我考虑到phi的取值范围比较小,就使用一种最笨的抽样方式来解决了。
我使用独立的[0°,45°]的均匀分布作为转移。也就是说,不管上一次的结果是什么,下一次的新值都从这个均匀分布中独立的产生。这个转移函数显然是对称的。然后用metropolis定义的概率进行接受或拒绝就行了。这个转移函数的缺点是不能利用上次成功的经验,但是优点是两个点之间相关性会减小(如果接受)。由于这个问题本身的困难程度不大,只要抽样的次数足够多(如抽10W次),就可以比较精确得得到后验分布了。
当然这里要怎么选取转移函数也是一个可以讨论的问题。但是我考虑到phi的取值范围比较小,就使用一种最笨的抽样方式来解决了。
我使用独立的[0°,45°]的均匀分布作为转移。也就是说,不管上一次的结果是什么,下一次的新值都从这个均匀分布中独立的产生。这个转移函数显然是对称的。然后用metropolis定义的概率进行接受或拒绝就行了。这个转移函数的缺点是不能利用上次成功的经验,但是优点是两个点之间相关性会减小(如果接受)。由于这个问题本身的困难程度不大,只要抽样的次数足够多(如抽10W次),就可以比较精确得得到后验分布了。
II、再讲一下抽样完毕以后如何做分析:
首先要砍掉前面没有达到平稳的点,称为burn-in。对于这种简单问题砍掉5000个点已经绰绰有余了。然后对于剩下的点,一维维的画直方图,或者做密度估计,都可以。一维问题总是很容易处理的,通过画图你一定会对你感兴趣的参数获得认识了。
END
--------------------------讨论-------------------
1、请问 怎么证明这个算法的正确性 比如说 他的不可约 平稳分布和极限分布是否相同等等
一个Markov chain收敛条件是不可约、非周期、正常返。这三个条件是容易成立的(只要模型可识别)。至于极限分布是不是想要的分布,只要验证目标分布是这个迭代算法的不动点就行了。对于大多数的问题这种验证都是平凡的,因此在统计文献中一般都省略此步骤。
2、在运用了MCMC后,你文中最后得到的结果应该是待估计参数的边缘后验分布吧,但实际运用中需要估计参数的确切值,这个确切的值是多少?比如A1的值估计出来为多少?
贝叶斯统计认为,参数的所有信息都包含在后验分布当中。同时,通过有限的样本是无法准确估计参数的,因此用后验分布来表示参数的不确定性。
如果需要点估计,可以根据实际情况选取后验均值或众数(posterior mode)。但是这样的方法会忽略参数的不确定性。因此我建议在后续计算中选用可以保留这种不确定性的方法。最近应用数学中有个很火的领域叫uncertainty quantification,与此相关。
如果需要点估计,可以根据实际情况选取后验均值或众数(posterior mode)。但是这样的方法会忽略参数的不确定性。因此我建议在后续计算中选用可以保留这种不确定性的方法。最近应用数学中有个很火的领域叫uncertainty quantification,与此相关。
................
—————MCMC资源————原文链接:http://asc.2dark.org/node/18
google:MCMC tutorial
http://en.wikipedia.org/wiki/Markov_chain_Monte_Carlo
http://www.soe.ucsc.edu/classes/cmps290c/Winter06/paps/mcmc.pdf
http://public.lanl.gov/kmh/talks/maxent00b.pdf
http://www.soe.ucsc.edu/classes/cmps290c/Winter06/paps/mcmc.pdf
http://public.lanl.gov/kmh/talks/maxent00b.pdf
MCMC preprint service:
http://www.statslab.cam.ac.uk/~mcmc/
http://www.statslab.cam.ac.uk/~mcmc/
David MacKay's book(electronic version availiable):
http://www.inference.phy.cam.ac.uk/mackay/itila/
http://www.inference.phy.cam.ac.uk/mackay/itila/
Review by Adrew Liddle:
Statistical methods for cosmological parameter selection and estimation
Statistical methods for cosmological parameter selection and estimation
2014年6月9日星期一
2014年6月4日星期三
RSA notes
1, Does RSA work for any message M?
Yes.
Even if message m and modulus n is not co-prime, ist gcd{m, n} !=1
2,
Common modulus attack on RSA when the 2 public exponents differ by a single bit
3, Common n attack
Multiple parties cannot share the same n. Because for one person, knowing {e1, d1} = knowing Phi{n} = can factor n = can find another d2 according to e2 = know plain text sent to person2
2014年6月1日星期日
Changing Fonts in Matlab Plots
Ref: http://stackoverflow.com/questions/8934468/changing-fonts-in-matlab-plots/8957573#8957573
more setting see HINT/mfiles/hist_custom.m
more setting see HINT/mfiles/hist_custom.m
It's possible to change default fonts, both for the axes and for other text, by adding the following lines to the
startup.m
file.% Change default axes fonts.
set(0,'DefaultAxesFontName', 'Times New Roman')
set(0,'DefaultAxesFontSize', 14)
% Change default text fonts.
set(0,'DefaultTextFontname', 'Times New Roman')
set(0,'DefaultTextFontSize', 14)
If you don't know if you have a
startup.m
file, runwhich startup
to find its location. If Matlab says there isn't one, run
userpath
to know where it should be placed.
订阅:
博文 (Atom)