# Persi Diaconis (1)

Persi Diaconis（维基）确实如此，他在搞魔术的时候，为了想研究如何防止被其他魔术师骗，买了本 William Feller 的 An Introduction to Probability Theory and Its Applications，但是里面涉及到了微积分等知识，看不懂，那年他18岁。他发誓要回学校学习，以此能够看得懂这本书。24岁重返校园（City College of New York）。他向《科学美国人》投稿介绍他两个有意思的洗牌方法。这样被一个马丁·葛登能的人看重，给他写了推荐信去哈佛大学，当时哈佛的统计学家 Fred Mosteller 正在研究魔术，于是就要了他（http://blog.sciencenet.cn/home.php?mod=space&amp;uid=1557&amp;do=blog&amp;id=418859）。

Persi Diaconis 做了几个很有意思的工作，如洗牌多少次能够洗得比较彻底（我希望在统计之美里面，能够有一篇来单独介绍洗牌问题）等。他还有个绝活，据刘军老师说，他总可以抛硬币时，抛出他想要的那一面。而他每次的学术报告之前，都会表演一番，很多人实际上不是去听他的报告的，而是看看他的绝活。 win your ex girlfriend backfree advice on how to make her want you back How To Get Your Ex Girlfriend Back how to get back with your ex girlfriendhow to get your girlfriend back

# MCMC案例学习

## 1. 问题的提出

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



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)
}


# 议员是如何投票的？

[latexpage]

## 一、议员投票 Nelson’s votes have tended to be more liberal than conservative. He has received high ratings from left-of-center groups such as Americans for Democratic Action, and low ones from right-of-center groups such as the Eagle Forum and the Club for Growth. According to ratings by the National Journal, Nelson’s votes have been liberal on economic matters, moderate on social issues, and liberal but close to the center on foreign policy.

While she(Cantwell) scores high on a progressive chart from ProgressivePunch.org, Cantwell has made several controversial votes during her time in the Senate that have created friction between her and members of the Democratic Party.

Enzi was ranked by National Journal as the sixth-most conservative U.S. Senator in its March 2007 conservative/liberal rankings. Despite his strong support of the War in Iraq, he was one of 14 U.S. Senators to vote against the Iraq War funding bill in May 2007 because he opposes the clauses of the bill which increase domestic spending.

## 二、图模型——罚极大似然

$f(x)=\frac{1}{(2\pi)^{p/2}|\Sigma|^{1/2}}\exp\Big\{-\frac{1}{2}(X-\mu)^T\Sigma^{-1}(X-\mu)\Big\}$

$S=\frac{1}{n}\sum_{i=1}^N(x_i-\bar{x})(x_i-\bar{x})^T$

$\rho_{ij|\\\{i,j\}}=-\frac{\omega_{ij}}{\sqrt{\omega_{ii}\omega_{jj}}}$

$\begin{split} \sum_{i=1}^n\log & f(x_i)=\frac{np}{2}\log(2\pi)-\frac{n}{2}\log|\Sigma| -\frac{1}{2}\sum_{i=1}^n(x_i-\bar{x})^T\Sigma^{-1}(x_i-\bar{x})\\ &=\frac{np}{2}\log(2\pi)-\frac{n}{2}\log|\Sigma| -\frac{1}{2}tr\Big(\sum_{i=1}^n(x_i-\bar{x})(x_i-\bar{x})^T\Sigma^{-1}\Big)\\ &\propto p\log(2\pi)-\log|\Sigma|-tr\big(S\Sigma^{-1}\big) \end{split}$

$tr(S\Omega)-\log|\Omega|$

\begin{equation}
\min_{\Omega\succ 0}\Big\{tr(S\Omega)-\log|\Omega|+\lambda\sum_{i\neq j}|\omega_{ij}|\Big\} \label{obj:main}
\end{equation}

求解（1）式，得到$\Omega$矩阵，它的非零元素便对应“投票图”中的一条边。

## 三、数据与代码

R数据下载——数据中包含三个矩阵（数据来自http://www.senate.gov）。

• idList矩阵为参议员信息。每行对应一个参议员。5列分别对应：议员id，last name, first name, state, party(Democrat, Republican or Independent)
• bilInfo为投票议案信息(bil info)。每行对应一项议案信息，4列分别为：议案编号、投票时间、议案描述、投票结果。
• voteMatrix矩阵为投票信息，每行对应bilInfo的议案，每列对于idList的一个议员。若voteMatrix的第i行第j列数值为1，表示议员j对议案i投了赞成票；-1表示反对票；NA表示未投票。

