标签归档:模拟

漫谈正态分布的生成

本文作者简介:王夜笙,就读于郑州大学信息工程学院,感兴趣的方向为逆向工程和机器学习,长期从事数据抓取工作(长期与反爬虫技术作斗争~),涉猎较广(技艺不精……),详情请见我的个人博客~
个人博客地址:http://bindog.github.io/blog/

感谢怡轩同学的悉心指导~

之前拜读了靳志辉(@rickjin)老师写的《正态分布的前世今生》,一直对正态分布怀着一颗敬畏之心,刚好最近偶然看到python标准库中如何生成服从正态分布随机数的源码,觉得非常有趣,于是又去查找其他一些生成正态分布的方法,与大家分享一下。

利用中心极限定理生成正态分布

设$X_1,X_2,\cdots ,X_n$为独立同分布的随机变量序列,均值为$\mu$,方差为$\sigma^2$,则

$$Z_n=\frac{X_1+X_2+\cdots+X_n-n\mu}{\sigma \sqrt n}$$

具有渐近分布$N(0,1)$,也就是说当$n \rightarrow \infty$时,

$$P\left \{ \frac{X_1+X_2+\cdots+X_n-n\mu}{\sigma \sqrt n} \leq x \right \} \rightarrow \frac{1}{\sqrt{2\pi} } \int_{-\infty }^{x} e^{ -\frac{t^2}{2} } \, dt$$

换句话说,$n$个相互独立同分布的随机变量之和的分布近似于正态分布,$n$越大,近似程度越好。当然也有例外,比如$n$个独立同分布的服从柯西分布随机变量的算术平均数仍是柯西分布,这里就不扩展讲了。

根据中心极限定理,生成正态分布就非常简单粗暴了,直接生成n个独立同分布的均匀分布即可,看代码

继续阅读漫谈正态分布的生成

嘿,朋友,抢红包了吗?

如果你有一台智能手机,如果你装了一个名叫微信的软件,那么你今年的春节很可能是在下面这样的场景中度过的(图片来自微信群):

lucky_money_3 lucky_money_2 lucky_money_1

这也使得众多的网络大V发出了下面的感慨:

lucky_money_4

而最近几天不少微信群里面又流行起来一种“红包接力”的玩法,大概的规则是:群里面先由一人发一个红包,然后大家开始抢,其中金额最大的那个人继续发新一轮的红包,之后不断往复循环。

这时候大家或许就会问了,一直这么玩下去会有什么结果呢?是“闷声赚大钱”了,还是“错过几个亿”了?是最终实现“共同富裕”了,还是变成“寡头垄断”了?要回答这些问题,我们不妨用统计模拟的方法来做一些随机实验,得到的结果或许会让你大跌眼镜呢。

继续阅读嘿,朋友,抢红包了吗?

MCMC案例学习

本文是R中mcmc包的一篇帮助文档,作者为Charles J.Geyer。经过knitr编译后的pdf文档可见此处,提供中文译稿的作者:
闫超,天津财经大学统计系2011级研究生,方向:非寿险准备金评估。
高磊,天津财经大学统计系2011级研究生,方向:非寿险准备金评估。

这个案例,我们不关心题目的具体意义,重点放在利用贝叶斯的观点来解决问题时,MCMC在后续的计算中所发挥的巨大作用。我们知道,贝叶斯的结果往往是一个后验分布。这个后验分布往往很复杂,我们难以用经典的方法求解其期望与方差等一系列的数据特征,这时MCMC来了,将这一系列问题通过模拟来解决。从这个意义上说,MCMC是一种计算手段。依频率学派看来,题目利用广义线性模型可以解决,在贝叶斯看来同样以解决,但是遇到了一个问题,就是我们得到的非标准后验分布很复杂。我们正是利用MCMC来解决了这个分布的处理问题。本文的重点也在于此。

在使用MCMC时作者遵循了这样的思路,首先依照贝叶斯解决问题的套路,构建了非标准后验分布函数。然后初步运行MCMC,确定合适的scale。继而,确定适当的模拟批次和每批长度(以克服模拟取样的相关性)。最后,估计参数并利用delta方法估计标准误。

1. 问题的提出

这是一个关于R软件中mcmc包的应用案例。问题出自明尼苏达大学统计系博士入学考试试题。这个问题所需要的数据存放在logit数据集中。在这个数据集中有五个变量,其中四个自变量x1、x2、x3、x4,一个响应变量y

对于这个问题,频率学派的处理方法是利用广义线性模型进行参数估计,下面是相应的R代码以及结果:

library(mcmc)
data(logit)
out <- glm(y ~ x1 + x2 + x3 + x4, data = logit, family = binomial(), x = T)
summary(out)

Call:
glm(formula = y ~ x1 + x2 + x3 + x4, family = binomial(), data = logit, 
    x = T)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-1.746  -0.691   0.154   0.704   2.194  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)    0.633      0.301    2.10   0.0354 * 