library(spaceExt)
library(igraph)
#移除投票缺失较多的议员
sel=which(!(colSums(is.na(voteMatrix))>100))
partyD=as.numeric(idList[sel,5]=="D")
partyI=as.numeric(idList[sel,5]=="ID" | idList[sel,5]=="I")
#用spaceExt做计算
res=glasso.miss(voteMatrix[,sel],20,rho=0.1,penalize.diagonal=FALSE)
#计算偏相关系数矩阵
p=-res\$wi d=1/sqrt(diag(res\$wi))
p=diag(d)%*%p%*%diag(d)
diag(p)=0
#igraph包生成图模型、作图
V(g)\$color=(partyD+2+partyI*2) V(g)\$label=idList[sel,3]
par(mar=c(0,0,0,0))
plot(g,layout=layout.fruchterman.reingold, edge.color=grey(0.5),vertex.size=10)

1. Yuan, M., and Lin, Y. (2007) Model selection and estimation in the Gaussian graphical model. Biometrika , 94, 1, pp. 19–35
2. Friedman, J. H., Hastie T, Tibshirani R. “Sparse inverse covariance estimation with the graphical lasso.” Biostat (2008) 9 (3): 432-441.
3. Banerjee, O., Ghaoui, L. E. and d’Aspremont, A. (2007), Model selection through sparse maximum likelihood Estimation, J. Machine Learning Research 101.
4. Wang, H. The Bayesian Graphical Lasso
5. Mladen Kolar, Le Song, Amr Ahmed, Eric P. Xing. Estimating time-varying networks. Annals of Applied Stat. arXiv:0812.5087v2

# Sweave后传：统计报告中的大规模计算与缓存

Sweave本身做不到缓存（或者很难做到），但由于R语言的扩展性以及一些疯子和天才的存在，我们可以通过R的一些附加包来完成缓存。这里的缓存有两种（不关心细节的读者可以略过）：

• 一种是对R对象的缓存，比如数据、模型（拟合好的模型）和其它计算结果等，这个缓存容易理解，我们可以把它简单想象为把计算过程中的对象存入某些特殊的缓存库（文件），下次计算的时候先查看缓存库中是否存在我们需要的对象，如果存在，那么就直接加载进来，否则就重新计算；如（伪代码）：
## 如果存在，就加载
if (file.exists('x.RData')) {
}
## 如果当前工作环境中不存在x，那么就重新创建它
if (!exists('x')) {
x = rnorm(100000)
}

R对象的缓存可以通过cacheSweave包（作者Roger D. Peng）实现，它真正采用的方法比上面的伪代码当然要高级（比如用代码的MD5值作数据库名等），不过用户可以不必关心这其中的细节。这种“当需要时才加载”的方式在术语上叫“延迟加载”，即Lazy Load，这在R里面也很常见，比如很多R包的数据都是采用延迟加载的方式（加载包的时候并不立刻加载数据，而是用data(name)在需要的时候加载）。

• 另一种是图形的缓存，这是一个不可思议的魔法。cacheSweave的帮助文档和vignette中特别提到了它只能缓存R对象，无法缓存图形和其它附属输出，而开源世界的魅力就是你总能找到一些奇妙的解决方案。R包pgfSweave借力于pgf实现了图形的缓存。pgfSweave包的两位作者Cameron Bracken和Charlie Sharpsteen都是水文学相关专业出身，却在R的世界里做了这些奇妙的工作，让我觉得有点意外，这是题外话。pgf是LaTeX世界的又一项重大发明，它提供了一套全新的绘图方式（语言），而且它的图形可以通过LaTeX直接生成PDF输出（这是pgfSweave能实现缓存的关键）。关于pgf，我们可以翻一翻它726页的手册，如果你能坚持看上10分钟，感叹超过50次，那么你一定是一位超级排版爱好者。它的图形质量之高、输出之漂亮，真的是让人（至少让我）叹为观止。同样，这里我也不详细介绍细节。pgf图形的“缓存”是通过“外置化”（externalization）来实现的，简言之，pgf图形有一种输出形式如下：
\beginpgfgraphicnamed{graph-output-1}
\input{graph-output-1.tikz}
\endpgfgraphicnamed


如果LaTeX在运行过程中遇到这样的代码段，那么它会把tikz图形（tikz是高层实现，pgf是基础设施）先编译为PDF图形，等下次重复运行LaTeX的时候，这段代码就不必重新运行了，LaTeX会跳过它直接引入PDF图形。这样就省去了编译tikz图形的时间，同时也能得到高质量的图形输出。

pgfSweave从R的层面上可以输出这样的代码段，并生成一个shell脚本来编译图形为PDF；而pgfSweave也会尝试根据R代码的MD5来判断一段画图的代码是否需要重复运行——如果R代码没有变化，那么图形就不必重制。这是自然而然的“缓存”。对了，pgfSweave对Sweave的图形输出作了扩展，新增了tikz输出，这是通过另一个R包tikzDevice实现的（同样两个作者），这个包可以将R图形录制为tikz绘图代码。pgfSweave依赖于cacheSweave包，把它缓存R对象的机制也引入进来，所以我们可以不必管cacheSweave，把精力放在pgfSweave上就可以了。

Gibbs抽样的最终目的是从一个多维分布$f(X_1,X_2,\ldots,X_n)$中生成随机数，而已知的是所有的条件分布$f(X_1|X_2,\ldots,X_n)$、$f(X_2|X_1,\ldots,X_n)$、……、$f(X_n|X_1,\ldots,X_{n-1})$，并且假设我们可以很方便从这些条件分布中抽样。剩下的事情就很简单了：任意给定初始值$x_1^{(0)},x_2^{(0)},\ldots,x_n^{(0)}$，括号中的数字表示迭代次数，接下来依次迭代，随机生成$x_1^{(k)} \sim f(X_1|x_2^{(k-1)},\ldots,X_n^{(k-1)})$、$x_2^{(k)} \sim f(X_2|x_1^{(k)},\ldots,x_n^{(k-1)})$、……、$x_n^{(k)} \sim f(X_n|x_1^{(k)},\ldots,x_{n-1}^{(k)})$（迭代次数$k=1,\cdots,N$），也就是，根据上一次的迭代结果并已知条件分布的情况下，从条件分布中生成下一次的随机数（这个算法很简单，你什么时候看明白上标下标了就明白了）。这是个串行的过程，所以无法避免循环。

$\left[\begin{array}{c}X\\ Y\end{array}\right]\sim\mathcal{N}\left(\left[\begin{array}{c}0\\ 1\end{array}\right],\left[\begin{array}{cc} 4 & 4.2\\ 4.2 & 9\end{array}\right]\right)$

A Simple Demo on Caching R Objects and Graphics with pgfSweave (PDF)

A Simple Demo on Caching R Objects and Graphics with pgfSweave (LyX)

• cache=TRUE 使得文档中的R对象可以被缓存；
• external=TRUE 使得图形可以被缓存；

1. pgfSweave对图形的缓存经常有点过度，即使你改动了代码，缓存仍然不会被替换，这种情况下你可以把生成的tikz和pdf图形删掉；
2. 中文问题：tikzDevice包好像还没有完全解决中文问题，作者正在尝试我提供的中文图形，也许在未来不久的新版本中会有完全的中文支持，但目前实际上我已经在我的自动配置脚本中解决了这个问题。我们只需要保证我们的中文LyX文档用UTF-8编码就可以了——在文档选项中使用CTeX类（以及UTF8选项），并设置语言为中文、编码为Unicode (XeLaTeX) (utf8)。这种设置至少我在Windows和MikTeX下可以成功编译中文Sweave文档。目前高亮代码功能在中文文档中还不能用，但我已经向相关作者发了补丁，等CRAN管理员回来发布新版parser包之后这个问题也会解决。所以，中文用户也不必担心，一切都（将）可以使用。
3. 依旧是样式问题：细心的读者也许能发现上面的PDF文档中的图形和这里网页中的图形略有区别，PDF文档中图形的字体和正文字体是一致的，这也是pgf的优势，前文第5节曾提到过。
4. 如果使用缓存，那么代码段中的对象都无法打印出来，即使用print(dat)也无法打印dat，这是因为缓存的代码段只管是否加载对象，而无法执行进一步的操作（包括打印），对于这种情况，我们应该设计好我们的Sweave文档结构，让某些代码段作计算（并缓存），某些代码段作输出（不要缓存）。

# 从中心极限定理的模拟到正态分布

## 二、中心极限定理的模拟

### 1、对称的钟形曲线不代表正态分布

data a;
do x1=1 to 6;
do x2=1 to 6;
output;
end;
end;  /*模拟掷骰子两次，生成36行数据*/
data a;
set a;
x=sum(x1,x2);
proc univariate data=a noprint;
var x;
histogram/normal (mu=est sigma=est);
run;

x = rowSums(expand.grid(1:6, 1:6))
hist(x, breaks = seq(min(x) - 0.5, max(x) + 0.5, 1), main = "")


d = c(1, 2, 3, 4, 5, 15)
x = rowSums(expand.grid(d, d))
hist(x, breaks = seq(min(x) - 0.5, max(x) + 0.5, 1), main = "")


### 3、给出拟合好坏的指标

if (!require(animation)) install.packages("animation")
library(animation)
ani.options(interval = 0.1, nmax = 100)
par(mar = c(4, 4, 1, 0.5))
clt.ani()

## 脚注

symbol;
goptions ftext= ctext= htext=;

%macro a(n);
data rv;
retain _seed_ 0;
do _i_ = 1 to &n;
binom1 = ranbin(_seed_,&n,0.8); /*ranbin函数返回以seed、n、0.8为参数的二项分布随机变量*/
output;
end;
drop _seed_ _i_;
run;
%mend;

%a(100);

proc univariate data=rv noprint;
var binom1;
histogram/normal( mu=est sigma=est);
inset normal ;
run;

R代码：

a = function(n) hist(rbinom(n, n, 0.8))