x1             0.739      0.362    2.04   0.0410 * 
x2             1.114      0.363    3.07   0.0021 **
x3             0.478      0.354    1.35   0.1766   
x4             0.694      0.399    1.74   0.0817 . 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 137.628  on 99  degrees of freedom
Residual deviance:  87.668  on 95  degrees of freedom
AIC: 97.67

Number of Fisher Scoring iterations: 6

但是,使用频率学派的分析方法解决这个问题并不是我们想要的,我们希望使用Bayesian分析方法。对于Bayesian分析而言,我们假定数据模型(即广义线性模型)与频率学派一致。同时假定,五个参数(回归系数)相互独立,并服从均值为0、标准差为2的先验正态分布。

定义下面的R函数来计算非标准的对数后验分布概率密度(先验密度与似然函数相乘)。我们为什么要定义成密度函数的对数形式?因为虽然我们是从分布中采样,但是MCMC算法的执行函数metrop()需要的一个参数正是这个分布的密度函数的对数形式。

x <- out$x
y <- out$y
lupost <- function(beta, x, y) {
    eta <- as.numeric(x %*% beta)
    logp <- ifelse(eta < 0, eta - log1p(exp(eta)), -log1p(exp(-eta)))
    logq <- ifelse(eta < 0, -log1p(exp(eta)), -eta - log1p(exp(-eta)))
    logl <- sum(logp[y == 1]) + sum(logq[y == 0])
    return(logl - sum(beta^2)/8)
}

继续阅读MCMC案例学习

需要相亲几次才能找到靠谱的对象?

谈到相亲就不得不提到著名的麦穗问题。说有一天,苏格拉底带领几个弟子来到一块成熟的麦地边。他对弟子们说:“你们去麦地里摘一个最大的麦穗,但要求只能摘一次,只许进不许退,我在麦地的尽头等你们。”可以看得出,相亲这种活动就有点类似于摘麦穗,在等待和决断之间达成平衡是解决问题的重点。

将上述的麦穗问题进一步抽象就是一个经典的概率问题。若一个袋子里有100个不同的球。每个球上标明了其尺寸大小。我们每次随机无放回的从袋中取一个球出来,观察其大小属性之后需决定要或是不要。如果要,取球就此停止。如果不要, 再继续取球,但不准再回头要原先的球。这样下去,直到100个球取完为止。目的就是取到那个最大的球。

对于这个概率问题,一种思路就是取1到100之间的某个数字n,以它作为分割点将整袋球划分为两组,第一组即从第1个到第n个球,第二组即从第n+1个到第100个球。我们以第一组为观察对象,找到第一组中最大的球M,记录其大小但并不行动。然后从第二组中寻找大于M的第一个球,取该球为最终选择。那么n应该设为多少,才能使取到最大球的概率尽可能的高呢?Ross在其《概率论基础教程》中已经给出了精确的解析(英文版第8版P345;或者陈木法的《随机过程导论》 P105),最优的n公式表达为 $n^*=\inf\{r:\sum_{n=r}^{N-1}\frac{1}{n}\le1\}$;在N充分大时,n应该取在1/e的比例处,也就是所有球数目的37%处。在解的过程中先运用条件概率得到全概率表达式,再用连续化的积分来近似离散化的概率,将积分求出后进行求导得到最终答案。看到这里各位是否有些惊讶,这个e可是无所不在啊。

继续阅读需要相亲几次才能找到靠谱的对象?

有边界区间上的核密度估计

一、一个例子

核密度估计应该是大家常用的一种非参数密度估计方法,从某种程度上来说它的性质比直方图更好,可以替代直方图来展示数据的密度分布。但是相信大家会经常遇到一个问题,那就是有些数据是严格大于或等于零的,在这种情况下,零附近的密度估计往往会出现不理想的情况。下面以一个指数分布的模拟数据为例(样本量为1000),R程序代码为:

set.seed(123);
x=rexp(1000,1);
plot(density(x,kernel="epanechnikov"),ylim=c(0,1.2));
lines(seq(0,8,by=0.02),dexp(seq(0,8,by=0.02),1),col="blue");
abline(v=0,col="red");

可以看出,理论上应该单调递减的密度函数在0附近有明显的“陡坡”,而且不应该有密度的小于零的区域也有着正的估计值。当样本量增大时,这种现象也不会得到明显好转,下图是将样本量改为10000时的情形。

set.seed(123);
x=rexp(10000,1);
plot(density(x,kernel="epanechnikov"),ylim=c(0,1.2));
lines(seq(0,8,by=0.02),dexp(seq(0,8,by=0.02),1),col="blue");
abline(v=0,col="red");

我们也许从平时看的书中了解到,当样本量趋于无穷时,核密度估计值将是收敛到真实的密度函数的,但我们可能不会特意去研究那些结论成立的条件。以上这两个简单的例子似乎给了我们一个直观的感觉,那就是当真实密度函数的支撑集(函数f(x)的支撑集指的是使得f(x)≠0的x的集合)有边界时,核密度估计值可能会出现一些不理想的情况。下面就简单地给出一些理论的结果。

二、理论分析

在一些必要的条件下(真实的密度函数f二阶导绝对连续,三阶导平方可积),核密度估计值$\hat{f}(x)$的偏差有表达式$Bias[\hat{f}(x)]=\frac{h^2\sigma_k^2f”(x)}{2}+O(h^4)$,其中h是带宽,$\sigma_k^2=\int u^2k(u)du$,k(u)是支集为[-1,1]的核函数(即在[-1,1]上有值,其余的地方取零)。可以看出这个偏差是随着带宽h的减小以$h^2$的速度趋于零的。

而假设密度函数以0为边界,那么上述表达式将不再成立,而是代之以
$E[\hat{f}_k(x)]=a_0(x)f(x)-ha_1(x)f'(x)+O(h^2)$
其中$a_i(x)=\int_{-1}^{x/h}u^ik(u)du$。可以看出,当$x \ge h$时,$a_0(x)=1$,$a_1(x)=0$,此时的偏差跟之前的那个表达式没有区别;但当$0 \le x<h$时,$a_0(x)$和$a_1(x)$都是非零的,于是偏差总是存在。

也许你会提议说,将估计值除以$a_0(x)$,偏差就可以减小了吧?的确,这样是一种改进的办法,但也要注意到,此时h的一次项不会消除,也就是说原来$h^2$的衰减速度放慢到了h,从效率上说相对于理想的情况是大打了折扣。

这时候一个巧妙的办法是,用另外一个核函数l(x)对f也做一次估计,那么就有
$E[\hat{f}_l(x)]=b_0(x)f(x)-hb_1(x)f'(x)+O(h^2)$
其中的$b_0$和$b_1$意义类似,只不过是针对l(x)的。

对以上两个式子进行线性组合,则会有
$b_1(x)*E[\hat{f}_k(x)]-a_1(x)*E[\hat{f}_l(x)]=[b_1(x)a_0(x)-a_1(x)b_0(x)]f(x)+O(h^2)$
然后把f(x)的系数移到等式左边,O(h)项的偏差就神奇地消失了。

通过观察核密度估计的表达式,我们可以将上面这个过程等价地认为是对f(x)用了一个新的核函数进行估计,这个新的核函数是
$p(x)=\frac{b_1(x)k(x)-a_1(x)l(x)}{b_1(x)a_0(x)-a_1(x)b_0(x)}$

特别地,如果将l(x)取为x*k(x),那么p(x)将有一个简单的形式
$p(x)=\frac{(a_2(x)-a_1(x)x)k(x)}{a_0(x)a_2(x)-a_1^2(x)}$

当$x \ge h$时,这个新的核函数p(x)就是k(x),而当$x \ge h$时(也就是在边界),它会对最初的核函数进行调整。当$x<0$时,不管算出来的估计值是多少,都只需将密度的估计值取为0即可。

三、程序实现

下面这段程序是对本文的第一幅图进行“整容”,代码及效果图如下:

k=function(x) 3/4*(1-x^2)*(abs(x)<=1);
a0=function(u,h)
{
	lb=-1;
	ub=pmin(u/h,1);
	0.75*(ub-lb)-0.25*(ub^3-lb^3);
}
a1=function(u,h)
{
	lb=-1;
	ub=pmin(u/h,1);
	3/8*(ub^2-lb^2)-3/16*(ub^4-lb^4);
}
a2=function(u,h)
{
	lb=-1;
	ub=pmin(u/h,1);
	0.25*(ub^3-lb^3)-0.15*(ub^5-lb^5);
}
kernel.new=function(x,u,h)
{
	k(x)*(a2(u,h)-a1(u,h)*x)/(a0(u,h)*a2(u,h)-a1(u,h)^2);
}
den.est=function(u,ui,h)
{
	sapply(u,function(u) ifelse(u<0,0,mean(kernel.new((u-ui)/h,u,h))/h));
}
set.seed(123);
dat=rexp(1000,1);
x=seq(0,8,by=0.02);
y=den.est(x,dat,2*bw.nrd0(dat));
plot(x,y,type="l",ylim=c(0,1.2),main="Corrected Kernel");
lines(x,dexp(x,1),col="red");

从中可以看出,边界的偏差问题已经得到了很好的改进。

如果真实的密度函数的支集不是[0,+∞],而是某一个闭区间[m,n],那么偏差修正的过程与上面类似,只不过是要将$a_i(x)$定义为$a_i(x)=\int_{(x-n)/h}^{(x-m)/h}u^ik(u)du$。在编程序的时候,也只需把积分的上下限进行相应的调整即可。

四、参考文献

Jeffrey S. Simonoff, 1998. Smoothing Methods in Statistics. Springer-Verlag

相关链接:http://pages.stern.nyu.edu/~jsimonof/SmoothMeth/