% $Header: /cvsroot/latex-beamer/latex-beamer/examples/beamerexample4.tex,v 1.4 2004/10/07 20:53:07 tantau Exp $

\documentclass[envcountsect,cjk]{beamer}
%\setbeamertemplate{theorems}[numbered]
% Copyright 2003 by Till Tantau <tantau@users.sourceforge.net>.
%
% This program can be redistributed and/or modified under the terms
% of the LaTeX Project Public License Distributed from CTAN
% archives in directory macros/latex/base/lppl.txt.
%
% Many thanks to Huang Yushuo for helping me with this file
% (I don't speak Chinese...)
%\setbeamercovered{dynamic}
\setbeamercovered{transparent}
%\usepackage{times}
%\usepackage{mathtime}
\usepackage{graphicx}
\usepackage{Picins}
\usepackage{indentfirst}
\usepackage{amsmath,amssymb,amsfonts,amsthm}
\usepackage{multirow}
\usepackage{textcomp}

\usepackage{mathrsfs}
\usepackage{latexsym}
\usepackage{makeidx}% 制作索引
\usepackage{accents}
\usepackage{slashbox}
\usepackage{picinpar}
%\makeindex
\usepackage{undertilde}
\usepackage{bm}

%\usepackage{textcomp}
\setlength{\abovedisplayskip}{2pt plus1pt minus1pt}     %公式前的距离
\setlength{\belowdisplayskip}{2pt plus1pt minus1pt}     %公式后面的距离
\setlength{\arraycolsep}{2pt}   %在一个array中列之间的空白长度, 因为原来的太宽了
\allowdisplaybreaks[4]  % \eqnarray如果很长, 影响分栏、换行和分页
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\renewcommand{\theequation}{\thesection.\arabic{equation}}
\newtheorem{remark}{Remark}
\makeatletter\@addtoreset{equation}{section}\makeatother
%\makeatletter\@addtoreset{section}{part}\makeatother
\newcommand{\D}{\scriptscriptstyle}
\setlength{\parindent}{20pt}                %Paragraph indent amounts.
\setlength{\parskip}{3pt plus1pt minus2pt}  %Vertical distance between par
\setlength{\baselineskip}{20pt plus2pt minus1pt}%Distance between lines in a par.
%\usepackage{graphicx}
\makeatletter
\def\widebar{\accentset{{\cc@style\underline{\mskip10mu}}}}
\def\Widebar{\accentset{{\cc@style\underline{\mskip8mu}}}}
\makeatother
\def\dl{\displaystyle}

\usepackage{CJK}


%\usetheme{Warsaw}
%\usetheme{Hannover}
%\usetheme{Darmstadt}
\usetheme{JuanLesPins}
\AtBeginSection[]{\frame{\tableofcontents[current]}}
\AtBeginPart{\frame{\partpage}}
\begin{document}
  \begin{CJK*}{GBK}{rm}

\renewcommand{\raggedright}{\leftskip=0pt \rightskip=0pt plus 0cm}

    \title[贝叶斯分析]{\bf 贝~~~~叶~~~~斯~~~~分~~~~析}
    \author{姓名：韦来生~~~张伟平}
\date[\today]{中~~国~~科~~大~~出~~版~~社 \\ \vskip 1cm
     第一版, ~2013}

    \frame{\titlepage}


\begin{frame}
  \titlepage
\end{frame}

\begin{frame}{第5章 ~贝叶斯计算方法}
\scriptsize
\tableofcontents
\end{frame}

\section{5.1 ~引言}

\frame{\frametitle{引言}

\begin{itemize}

\item
在Bayes分析中常常需要计算后验分布的期望、方差、分位数或众数等数字特征.
如果先验分布不是共轭先验，那么后验分布往往不再是标准的分布.
因此需要计算的后验分布数字特征往往没有显式表达，这就需要一些特殊的计算方法.
\medskip\pause

\item {\bf 例~5.1.1} ~假设随机变量~$X\sim N(\theta,\sigma^2)$，其中~$\sigma^2$
已知. 若选取~$\theta$ 的先验分布为柯西分布 ~$C(\mu,\tau),$
则~$\theta$ 的后验分布为\vspace{-2mm}
$$
 \pi(\theta|x)\propto exp\{-(\theta-x)^2/(2\sigma^2)\}(\tau^2+(\theta-\mu)^2)^{-1}.\vspace{-2mm}
$$
因此后验期望和后验方差为\vspace{-2mm}
$$
E^\pi(\theta|x)=\frac{\dl\int_{-\infty}^\infty \theta
exp\{-(\theta-x)^2/(2\sigma^2)\}(\tau^2+(\theta-\mu)^2)^{-1}d\theta}
{\dl\int_{-\infty}^\infty
exp\{-(\theta-x)^2/(2\sigma^2)\}(\tau^2+(\theta-\mu)^2)^{-1}d\theta},
$$
\end{itemize}
}

\frame{\frametitle{引言}

\vspace{-7mm}

\begin{eqnarray*}
V^\pi(\theta|x)\!\!\!&&=\frac{\dl\int_{-\infty}^\infty \theta^2
exp\{-(\theta-x)^2/(2\sigma^2)\}(\tau^2+(\theta-\mu)^2)^{-1}d\theta}
{\dl\int_{-\infty}^\infty
exp\{-(\theta-x)^2/(2\sigma^2)\}(\tau^2+(\theta-\mu)^2)^{-1}d\theta}\\
&&~~-(E^\pi(\theta|x))^2.
\end{eqnarray*}\\\vspace{-3mm}
\begin{itemize}\pause

\item 显然上面两个积分没有显式解，各种数值积分方法
可被用来逼近这两个积分. 下面的例子给出了一个更困难的问题.
\medskip\pause

\item {\bf 例~5.1.2}
~假设~$X_1,\ldots,X_k$ 为独立的普阿松~(poisson) 随机变量，且
~$X_i\sim Poisson(\theta_i), i=1,\ldots,k.$ 如果~$\theta_i$
的先验分布为通过其对数的联合分布给出:\vspace{-2mm}
$$
\bm\nu=(\log(\theta_1),\ldots,\log(\theta_k))'\sim N(\mu
\mathbf{1}_k,
 \tau^2\{(1-\rho)\bm{I}_k +\rho \bm{J}_k \}),\vspace{-2mm}
$$
其中~$\mathbf{1}_k$ 为元素是~1 的~$k$ 维列向量，$\bm{I}_k$ 为~$k$
阶单位阵，$\bm{J}_k$ 的~$k$ 阶元素全为~1 的方阵，$\mu, \tau^2, \rho$
为已知的常数.
\end{itemize}
}

\frame{\frametitle{引言} %\vspace{-2mm}

\begin{itemize}

\item 则由样本分布\vspace{-3mm}
$$
f(\mathbf{x}|\bm\nu)=exp\Big\{-\sum_{i=1}^k(e^{\nu_i}-\nu_ix_i)\Big\}\Big/\prod_{i=1}^kx_i!\vspace{-3mm}
$$
以及先验分布\vspace{-2mm}
$$ \pi(\bm\nu)\propto
exp\Big\{-\frac{1}{2\tau^2}(\bm\nu-\mu\mathbf{1}_k)'[(1-\rho)\bm{I}_k
+\rho \bm{J}_k ]^{-1}(\bm\nu-\mu\mathbf{1}_k)\Big\}\vspace{-4mm}
$$\pause

\item
可以得到后验分布\vspace{-3mm}
\begin{eqnarray*}
\pi(\bm\nu|\mathbf{x})\!\!\!\!&&\propto g(\bm\nu|\mathbf{x})
=exp\bigg\{-\sum_{i=1}^k(e^{\nu_i}-\nu_ix_i)\\
&&~~-\frac{1}{2\tau^2}(\nu-\mu\mathbf{1}_k)'[(1-\rho)\bm{I}_k +\rho
\bm{J}_k ]^{-1}(\nu-\mu\mathbf{1}_k)\bigg\}.\vspace{-2mm}
\end{eqnarray*}
\end{itemize}
}

\frame{\frametitle{引言}

\begin{itemize}

\item 因此如果感兴趣的是~$\theta_j$ 的后验均值，则需要计算\vspace{-2mm}
\[
E^\pi(\theta_j|\mathbf{x})=E^\pi(e^{\nu_j}|\mathbf{x})
=\frac{\dl\int_{\mathbb{R}^k}
e^{\nu_j}g(\bm\nu|\mathbf{x})d\nu}{\dl\int_{\mathbb{R}^k}
g(\bm\nu|\mathbf{x})d\nu},\vspace{-2mm}
\]
\pause

\item 这是两个~$k$ 重积分的比值. 当~$k$ 越大时就越难处理，而数值积分方法在这种场合下不再是有效的方法，
这种问题也就是常称的
维数灾难问题. %其是因为与积分计算无关的空间部分大小随着维数的增加快速增加，因而使得
数值逼近中的误差随着维数~$k$ 的幂次增加，最终导致算法失效.
因此数值积分方法在一维和二维积分以外的场合下不是优先使用的方法.
\medskip\pause

\item 近年来贝叶斯统计方法的快速普及发展，其中
包括~E-M 方法、马尔科夫链蒙特卡洛~(MCMC) 抽样方法中的~Gibbs
抽样方法和~Metropolis-Hastings 算法等等.
\end{itemize}
}

\section{5.2 ~分析逼近方法}

\frame{\frametitle{Laplace 逼近}
\begin{itemize}

\item 在计算一些积分时，
可以使用~Laplace 逼近方法来逼近所要计算的积分.
具体的，假设\vspace{-3mm}
$$\qquad\quad
E^\pi(g(\bm\theta)|x)=\frac{\dl\int_{\mathbb{R}^k}g(\bm\theta)f(x|\bm\theta)
\pi(\bm\theta)d\bm\theta}{\dl\int_{\mathbb{R}^k}f(x|\bm\theta)\pi(\bm\theta)d\bm\theta},
\qquad\quad(5.2.1)\vspace{-3mm}
$$
为感兴趣的量，其中~$g,~f,~\pi$
均为~$\theta$~的光滑可微函数.\medskip\pause

\item 将上述分子和分母上的积分统一考虑为如下形式的积分\vspace{-2mm}
$$
I=\int_{\mathbb{R}^k}q(\theta)exp\{-nh(\theta)\}d\theta,
\vspace{-2mm}
$$
其中~$-nh(\bm\theta)=ln[f(x|\bm\theta)\pi(\bm\theta)],$ 而~$h$
为一光滑函数且~$-h$ 有唯一的最大值点~$\hat\theta.$ 则~Laplace
逼近方法需要对~$q$ 和~$h$ 在~$\hat\theta$ 处进行~Talyor 展开.
\end{itemize}
}

\frame{\frametitle{$k=1$ 情形时的~Laplace 逼近}
\begin{itemize}

\item 我们首先考虑$k=1$ 的情形. 积分为\vspace{-1mm}
$$I=\int_{-\infty}^{\infty}q(\theta)exp\{-nh(\theta)\}d\theta,\vspace{-2mm}
$$\pause

\item 令~$\hat\theta$ 为~$-h(\theta)$ 的唯一极大似然估计，可将积分视为
~$\hat\theta\pm\delta$ 上的积分，即\vspace{-2mm}
$$I\approx\int_{\hat\theta-\delta}^{\hat\theta+\delta}q(\theta)exp\{-nh(\theta)\}d\theta.\vspace{-2mm}
$$\pause

\item 对~$q$ 和~$h$ 在~$\hat\theta$ 处进行~Taylor 展开\vspace{-2mm}
\begin{eqnarray*}
q(\theta)\!\!\!\!&&=q(\hat\theta)+(\theta-\hat\theta)q'(\hat\theta)
+\frac{1}{2}(\theta-\hat\theta)^2q''(\hat\theta)+\mbox{余项}\\
h(\theta)\!\!\!\!&&=h(\hat\theta)+(\theta-\hat\theta)h'(\hat\theta)
+\frac{1}{2}(\theta-\hat\theta)^2h''(\hat\theta)+\mbox{余项}\\
\end{eqnarray*}
\end{itemize}
}

\frame{\frametitle{$k=1$ 情形时的~Laplace 逼近}

\noindent 注意到~\!$\hat\theta$\! 是~\!$-h(\theta)$\!
的极值点，故上述~\!$h(\hat\theta)$\! 展开项中
~\!$(\theta-\hat\theta)h'(\hat\theta)$\! $=0,$ 因此\vspace{-2mm}
\begin{eqnarray*}
&&q(\theta)\exp\{-nh(\theta)\}\!=\!\Big[q(\hat\theta)\!+\!(\theta\!-\!\hat\theta)q'(\hat\theta)
\!+\!\frac{1}{2}(\theta\!-\!\hat\theta)^2q''(\hat\theta)\!+\!\cdots\Big]\\
&&~~~~~~~\times
\exp\Big\{-n\Big[h(\hat\theta)+\frac{1}{2}(\theta-\hat\theta)^2h''(\hat\theta)+\cdots\Big]\Big\}\\
&&\approx
e^{-nh(\hat\theta)}q(\hat\theta)\Big[1+(\theta-\hat\theta)\frac{q'(\hat\theta)}{q(\hat\theta)}
+\frac{1}{2}(\theta-\hat\theta)^2\frac{q''(\hat\theta)}{q(\hat\theta)}\Big]\\
&&~~~~~~~\times\exp\Big\{-\frac{n}{2}[(\theta-\hat\theta)^2h''(\hat\theta)]\Big\}
\end{eqnarray*}
}

\frame{\frametitle{$k=1$ 情形时的~Laplace 逼近}
%\begin{itemize}

\noindent 将上式代入到积分式中，得到
\begin{eqnarray*}
I\!\!\!\!&&\approx\int_{\hat\theta-\delta}^{\hat\theta+\delta}\Big[q(\hat\theta)+(\theta-\hat\theta)q'(\hat\theta)
+\frac{1}{2}(\theta-\hat\theta)^2q''(\hat\theta)\Big]\\
&&~~~~\times
\exp\Big\{-n\Big[h(\hat\theta)+\frac{1}{2}(\theta-\hat\theta)^2h''(\hat\theta)\Big]\Big\}d\theta\\
&&=e^{-nh(\hat\theta)}q(\hat\theta)\int_{\hat\theta-\delta}^{\hat\theta+\delta}
\Big[1+(\theta-\hat\theta)\frac{q'(\hat\theta)}{q(\hat\theta)}
+\frac{1}{2}(\theta-\hat\theta)^2\frac{q''(\hat\theta)}{q(\hat\theta)}\Big]\\
&&~~~~\times\exp\Big\{-\frac{n}{2}[(\theta-\hat\theta)^2h''(\hat\theta)]\Big\}d\theta
\end{eqnarray*}

%\end{itemize}
}

\frame{\frametitle{$k=1$ 情形时的~Laplace 逼近}
\begin{itemize}

\item 令~$c=h''(\theta),~t=\sqrt{nc}(\theta-\hat\theta),~d\theta=\frac{dt}{\sqrt{nc}},$
故有\vspace{-2mm}
\begin{eqnarray*}
I\!\!\!\!&&\approx
e^{-nh(\hat\theta)}q(\hat\theta)\sqrt{\frac{2\pi}{nc}}
\int_{-\delta\sqrt{nc}}^{\delta\sqrt{nc}}\Big[1+\frac{t}{\sqrt{nc}}\frac{q'(\hat\theta)}{q(\hat\theta)}
+\frac{t^2}{2nc}\cdot\frac{q''(\hat\theta}{q(\hat\theta)} \Big]\\
&&~~~\times\frac{1}{\sqrt{2\pi}}e^{-\frac{t^2}{2}}dt.
\end{eqnarray*}\\\vspace{-2mm}\pause

\item 由于$\dl\int_{-\delta\sqrt{nc}}^{\delta\sqrt{nc}} t
e^{-\frac{t^2}{2}}dt=0,~\dl\int_{-\delta\sqrt{nc}}^{\delta\sqrt{nc}}t^2
e^{-\frac{t^2}{2}}dt=1$, 故有\vspace{-2mm}
\begin{eqnarray*}
I\!\!\!\!&&\approx
e^{-nh(\hat\theta)}q(\hat\theta)\sqrt{\frac{2\pi}{nc}}
\Big[1+\frac{1}{2nc}\cdot\frac{q''(\hat\theta)}{q(\hat\theta)}\Big]\\
&&=e^{-nh(\hat\theta)}q(\hat\theta)\sqrt{\frac{2\pi}{nc}}\big[1+O(n^{-1})\big].
\end{eqnarray*}
\end{itemize}
}

\frame{\frametitle{一般情形的~Laplace 逼近}
\begin{itemize}

\item {\bf 例~5.2.1} ~用~Laplace 逼近方法计算积分\vspace{-2mm}
$$n!=\int_0^\infty x^n
e^{-x}dx=\int_0^\infty x^n e^{-n(\frac{x}{n}-\ln x)}dx.\vspace{-3mm}
$$\pause

\item 将上述~\!$k=1$\! 时的结论推广到~$\theta$ 的维数~\!$k>1$\! 的情形, 此时\vspace{-1mm}
$$
I=\int_{\mathbb{R}^k}q(\bm\theta)exp\{-nh(\bm\theta)\}d\bm\theta.\vspace{-2mm}
\vspace{-1mm}
$$
记~$\hat{\bm\theta}$ 为~$-h(\bm\theta)$ 的~MLE,
~$\bigtriangleup_k(\bm\theta)=\Big(\frac{-\partial^2 h(\bm\theta)}
{\partial\theta_i\partial\theta_j}\Big)$ 为~Hassian 阵. 用于~$k=1$
类似的方法对~$q(\bm\theta)$ 和~$h(\bm\theta)$ 作~Taylor
展开，代入到积分中得到\vspace{-1mm}
$$I\approx
e^{-nh(\hat{\bm\theta})}q(\hat{\bm\theta})(2\pi)^{-k/2}n^{-k/2}|\bigtriangleup_k(\hat{\bm\theta})|^{-1/2}
\big[1+O(n^{-1})\big].\vspace{-4mm}
$$\pause

\item 其中~$\bigtriangleup_k(\bm\theta)$
为由所有二阶偏导数构成的~$k\times k$
维的方阵，$|\bigtriangleup_k(\bm\theta)|$
为此方阵的行列式.%\medskip
\end{itemize}
}

\frame{\frametitle{一般情形的~Laplace 逼近}
\begin{itemize}

\item 因此积分~(5.2.1) 为\vspace{-2mm}
\begin{eqnarray*}
&&E^\pi(g(\bm\theta)|x)=\frac{\int_{\mathbb{R}^k}g(\bm\theta)f(x|\bm\theta)
\pi(\bm\theta)d\bm\theta}{\int_{\mathbb{R}^k}f(x|\bm\theta)\pi(\bm\theta)d\bm\theta}\\
&&=\frac{\int_{\mathbb{R}^k}g(\bm\theta)\exp\{-nh(\bm\theta)\}d\bm\theta}
{\int_{\mathbb{R}^k} 1\cdot \exp\{-nh(\bm\theta)\}d\bm\theta}\\
&&=\frac{e^{-nh(\hat{\bm\theta})}g(\hat{\bm\theta})(2\pi)^{-k/2}n^{-k/2}
|\bigtriangleup_k(\hat{\bm\theta})|^{-1/2}\big[1+O(n^{-1})\big]} {
e^{-nh(\hat{\bm\theta})}~1\cdot(2\pi)^{-k/2}n^{-k/2}
|\bigtriangleup_k(\hat{\bm\theta})|^{-1/2}\big[1+O(n^{-1})\big]}\\
&&=g(\hat{\bm\theta})\big[1+O(n^{-1})\big] %\qquad\qquad (5.2.4)
\end{eqnarray*}\\\vspace{-2mm}
因此通过~Laplace 逼近获得了近似计算的结果.\medskip\pause

\item 如果在~Taylor
展开式中有更多的项被保留，则二阶逼近也是容易得到的.

\end{itemize}
}


\section{5.3 ~E-M 方法}

\frame{\frametitle{引言}

\begin{itemize}
\item E-M~~\!(Expectation-Maximization)\!~ 算法是由~Dempster et al.
(1977)~提出的一种最优化方法，常用于不完全数据求极大似然估计.
\medskip\pause

\item 设~$Y|\theta$ 有密度~$f(y|\theta),$ 令~$\theta$~的先验分布为~$\pi(\theta)$，
由此得到的后验分布记为~$\pi(\theta|y).$ 当计算~$\pi(\theta|y)$
数字特征非常困难时，有一些“数据扩张”的方法或许可以用来解决此类困难.
\medskip\pause

\item 其想法是将观测到的数据~$y$ 与缺失数据或者隐变量数据~$z$ 扩张为“完全”数据~$\bm x=(y,z),$
使得扩张后的后验分布~$\pi(\theta|x)$ $=\pi(\theta|y,z)$
在计算上是容易处理的. \medskip\pause

\item 在贝叶斯计算中，E-M 算法只能计算后验众数. 但是，如果数据扩张方法
能够导致一个易于计算的后验分布，则可以使用更有效的计算工具来计算后验分布的各种数字特征.
\end{itemize}
}

\frame{\frametitle{E-M 方法的步骤}
\begin{itemize}

\item
记~$p(z|y,\hat\theta)$ 为~$Z$ 在给定~$y,\hat\theta$~
时候的预测分布，$\hat\theta^{(i)}$ 为在第~$i$ 次迭代时候 ~$\theta$
的估计值，则~E-M 算法的基本步骤如下：
\medskip\pause

\item (1) 计算~$z^{(i)}=E[Z|y,\hat\theta^{(i)}];$
\medskip\pause

\item (2) 扩充观测数据~$y$ 为~$(y,z^{(i)}),$ 最大化~$\pi(\theta|y,z^{(i)}),$
记其最大值为~$\hat\theta^{(i+1)};$
\medskip\pause

\item(3)
使用~$\hat\theta^{(i+1)}$ 和~(1)，获得~$z^{(i+1)},$ 然后再代入~(2),
如此重复下去直至达到收敛要求.
\medskip\pause

\item 这种先求期望，然后求最大值的步骤即为~E-M 算法名称的来源.
\end{itemize}
}

\frame{\frametitle{E-M 算法的具体实施}
\begin{itemize}

\item 注意到~$\pi(\theta|y)=\pi(\theta,z|y)/p(z|y,\theta),$
我们有\vspace{-2mm}
$$
\ln\pi(\theta|y)=\ln\pi(\theta,z|y)-\ln p(z|y,\theta).\vspace{-5mm}
$$\pause

\item 两边同时对~$Z|y,\hat\theta^{(i)}$ 取期望，则有\vspace{-2mm}
\begin{eqnarray*}
\ln\pi(\theta|y)\!\!\!&&=\int\ln\pi(\theta,z|y)p(z|y,\hat\theta^{(i)})dz\\
&&~~~-\int\ln p(z|y,\theta)p(z|y,\hat\theta^{(i)})dz\\
&&=Q(\theta,\hat\theta^{(i)})-H(\theta,\hat\theta^{(i)}),\hspace{3cm}
(5.3.1)
\end{eqnarray*}\\\vspace{-1mm}\pause


\end{itemize}
}

\frame{\frametitle{算法的具体实施}
\begin{itemize}

\item 则~E-M 算法在第~$i$ 步需要计算如下两个步骤：\pause\smallskip

\textbf{E 步} ~计算~$Q(\theta,\hat\theta^{(i)});$\pause

\smallskip
\textbf{M 步} ~对~$\theta$ 最大化~$Q(\theta,\hat\theta^{(i)})$
以获得~$\hat\theta^{(i+1)},$ 即\vspace{-2mm}
  $$Q(\hat\theta^{(i+1)},\hat\theta^{(i)})=\max_\theta Q(\theta,\hat\theta^{(i)}).\vspace{-4mm}
  $$\pause

\item
由~(5.3.1)，注意到\vspace{-2mm}
\begin{eqnarray*}
&&\ln\pi(\hat\theta^{(i+1)}|y)-\ln\pi(\hat\theta^{(i)}|y)=\big\{Q(\hat\theta^{(i+1)},\hat\theta^{(i)})
-Q(\hat\theta^{(i)},\hat\theta^{(i)})\big\}\nonumber\\
&&~~~~-\big\{H(\hat\theta^{(i+1)},\hat\theta^{(i)})-H(\hat\theta^{(i)},\hat\theta^{(i)})\big\}.
\hspace{2.6cm}(5.3.2)
\end{eqnarray*}\\\vspace{0mm}\pause

\item 由~M 步知道~$Q(\hat\theta^{(i+1)},\hat\theta^{(i)})\ge
Q(\hat\theta^{(i)},\hat\theta^{(i)}).$
\medskip\pause

\end{itemize}
}

\frame{\frametitle{算法的具体实施}
\begin{itemize}

\item 进一步，对任何~$\theta$, 有\vspace{-2mm}
\begin{eqnarray*}
&&H(\theta,\hat\theta^{(i)})-H(\hat\theta^{(i)},\hat\theta^{(i)})\\
&&=\int\ln p(z|y,\theta)p(z|y,\hat\theta^{(i)})dz-\int\ln p(z|y,\hat\theta^{(i)})p(z|y,\hat\theta^{(i)})dz \\
&&=\int\ln \frac{p(z|y,\theta)}{
p(z|y,\hat\theta^{(i)})}p(z|y,\hat\theta^{(i)})dz\le 0.
\end{eqnarray*}\\\vspace{0mm}\pause

\item 最后的不等式是由于~$\log x\le x-1.$ 因此\vspace{-2mm}
$$
H(\hat\theta^{(i+1)},\hat\theta^{(i)})\le
H(\hat\theta^{(i)},\hat\theta^{(i)}).\vspace{-5mm}
$$\pause

\item 由~(5.3.2) 可见对任何~$i$ 有\vspace{-2mm}
$$\ln\pi(\hat\theta^{(i+1)}|y)\ge \ln\pi(\hat\theta^{(i)}|y).\vspace{-2mm}
$$
\end{itemize}
}

\frame{\frametitle{例子}
\begin{itemize}

\item 因此从任何初始值出发， E-M
算法一般都能达到一个局部最大值.\medskip\pause

\item {\bf 例~5.3.1} ~(基因连锁模型) ~考虑~Rao~(1973)
关于特定基因组合率的数据~(下表), 详细的介绍 可以参看~Sorensen 和
~Gianola (2002). 这里~197 个观测值被分为~4
个类，以及相应的概率.\vspace{1mm} \pause
\begin{table}[ht]
\centering {\bf \it 表~5.3.1 基因连锁数据}\vspace{2mm}\\
% \vspace{1pt}
 \begin{tabular}{l|l}
   \hline
   % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
  ~~~~ 计数\qquad\qquad~~~~& ~~~~ 概率\qquad\qquad~~~~ \\
   \hline
   ~~~$y_1=125$& ~~~$\frac12+\frac{\theta}{4}$ \\
   ~~~$y_2=18$ & ~~~$\frac14(1-\theta)$ \\
   ~~~$y_3=20$ & ~~~$\frac14(1-\theta)$ \\
   ~~~$y_4=34$ & ~~~$\frac{\theta}{4}$ \\
   \hline
 \end{tabular}
\end{table}
\end{itemize}
}

\frame{\frametitle{基因连锁数据的~E-M 算法迭代结果}
\begin{itemize}

\item 表~(5.3.2) 给出从初始~$\hat\theta^{(0)}=0.5$ 出发的迭代结果.\pause

 \begin{table}[ht]
 \centering{\bf\it 表~5.3.2 基因连锁数据的~E-M 算法迭代结果}\\
 \vspace{6pt}
 \begin{tabular}{c|c}
 \hline
 ~~\qquad 迭代次数\qquad~~ & ~~\qquad$\hat\theta^{(i)}$\qquad~~ \\
    \hline
   ~~~~~~~ 1 ~~~~~~& ~~~0.60825~~~~ \\
   ~~~~~~~ 2 ~~~~~~& ~~~0.62432~~~~ \\
   ~~~~~~~ 3 ~~~~~~& ~~~0.62648~~~~ \\
   ~~~~~~~ 4 ~~~~~~& ~~~0.62678~~~~ \\
   ~~~~~~~ 5 ~~~~~~& ~~~0.62682~~~~ \\
   ~~~~~~~ 6 ~~~~~~& ~~~0.62682~~~~ \\
\hline
\end{tabular}
\end{table}\pause

\item 即所得的估计~$\hat\theta=0.62682.$




\end{itemize}
}


\section{5.4 ~蒙特卡洛抽样方法}
\frame{\frametitle{ 引言}

\begin{itemize}

\item 如果所求的期望没有显式表达，那么除了可以使用分析逼近方法和数值积分方法之外，
蒙特卡洛抽样方法也是一个可选的计算方法.
%这种概率化的技巧在统计推断中是常用的工具.
为了估计总体均值或者总体分位数， 则从总体中抽取~(产生)
足够多的样本，然后使用样本均值或者样本分位数来估计相应的总体特征.
大数定律保证了 所得估计量是相合估计. \medskip\pause

\item 假设~$f$ 为一概率函数~(连续情形为密度函数，离散情形为概率分布), 以及
感兴趣的量为如下形式的有限期望\vspace{-2mm}
$$\qquad\qquad\quad E_fh(X)=\int_{\mathscr X}h(x)f(x)dx.\qquad\qquad\quad (5.4.1)\vspace{-2mm}
$$
\end{itemize}
}

\frame{\frametitle{引言}
\begin{itemize}

\item 如果可以从~$f$~中产生~$i.i.d$ 观测~$X_1,X_2,\cdots,
X_m$，则\vspace{-3mm}
$$\hspace{2.5cm}
\bar{h}_m=\frac{1}{m}\sum_{i=1}^mh(X_i)\hspace{3cm}(5.4.2)\vspace{-3mm}
$$
几乎处处收敛到感兴趣的量~$E_fh(X).$ 这一结果保证了在样本量~$m$
足够大时可以使用 ~$\bar{h}_m$ 作为~$E_fh(X)$~的估计.
这种用估计量~(5.4.2) 去逼近蒙特卡洛积分~(5.4.1)
的方法被称为蒙特卡洛抽样方法.\medskip\pause

\item
为给出估计精度或者逼近误差程度，我们可以使用类似的办法计算 标准差.
如果~$Var_f h(X)$~是有限的，则样本标准差是\vspace{-2mm}
$$s_m=\bigg\{\frac{1}{m-1}\sum\limits_{i=1}^m(h(X_i)-\bar{h}_m)^2\bigg\}^{1/2}.\vspace{-2mm}
$$
可以作为标准差的估计. 因此可以用~$s_m/\sqrt{m}$ 作为~$\bar{h}_m$
标准差的一个估计.
\end{itemize}
}

\frame{\frametitle{引言}
%\vspace{-2mm}

\begin{itemize}

\item 如果需要得到~$E_fh(X)$ 的置信区间，则由\vspace{-2mm}
$$
\frac{\sqrt{m}(\bar{h}_m-E_fh(X))}{s_m}\rightarrow
N(0,1)\vspace{-2mm}
$$
可以容易得到一个渐近水平为~$1-\alpha$ 的置信区间\vspace{-2mm}
$$\big[\bar{h}_m-z_{\alpha/2}s_m/\sqrt{m},\bar{h}_m+z_{\alpha/2}s_m/\sqrt{m}\big].\vspace{-2mm}
$$
这里~$z_{\alpha}$ 表示标准正态分布的上~$\alpha$
分位数.\medskip\pause

\item  以上的讨论表明，如果我们需要计算后验均值，则可以通过从后验分布中产生~$i.i.d$ 样本，
然后计算相应的样本均值作为估计.
但这种方法很少能直接使用，因为大多数情况下后验分布都不是标准分布而难以从中抽样.
\end{itemize}
}

\frame{\frametitle{一个例子}
%\begin{itemize}

%\item
\noindent{\bf 例~5.4.1~(续例~5.1.1)} ~由于\vspace{-2mm}
$$
E^\pi(\theta|x)=\frac{\dl\int_0^\infty \theta
exp\{-(\theta-x)^2/(2\sigma^2)\}(\tau^2+(\theta-\mu)^2)^{-1}d\theta}
{\dl\int_0^\infty
exp\{-(\theta-x)^2/(2\sigma^2)\}(\tau^2+(\theta-\mu)^2)^{-1}d\theta},\vspace{-2mm}
$$
因此~$E^\pi(\theta|x)$~为在正态分布~$N(x,\sigma^2)$~下\vspace{-1mm}
$$h_1(\theta)=\theta(\tau^2+(\theta-\mu)^2)^{-1}~~\mbox{和}~
~h_2(\theta)=(\tau^2+(\theta-\mu)^2)^{-1}\vspace{-1mm}
$$
两者的期望之比.
于是由前面的讨论知道，如果~$\theta_1,\cdots,\theta_m$
为从正态分布~$N(x,\sigma^2)$ 中产生的~$i.i.d$
样本，则~$E^\pi(\theta|x)$ 的估计量为\vspace{-2mm}
$$
\widehat{E^\pi(\theta|x)}=\frac{\dl\sum_{i=1}^m\theta_i(\tau^2+(\theta_i-\mu)^2)^{-1}}
{\dl\sum_{i=1}^m(\tau^2+(\theta_i-\mu)^2)^{-1}}.\vspace{-2mm}
$$
%\end{itemize}
}


\frame{\frametitle{一个例子}
\begin{itemize}

\item 问题并没有被完美解决. 由于从~$N(x,\sigma^2)$~中抽取的~$\theta$~集中在~$x$~附近，
并没有充分反应出柯西先验分布对后验分布的贡献，
而应当有显著的样本比例来自于后验分布的尾部.\medskip\pause

\item 因此如果把 ~$E^\pi(\theta|x)$~视为\vspace{-2mm}
$$\theta exp\{-(\theta-x)^2/(2\sigma^2)\}~~
\mbox{和}~~exp\{-(\theta-x)^2/(2\sigma^2)\}\vspace{-2mm}
$$
在柯西分布下的期望之比，则一个合适的估计量为\vspace{-2mm}$$
\widehat{E^\pi(\theta|x)}=\frac{\dl\sum_{i=1}^m\theta_i
exp\{-(\theta_i-x)^2/(2\sigma^2)\}}{\dl\sum_{i=1}^m
exp\{-(\theta_i-x)^2/(2\sigma^2)\}},\vspace{-2mm}
$$
其中~$\theta_1,\cdots,\theta_m$~为从~$C(\mu,\tau)$~中抽取的~$i.i.d$~样本.
\end{itemize}
}

\frame{\frametitle{蒙特卡洛重要性抽样}
\begin{itemize}

\item %这样也还是没有令人满意的解决问题.
由于后验分布没有像柯西分布那样的尾部，因此相对于后验分布
的中心而言，这种做法从尾部抽取了过多的样本.
这样就导致收敛速度变慢且在固定~$m$ 时逼近误差增大. 为了达到
满意的逼近应该直接从后验分布本身抽样.
为此，上述抽样方法的一个变种---{\bf
蒙特卡洛重要性抽样}方法被提出.\medskip\pause

\item 考虑式~(5.4.1).
假设从~$f$ 中直接抽样很困难，而从与~$f$ 很靠近的一个分布~$g$
中抽样比较容易. 那么我们可以将式~(5.4.1) 表示为\vspace{-3mm}
\begin{eqnarray*}
E_fh(X)\!\!\!&&=\int_{\mathscr X}h(x)f(x)dx
=\int_{\mathscr X}h(x)\frac{f(x)}{g(x)}g(x)dx\\
&&=\int_{\mathscr X}\{h(x)w(x)\}g(x)dx =E_g[h(x)w(x)],\vspace{-4mm}
\end{eqnarray*}
其中~$w(x)=f(x)/g(x).$
\end{itemize}
}

\frame{\frametitle{蒙特卡洛重要性抽样的例子} %\vspace{-3mm}

\begin{itemize}

\item 从~$g$ 中产生~$i.i.d$ 样本~$X_1,\cdots,X_m,$ 则一个合适的估计量为\vspace{-3mm}
$$
\widehat{hw}_m=\frac1m\sum_{i=1}^mh(X_i)w(X_i).\vspace{-4mm}
$$
抽样分布~$g$ 称为{\bf 重要性函数}.
我们用下例来说明重要性抽样方法的应用.\medskip\pause

\item {\bf 例~5.4.2} ~假设~$X_1,\cdots,X_n$~为从~$N(\theta,\sigma^2)$
中抽取的~$i.i.d$ 样本，其中~$\theta,~\sigma^2$ 均未知. 取 ~$\theta$
和~$\sigma^2$ 的先验为独立先验，其中~$\theta$
服从双指数先验分布，密度为~$\frac{1}{2}e^{-|\theta|/2}$，$\sigma^2$
有 先验密度~$(1+\sigma^2)^{-2}.$
这两个先验分布都不是标准先验，但都是从稳健性考虑而选取的. $\theta$
的后验均值为感兴趣的量，用蒙特卡洛重要性抽样方法求积分\vspace{-2mm}
$$
E^\pi(\theta|\bm{x})=\int_{-\infty}^\infty\int_0^\infty\theta\pi(\theta,\sigma^2|\bm{x})d\theta
d\sigma^2.\vspace{-2mm}
$$%\pause
\end{itemize}
}

\frame{\frametitle{蒙特卡洛重要性抽样的例子} %\vspace{-2mm}

\begin{itemize}

\item {\bf 解}
~由于~$\pi(\theta,\sigma^2|\bm{x})$
不是一个标准分布，我们先来寻找一个离它较近的分布.
记~$\bar{x}={\dl\frac{1}{n}}\sum\limits_{i=1}^nx_i,~s^2_n={\dl\frac{1}{n}}\sum\limits_{i=1}^n(x_i-\bar{x})^2.$
\smallskip\pause

\item 因此有\vspace{-2mm}
\begin{eqnarray*}
\pi(\theta,\sigma^2\!\!\!&&|\bm{x})\!\propto\!
(\sigma^2)^{-n/2}\exp\left\{\!-\frac{n}{2\sigma^2}\left[(\theta-\bar{x})^2\!+\!s_n^2\right]\right\}
\!\cdot\!\frac{\exp\{-|\theta|\}}{(1+\sigma^2)^{2}}\\
&&=\frac{[(\theta-\bar{x})^2+s_n^2]^{n/2+1}}{(\sigma^2)^{(n/2+2)}}
\exp\left\{-\frac{n}{2\sigma^2}\left[(\theta-\bar{x})^2+s_n^2\right]\right\}\\
&&~~~\times\left[(\theta-\bar{x})^2+s_n^2\right]^{-(n/2+1)}\exp\{-|\theta|\}
\Big(\frac{\sigma^2}{1+\sigma^2}\Big)^2\\
&&\propto g_1(\sigma^2|\theta,\bm{x})\cdot
g_2(\theta|\bm{x})\cdot\exp\{-|\theta|\}\Big(\frac{\sigma^2}
{1+\sigma^2}\Big)^2,
\end{eqnarray*}
\end{itemize}
}

\frame{\frametitle{蒙特卡洛重要性抽样的例子}
\begin{itemize}

\item 其中~$g_1$
为逆伽玛分布~$\Gamma^{-1}(\alpha,\beta),$ 其中~$
\beta=n[(\theta-\bar{x})^2+s_n^2]\big/2,$ $\alpha=n/2+1;$ $g_2$
为广义~$t$ 分布~$\mathscr{T}_1\big(n+1,\bar{x},s_n^2/(n+1)\big).$
\medskip\pause

\item 注意到~$e^{-|\theta|/2}\big[\sigma^2/(1+\sigma^2)\big]^2$ 的尾部
对~$g_1\cdot g_2$
没有太大影响，故可选择~$g(\theta,\sigma^2|\bm{x})=g_1(\sigma^2|\theta,\bm{x})
\cdot g_2(\theta|\bm{x})$ 作为重要性函数.\medskip\pause

\item 因此抽样时需要先从~$g_2$
抽取一个~$\theta,$ 然后在给定此~$\theta$ 的条件下，再从~$g_1$
抽取一个~$\sigma^2,$ 合在一起 组成一次抽样~$(\theta,\sigma^2).$
在抽取了~$m$ 对~$(\theta,\sigma^2)$ 后得到~$E^\pi(\theta|x)$
的一个估计\vspace{-2mm}
$$
\widehat{E^\pi(\theta|\bm{x})}=\dl\sum_{i=1}^m\theta_iw(\theta_i,\sigma^2_i|\bm{x})\bigg
/\dl\sum_{i=1}^mw(\theta_i,\sigma^2_i|\bm{x}),\vspace{-2mm}
$$
其中~$w(\theta,\sigma^2|\bm{x})=f(\bm{x}|\theta,\sigma^2)\pi(\theta,\sigma^2)/g(\theta,\sigma^2|\bm{x}).$
\end{itemize}
}


\frame{\frametitle{蒙特卡洛重要性抽样的例子}

\begin{itemize}

\item 关于例~5.4.1 中积分如何达到满意的逼近效果的问题还未解决.
对此问题使用蒙特卡洛重要性抽样方法,
很难找到离目标后验分布比较接近的重要性函数. 因此,
应从后验分布本身抽样. 由于后验分布不是标准分布，难以从中直接抽样.
我们可以利用将在下一节介绍的~MCMC
方法，以目标后验分布~$\pi(\theta|x)$
作为马氏链的稳定分布，生成随机数，获得积分较为满意的模拟结果.
这一问题留给读者作为练习.\medskip\pause

\item 在一些高维问题里，组合使用数值积分、Laplace
逼近和蒙特卡洛方法可能会给出很好的结果.
关于蒙特卡洛统计方法的非常好的介绍可以参考~Robert 和~Casella~(1999).
\end{itemize}
}



\section{5.5 ~马尔科夫链蒙特卡洛方法}

\subsection{5.5.0 ~引言}

\frame{\frametitle{引言}
\begin{itemize}

\item
上一节介绍的蒙特卡洛抽样方法主要使用于后验分布是一维，且是标准分布的情形，
此时从中抽样比较容易.\medskip\pause

\item  它存在如下缺点：(1)
如后验分布不是标准分布使用蒙特卡洛方法中的~\!``重要性函数", 选择一个
合适的~\!``重要性函数", 如果选择不合适, 则蒙特卡洛逼近的效果会很差.
特别在后验分布是高维情形, 很难找到 合适的
~\!``重要性函数".\smallskip\pause

(2)~标准的蒙特卡洛方法一个严重缺点是在实施中后验分布的完全形式必须已知.
对那些后验分布不完全指定或者不直接指定的场合就不能处理了.
\medskip\pause

\item 在上述情形下如何从后验分布获取样本？此时马尔科夫链~(Markov Chain, 简记~MC) 将被用来解决这个问题.
以目标后验分布作为其平稳分布的马尔科夫链生成随机数,
%代替从后验分布中抽取的样本，
利用蒙特卡洛~(也简记为~MC) 逼近方法，获得相应积分的模拟结果.
\end{itemize}
}

\subsection{5.5.1 ~MCMC 中的马尔科夫链}
\frame{\frametitle{马尔科夫链的定义}
\vspace{-2mm}

\begin{itemize}

\item {\bf 定义~5.5.1}
~设~$\{X_n,~n\ge0\}$ 是取有限个或可列个值的随机过程，若~$X_n=i,$
表示过程在时刻~$n$ 状态处于~$i,~S=\{0,1,\cdots\}$ 为状态集.
若对一切~$n$ 有\vspace{-1mm}
\begin{eqnarray*}
&&P\big(X_{n+1}=j~\big|~\!X_0=i_0,X_1=i_1,\cdots,X_{n-1}=i_{n-1},X_n=i\big)\\
&&=P\left(X_{n+1}=j~\big|~\!X_n=i\right)\vspace{-2mm}
\end{eqnarray*}
则称~$\{X_n,n\ge0\}$
是离散时间马尔科夫链，常简称为马氏链.\medskip\pause

\item 由定义可知，对过程~$\{X_n,~n\ge0\}$ 的
将来状态~$\{X_{n+1}\}$ 只与现在状态~$\{X_n\}$ 有关,
而与过去状态~\!$\{X_k,k\le n-1\}$\! 无关.\medskip\pause

\item 条件概率~$P(X_{n+1}=j|X_n=i)$ 称为马氏链的一步转移概率，
若转移概率与~$n$
无关，为固定值，则称马氏链有平稳转移概率，记为~$p_{ij}.$
具有平稳转移概率的马氏链它也称为时间齐性马氏链.
~$P=\big(p_{ij}\big)~\forall~i,j\in S$ 称为马氏链的转移概率矩阵,
满足条件~$p_{ij}\geq 0,~~\mbox{且}~ \sum\limits_{i=0}^\infty
p_{ij}=1.$
\end{itemize}
}

\frame{\frametitle{ 1. ~马氏链的平稳性}
\begin{itemize}

\item  {\bf 定义~5.5.2} ~设马尔科夫链有
 转移概率阵~$P=(p_{ij}),$ 一个概率分布~$\pi=\{\pi_i,~i\geq 0\}$ 如果满足
~$\pi_j=\sum\limits_i \pi_i p_{ij},$
则称之为此马尔科夫链的平稳分布.\medskip\pause

\item 易看出，如果过程初始状态~$X_0$ 有平稳分布~$\pi=\{\pi_i,~i\geq 0\},$
即~$P(X_0=j)=\pi_j.$ 则有\vspace{-3mm}
$$
P(X_1=j)=\sum_iP(X_1=j|X_0=i)P(X_0=i) =\sum_i\pi_ip_{ij}%=P(X_0=j)
=\pi_j.\vspace{-5mm}
$$\pause

\item 由归纳法可得\vspace{-2mm}
$$P(X_n=j)=\sum_i P(X_n=j|X_{n-1}=i)P(X_{n-1}=i)
=\sum_i \pi_i p_{ij}.\vspace{-3mm}
$$
于是对所有的~$n,~X_n$ 有相同的分布~$\pi$. 即~$\{X_n,n\ge0\}$
作为随机过程是平稳的.
\end{itemize}
}

\frame{\frametitle{ 2. ~不可约性}
\begin{itemize}

\item {\bf 定义~5.5.3} ~一个具有可数状态空间~$S$ 和转移概率矩阵~$P=$ $(p_{ij})$ 的马氏链~$\{X_n\}$ 称为是
不可约的, 如果对任意两个状态~$i,~j\in S$，此链从状态~$i$
出发转移到状态~$j$ 的概率为正的.\medskip\pause

\item 即对某个~$n\ge1$ 有\vspace{-2mm}
$$
p_{ij}^{(n)}=P(X_n=j~|X_0=i)>0.\vspace{-4mm}
$$\pause

\item 由定义可知，具有“不可约性”的马氏链意味着从任一状态出发总可到达任一其它状态.
\end{itemize}
}


\frame{\frametitle{ 3. ~非周期性和正常返性}
\begin{itemize}

\item {\bf 定义~5.5.4} ~一个马氏链的一个状态~$i$ 有周期~$k,$
如果经过~$k$ 的倍数步后一定可以返回到状态~$i,$ 即\vspace{-2mm}
$$
k(i)=gcd\{n:~P\big(X_n=i~\big|~\! X_0=i)>0\big),\vspace{-2mm}
$$
其中~$gcd$ 表示“最大公约数”.
如果返回任一状态的次数的最大公约数是~$1,$
则称此马氏链是非周期的.\medskip\pause


\item 非周期的马氏链可以保证其不会陷入循环当中.\medskip\pause

\item {\bf 定义~5.5.5} ~对常返状态~$i,$ 令~$T_i=\inf\big\{n\geq
1:~X_n\!=\!i~\big|~\!X_0\!=\!i\big\}$ 为首次返回状态~$i$
的时刻，如果\vspace{-2mm}
$$
\mu_i=E(T_i)<\infty,\vspace{-2mm}
$$
则称状态~$i$ 是正常返的；若~$\mu_i=\infty$ 时称状态~$i$ 是零常返的.
\end{itemize}
}

\frame{\frametitle{ 4. ~遍历性}
\begin{itemize}

\item {\bf 定义~5.5.6} ~一个马氏链的状态称为遍历的,
如果它是非周期且正常返的. 如果马氏链的所有状态都是遍历的，则称此
马氏链是遍历的.\medskip\pause

\item ~综上所述，由马氏链的基本理论可知，我们需要构造的马氏链必须是不可约、正常返和非周期的.
满足这些正则条件的马氏链存在唯一的平稳分布.

\end{itemize}
}


%\subsection{5.5.2 马氏链的极限定理}
\frame{\frametitle{马氏链的极限定理}
\begin{itemize}

\item {\bf 定理~5.5.1} ~设~$\{X_n,~n\geq
0\}$ 为一具有可数状态空间~$S$ 的马氏链，其转移概率矩阵为~$P$. 进一步
假设它是不可约、非周期，有平稳分布~$\pi=(\pi_i: i\in S\},$
则有\vspace{-2mm}
$$
\sum_{j\in S}~\!\Big|P(X_n=j)-\pi_j\Big|\rightarrow 0,~~
n\rightarrow\infty, \vspace{-2mm}
$$
对~$X_0$ 的任意初始分布~$\pi$.\medskip\pause

\item 换言之，对比较大的~$n$，$X_n$ 的分布将会接近~$\pi$.
对一般的状态空间，类似的结果
也存在：在合适的条件下，当~$n\rightarrow\infty$ 时~$X_n$
的分布将收敛到~$\pi$.
\end{itemize}
}

\frame{\frametitle{马氏链的大数定律}
\begin{itemize}

\item  {\bf 定理~5.5.2}~(马氏链的大数定律) ~假设~$\{X_n,~n\geq
0\}$ 为一具有可数状态空间~$S$ 的马氏链，其转移概率矩阵为~$P$. 进一步
假设它是不可约的且有平稳分布~$\pi=(\pi_i: i\in S\}$.
则对任何有界函数~$h: S\rightarrow \mathbb{R}$ 以及初始值~$X_0$
的任意初始分布有\vspace{-2mm}
$$
\frac1n\sum_{i=0}^{n-1}h(X_i)\rightarrow \sum_j h(j)\pi_j,\quad
n\rightarrow\infty,\vspace{-2mm}
$$
依概率成立.\medskip\pause

\item 当状态空间为不可数，马氏链~$\{X_n,~\!n\geq
0\}$ 为不可约且有平稳分布~$\pi$ 时，也有\vspace{-2mm}
$$\frac1n\sum_{i=0}^{n-1}h(X_i)\rightarrow \int_S h(x)d\pi(x),\quad
n\rightarrow\infty,%\vspace{-2mm}
$$
\end{itemize}
}

\frame{\frametitle{极限定理的应用}
\begin{itemize}

\item 这个定理结论是非常有用的.
比如给定集合~$S$ 上的概率分布~$\pi$，以及~$S$
上的实函数~$h(\theta).$\medskip\pause


\item 设要计算积分~$\mu=\dl\int_S
h(\theta)\pi(\theta|\mathbf{x})d\theta,$
当从后验分布~$\pi(\theta|\mathbf{x})$ 中难以直接抽样时，
则可以构造一个马氏链，使得其状态空间为~$S$ 且其平稳 分布~$\pi$
就是目标后验分布~$\pi(\cdot|\mathbf{x}),$ 从一初始值~$\theta_0$
出发，将此链运行一段时间，如~$0,1,2,\ldots,n-1,$ 生成随机数(样本)
$\theta_0,\theta_1,\cdots,\theta_{n-1},$\medskip\pause


\item 由马氏链的大数定律可知\vspace{-3mm}
$$
\bar{\mu}_n=\frac1n\sum_{j=0}^{n-1}h(\theta_j)\vspace{-3mm}
$$
为所要求积分~$\mu$ 的一相合估计.
这种技术称为马尔科夫链蒙特卡洛~(MCMC) 方法.
\end{itemize}
}


\subsection{5.5.2 ~MCMC 的实施}

\frame{\frametitle{1. ~MCMC 的实施的一般步骤}

\noindent 当从后验分布~$\pi(\theta|\mathbf{x})$
中难以直接抽样时，我们需要构造一个不可约、非周期的马氏链使其平稳
分布就是目标后验分布~$\pi(\cdot|\mathbf{x}),$
按下列步骤生成所需要的样本：\smallskip\pause

(1) ~选择初始值$\theta_0.$\smallskip\pause

(2) ~产生~$J$ 值, 迭代总次数，使马氏链达到平稳分布.\smallskip\pause

(3) ~监控算法的收敛性，若收敛性诊断失败， 则加大~$J$
生成更多的观测值.\smallskip\pause

(4) ~去除链开头的~B 个观测值~(即预烧期~\!), 将链后面的~$T=$ $J-B$
个观测值~$\theta_{B+1},\theta_{B+2},\cdots, \theta_J$
作为后验分析的样本.\smallskip\pause

(5) ~作后验分布图~(通常作一元边缘分布图).\smallskip\pause

(6) ~输出后验分布的~Summary
~(均值、中位数、标准差、分位数，相关系数等),
获得相应的蒙特卡洛积分的模拟结果. }

\frame{\frametitle{2. ~MCMC 实施中的若干术语--初始值}
\begin{itemize}

\item {\bf 初始值} ~初始值
被用来初始化一个马氏链. 如果初始值远离
后验密度的最高区域，且算法的迭代次数大小不足以消除初始值的影响，则它对后验推断可能会造成影响.
\medskip\pause

\item 我们可以通过一些方式降低或者避免初始值的影响，比如去掉开始一段时间的迭代值、或者从不同的初始值出发
获得抽样等等.\medskip\pause

\item 合理的初始值可以是靠近后验分布的中心位置或者似然函数的最大值点，但是靠近似然
函数的最大值点在一些场合下已经被证明不是一个很好的选择.\medskip\pause

\item 如果先验分布
是有信息的，则也可以选择先验分布的期望或者众数作为初始值.
一般地，选择多个从不同初始值开始的链仍然是 最推荐的做法.

\end{itemize}
}

\frame{\frametitle{MCMC 实施中的若干术语--预烧期和抽样步长}
\begin{itemize}

\item 在~MCMC 机制中用以保证链达到平稳状态所运行的时间称为是~{\bf 预烧期}, 其迭代数记为~$B.$
\medskip\pause

\item 为避免初始值的影响，预烧期中~$B$ 个迭代值将被从样本中去除.
只要链运行的时间足够长，去除预烧期对后验推断是几乎没有影响的.\medskip\pause

\item 显然马氏链产生的样本并不是相互独立的，如果需要独立样本，则我们可以通过监视
产生样本的自相关图，然后选择~{\bf 筛选间隔或抽样步长}~$L>1$ 使得~$L$
步长以后的自相关性很低.\medskip\pause

\item 我们可以通过每间隔~$L$ 个样本抽取一个来获得~(近似) 独立样本.
\end{itemize}
}

\frame{\frametitle{MCMC 实施--迭代保持数和算法的收敛性}
\begin{itemize}

\item {\bf 迭代保持数~$T$} ~设迭代总次数为~$J,$
迭代保持数为迭代总次数去掉预烧期迭代次数后，提供给后验贝叶斯分析用的实际样本数，故有
~$T=J-B.$ 如果考虑一个抽样步长~$L,$
则迭代保持数为去掉预烧期后，最终的~(近似) 独立样本数，此时~$T=$
$(J-B)/L.$\medskip\pause

\item {\bf 算法的收敛性} ~是指所得到的链是否达到了平稳状态.
如果达到了平稳分布，则我们得到的样本可以认为是从目标分布中
抽取的样本.\medskip\pause

\item 一般而言，我们并不清楚必须运行算法多长时间才能认为所
得到的链达到了平稳分布. 因此监视链的收敛性是~MCMC
计算方法中的本质问题.
\end{itemize}
}

\frame{\frametitle{MCMC 实施中的若干术语--蒙特卡洛~(MC) 误差}
\begin{itemize}

\item 在~MCMC 输出结果分析中，一个必须报告和监控的量就是{\bf 蒙特卡洛误差}. 蒙特卡洛误差
度量了每一个估计因为随机模拟而导致的波动性.\medskip\pause

\item 在计算感兴趣的参数时，
其精度应该随着样本量递增，
因而蒙特卡洛误差必然很低，它和样本量大小成反比并且用户自己可以控制.
因此增加迭代次数，感兴趣的量的估计精度也会增加.\medskip\pause

\item
常用的估计蒙特卡洛误差的方法有两种：组平均~(batch mean) 方法和
窗口估计量~(window estimator) 方法.
第一种方法简单容易操作，但是第二种方法更精确.
\end{itemize}
}



\frame{\frametitle{3. ~利用~MCMC 的输出结果描述目标后验分布}
\begin{itemize}

\item MCMC 输出结果给我们提供了一类随机样本
$$\theta_1,\theta_2,\cdots,\theta_t,\cdots,\theta_T.$$
这组样本可以看成从目标后验分布~$\pi(\theta|\mathbf{x})$
中生成的样本.\medskip\pause

\item 设~$g(\theta)$ 是感兴趣的参数~$\theta$ 的函数，从这组样本我们可以进行下列工作：\medskip\pause

\item (1) ~获得~$g(\theta)$ 的后验~Summary: ~后验均值、后验标准差、
后验中位数、后验分位数等.
例如，利用传统的估计方法和上述样本可分别获得~$g(\theta)$
的后验期望估计\vspace{-2mm}
$$
\widehat{E}\big(g(\theta)\big|\mathbf{x}\big)=\overline{g(\theta)}=\frac{1}{T}
\sum_{t=1}^T g\big(\theta_t\big),\vspace{-2mm}
$$

\end{itemize}
}

\frame{\frametitle{利用~MCMC 的输出结果描述目标后验分布}
\begin{itemize}

\item 以及后验标准差的估计\vspace{-2mm}
$$\widehat{SD}\big(g(\theta)\big|\mathbf{x}\big)=\bigg\{\frac{1}{T-1}
\sum_{t=1}^T
\left[g\big(\theta_t\big)-\widehat{E}\big(g(\theta)\big|\mathbf{x}\big)\right]^2\bigg\}^{1/2}.\vspace{-2mm}
$$
最后，还可以获得后验众数的估计，即在后验密度图中使后验密度达到最大的点.

\item (2)
~利用后验分位数，例如~$2.5\%$ 和~$97.5\%$
分位数，构造可信水平~$95\%$ 的可信区间.\medskip\pause

\item (3) ~计算和监测参数间的相关性.\medskip\pause

\item (4) ~产生边缘后验分布图~(直方图、密度函数图等).

\end{itemize}
}

\frame{\frametitle{4. ~MCMC 方法收敛性的诊断}
\begin{itemize}

\item 无论使用哪一种抽样方法，都要确定所得到的马氏链的收敛性，
即需要确定马氏链达到收敛状态时迭代的次数~(达到
收敛状态前的那一段链称为“预烧期”样本).\medskip\pause

\item 通常没有一个全能的方法确定马氏链的收敛性，监视链的收敛性有许多方法.
但是每种方法都是针对收敛性问题的不同方面提出的.\medskip\pause

\item 因此，在绝大数情况下，为了保证链的收敛性必须应用几种不同的方法去诊断.
下面将介绍几种常用的诊断方法.

\end{itemize}
}

\frame{\frametitle{ (1) ~蒙特卡洛~(MC) 误差}
\begin{itemize}

\item 诊断马氏链收敛性的最简单的方法就是监视~MC 误差.\medskip\pause

\item 因为较小的~MC
误差表明在计算感兴趣的量时精度较高.\medskip\pause

\item 因此~MC
误差越小表明马氏链的收敛性越好.

\end{itemize}
}

\frame{\frametitle{(2) ~样本路径图}
\begin{itemize}

\item 另外一种监控方式是使用样本路径图~(trace
plot)： 马氏链迭代次数对生成的值作图.
如果所有的值都在一个区域里且没有明显的周期性和趋势性，那么我们
可以假设收敛性已经达到.\medskip\pause

\item 为避免链陷入目标分布的某个局部区域，通常作几个平行的链，
它们的初始值非常分散.
在经过一段时间后，如果他们的样本路径图都稳定下来，而且混合在一起无法区别，这时可以
判定收敛性已经达到.\medskip\pause

\item 这个想法可以通过将多个链的样本路径图画在同一个图上来检查.
图~5.5.1 给出了一个明显没有达到收敛的例子. 图~5.5.2
则看起来更令人相信链达到了平稳分布，波动比较稳定，没有明显的
周期性和趋势性.
\end{itemize}
}

\frame{\frametitle{图~5.1} %\vspace{-2mm}

\begin{figure}[h,p,t]\vspace{-2mm}
\centering{
\includegraphics[scale=0.32,angle=0]{fig51.jpg}}\\ \vspace{0mm}
\centerline{\bf\it 图~5.1
~没有达到收敛的链的路径图}\label{fig51}%\vspace{5mm}
\end{figure}
}



\frame{\frametitle{图~5.2} \vspace{-2mm}

\begin{figure}[h,p,t]\vspace{-2mm}
\centering{
\includegraphics[scale=0.32,angle=0]{fig52.jpg}}\\ \vspace{0mm}
\centerline{\bf\it 图~5.2 ~达到收敛的链的路径图}\label{fig2}%\vspace{5mm}
\end{figure}
}

\frame{\frametitle{ (3) ~遍历均值图}
\begin{itemize}

\item 还有一种很有用的图方法是将马氏链的累积均值对迭代次数作图就得到此链的遍历均值图
~(ergodic mean plot).\medskip\pause

\item 这里累积均值
是指此量直至当前迭代的平均值.\medskip\pause

\item 如果累积均值在经过一些迭代后基本稳定，则表明算法
已经达到收敛~(见图~5.5.3).
\end{itemize}
}

\frame{\frametitle{图~5.3} \vspace{-2mm}

\begin{figure}[h,p,t]\vspace{-2mm}
\centering{
\includegraphics[scale=0.32,angle=0]{fig53.jpg}}\\ \vspace{-2mm}
\centerline{\bf\it 图~5.3 ~达到收敛的链的路径图}\label{fig53}%\vspace{5mm}
\end{figure}
}


\frame{\frametitle{(4) ~自相关函数~(ACF) 图和~Gelman-Rubin 方法}
\begin{itemize}

\item 诊断马氏链收敛性通过监视自相关函数图
~(Autocorrelations function plot) 也是很有用的.\medskip\pause

\item 链的迭代次数对~ACF 作图，因为较低或者
较高的自相关性分别表明了马氏链的快或慢的收敛性.\medskip\pause

\item Gelman-Rubin 方法是一种常用的收敛性检验方法.
\end{itemize}
}

\frame{\frametitle{(5) ~Gelman-Rubin方法}
\begin{itemize}

\item 许多统计检验工具也被开发出来用于收敛诊断~(Cowles and Carlin, 1996;
Brooks and Roberts, 1998).  CODA~(Best et al., 1996) 和~BOA~(Smith,
2005) 软件程序也被开发用于实施这些工具.\medskip\pause

\item Gelman-Rubin 方法是一种常用的收敛性检验方法. 详细介绍参看~Gelman \& Rubin
(1992) 或~Rizzo (2008). 参看本教材~203-206 页的介绍.
\end{itemize}
}


\subsection{5.5.3 ~M-H 算法}
\frame{\frametitle{引言}
\begin{itemize}

\item 从一般的后验分布中进行抽样的流行做法是采用马尔科夫蒙特卡洛~(MCMC) 方法.
MCMC~抽样策略是建立一个不可约的、非周期的马氏链，并且其平稳分布即为感兴趣的目标后验分布.
\medskip\pause

\item 核心问题是确定从当前值转移到下一个值的规则.
其主要任务是生成满足上述正则条件的一个
马氏链~$\big\{x_i,~i=1,\cdots,J\big\}. $\medskip\pause

\item Metropolis-Hastings 算法是一类最常用的~MCMC 抽样方法, 它首先由~Metropolis~(1953)
提出, 后来由~Hastings~(1970) 加以推广的.\medskip\pause

\item 本节我们介绍~Metropolis-Hastings
抽样方法及其几个变种，它们在贝叶斯推断中有着广泛的应用.
\end{itemize}
}

\frame{\frametitle{Metropolis-Hastings~(Ｍ-Ｈ) 抽样方法}
\begin{itemize}

\item 为简化记号,
在本节以下用~$f(\cdot)$ 表示目标分布，以~$g(\cdot)$
表示提议分布.\medskip\pause

\item 设我们希望从目标分布~$f(\cdot)$ 中抽样，M-H
算法从初始值~$x_{\D 0}$ 出发，指定一个从当前值~$x_t$
转移到下一个值~$x_{t+1}$
的规则，从而产生马氏链~$\{x_i,~i=0,1,\cdots\}.$\medskip\pause

\item 具体来说， 在给定当前值~$x_t,$ 从一个提议分布~$g(\cdot|x_t)$
产生一个候选点~$x',$ 若候选点~$x'$ 被接受，则链从状态~$x'$
转移到链的下一时刻~$t+1,$ 令~$x_{t+1}=x';$ 否则链停留在状态~$x_t,$
令~$x_{t+1}=x_t.$
\end{itemize}
}

\frame{\frametitle{~M-H 算法}
\begin{itemize}

\item 候选点~$x'$ 是否被接受为链的下一个值，根据下列接受概率来确定：$\alpha(x_t,x')=\min\{1,A\},$
其中\vspace{-2mm}
$$
A=\frac{f(x')g(x_t|x')}{f(x_t)g(x'|x_t)}. \vspace{-2mm}
$$\pause

\item 需要注意的一点是上式分子和分母中的密度函数可分别用“密度函数的核”代替，故密度函数中的正则化常数
因子可省略掉, 以便简化计算.\medskip\pause

\item ~M-H 算法生成满足正则条件的马氏链迭代过程如下：

\end{itemize}
}

\frame{\frametitle{~M-H 算法生成马氏链迭代过程}
\begin{itemize}

\item (1) ~选择一个提议分布~$g(\cdot|x_t);$\medskip\pause

\item (2) ~从提议分布中生成初始值~$x_{\D 0};$\medskip\pause

\item (3) ~对~$t=1,2,\cdots$ 重复下列步骤\medskip\pause

~~~~~(a) ~从提议分布~$g(\cdot|x_t)$
中产生一个候选值~$x'.$\medskip\pause

~~~~~(b) ~从均匀分布~$U(0,1)$ 中生成随机数~$U.$\medskip\pause

~~~~~(c) ~计算接受概率, 若~$U\leq A$ 则接受~$x'$ 且令~$x_{t+1}=x',$
否则令~$x_{t+1}=x_t.$\medskip\pause

~~~~~(d) ~增加~$t,$ 返回到~(a)
\end{itemize}
}

\frame{\frametitle{~M-H 方法的一个注释}
\begin{itemize}

\item Robert 和~Casella (1999)
指出，接受概率并非越大越好，因为这可能导致较慢的收敛性. Gelman 等
~(1996) 建议当参数维数是~$1$ 时，接受概率应略小于~$0.5$
是最优的，当维数大于~$5$ 时，接受概率应降至~$0.25$
左右.\medskip\pause

\item 为说明~M-H
抽样方法产生的马氏链具有平稳分布~$f$，则可以通过说明此链的转移核~(或转移概率)
和~$f$ 一起满足细致平衡方程. 此结论的证明见~Robert 和~Casella~(2004)
定理~7.2, 或参看本教材~207-208 页的内容.
\end{itemize}
}

\frame{\frametitle{提议分布应满足的条件}
\begin{itemize}

\item 提议分布~$g$ 的选择除了使得产生的马氏链满足不可约, 正常返, 非周期
且具有平稳分布~$f$ 等正则化条件外，还应满足：\medskip\pause

\item (1) ~提议分布的支撑集包含目标分布的支撑集.\medskip\pause

\item (2) ~容易从中抽样, 常取为已知的分布，如正态或~$t$ 分布等.\medskip\pause

\item (3) ~提议分布应使接受概率容易计算.\medskip\pause

\item (4) ~提议分布的尾部要比目标分布的尾部厚.\medskip\pause

\item (5) ~新的候选点被拒绝的频率不高.\medskip\pause
\end{itemize}
}

\frame{\frametitle{目标分布为后验分布时的~M-H 算法}
\vspace{-2mm}

\begin{itemize}

\item 前面所述的~M-H
算法可直接在贝叶斯框架下实现，只要将~$x$ 用感兴趣的参数~$\theta$
代替，将目标分布~$f(x)$ 用后验分布~$\pi(\theta|\mathbf{x})$ 代替.
%因此，在贝叶斯分析中，M-H算法产生马氏链~$\big\{\theta_i,~i=1,2,\cdots\big\}$ 的
方法概括如下：\medskip\pause

\item (1) ~选择一个提议分布~$g(\cdot|\theta_t);$\medskip\pause

\item (2) ~从提议分布中生成初始值~$\theta_0;$\medskip\pause

\item (3) ~对~$t=1,2,\cdots$ 重复下列步骤\medskip\pause

~~~~~(a) ~从提议分布~$g(\cdot|\theta_t)$
中产生一个候选值~$\theta'.$\medskip\pause

~~~~~(b) ~从均匀分布~$U(0,1)$~中生成随机数~$U.$\medskip\pause

~~~~~(c) ~计算接受概率, 若\vspace{-2mm}
   $$
   U\leq\frac{\pi(\theta'|\mathbf{x})g(\theta_t|\theta')}{\pi(\theta_t|\mathbf{x})g(\theta'|\theta_t)},
   \vspace{-1mm}
   $$
\hspace{1.5cm}则接受~$\theta'$ 且令~$\theta_{t+1}=\theta',$
否则令~$\theta_{t+1}=\theta_t.$\medskip\pause

~~~~~(d) ~增加~$t,$ 返回到~(a).
\end{itemize}
}

\frame{\frametitle{M-H 算法的例子}
\begin{itemize}

\item 下例说明如何利用~M-H 算法生成~Rayleigh (瑞利) 分布的样本.\medskip\pause

\item {\bf 例~5.5.2}
~使用~Metropolis-Hastings 抽样方法从瑞利分布中抽样.
瑞利分布的密度为\vspace{-2mm}
$$
f(x)=\frac{x}{\sigma^2}e^{-x^2/(2\sigma^2)}, ~x\ge 0,
~\sigma>0.\vspace{-2mm}
$$
瑞利分布适用于使用寿命受到快速老化的建模.\medskip\pause

\item {\bf 解} ~取自由度~(df) 为~$X_t$ 的~$\chi^2$ 分布为提议分布,
则使用~M-H 抽样方法的步骤如下:\medskip\pause

\item (1) ~令~$g(\cdot|X)$ 为~$\chi^2(df=X)$.
\end{itemize}
}

\frame{\frametitle{续例~5.5.2}
\begin{itemize}

\item (2) ~从~$\chi^2(1)$ 中产生~$X_0$, 并存在~x[1] 中.\medskip\pause

\item (3) ~对~$i=2,\cdots,N,$ 重复\medskip\pause

~~~~(a) ~从~$\chi^2(df\!=\!X_t)=\chi^2(df\!=\!\mbox{x}[i-1])$
中产生~$Y$.\medskip\pause

~~~~(b) ~产生~$U\sim U(0,1)$.\medskip\pause

~~~~(c) ~由~$X_t=\mbox{x}[i-1]$, 计算\vspace{-2mm}
  $$
  r(X_t,Y)=\frac{f(Y)g(X_t|Y)}{f(X_t)g(Y|X_t)},\vspace{-2mm}
  $$
  其中~$f$ 为瑞利密度. $g(Y|X_t)$ 为~$\chi^2(df=X_t)$ 的密度在~$Y$ 处的值, $g(X_t|Y)$
  为~$\chi^2(df=Y)$ 的密度在~$X_t$ 处的值. 若~$U\le r(X_t,Y),$ 则接受~$Y,$
  令~$X_{t+1}=Y;$
  否则令~$X_{t+1}=X_t$. 将~$X_{t+1}$ 存在~x[i] 里.\medskip\pause

~~~~(d) ~增加~$t,$ 返回到~(a).
\end{itemize}
}

\frame{\frametitle{续例~5.5.2}
\vspace{-2mm}

\begin{itemize}

\item
在密度~$f$ 中的常数可以在计算~$r$ 中抵消, 因此\vspace{-2mm}
$$
r(x_t,y)\!=\!\frac{f(y)g(x_t|y)}{f(x_t)g(y|x_t)}\!=\!\frac{ye^{-y^2/2\sigma^2}}{x_te^{-x_t^2/2\sigma^2}}\!\times\!
\frac{\Gamma(x_t/2)2^{x_t/2}x_t^{y/2-1}e^{-x_t/2}}{\Gamma(y/2)2^{y/2}y^{x_t/2-1}e^{-y/2}}.\vspace{-2mm}
$$
在此例中, 我们还是通过计算整个密度在某点的值来计算~$r$.
下面的代码计算瑞利密度在某点的值:\medskip\pause

\item
下面我们产生~$\sigma=4$ 的瑞利分布随机数.
使用的提议分布为自由度是~$x_t=\mbox{x}[i-1]$ 的 ~$\chi^2(df=x_t)$
分布:\medskip\pause

\item 在计算~$r(X_{t},Y)$ 中, 分子和分母分别用变量~num 和~den 表示. 记数变量~k 记录了
候选点被拒绝的次数~$k=4070,$ 大约~40\% 的候选点被拒绝了.
因此这个方法产生链的效率不高. 我们使用样本对时间作图~(称为~trace
plot), 来观测其样本路径图~(见图~5.5.1).
注意在候选点被拒绝的时间点上链没有移动, 因此图中有很多短的水平平移.
\end{itemize}
}

\frame{\frametitle{图~5.5.1}

\begin{figure}[h,p,t]\vspace{-7mm}
\centering{
\includegraphics[scale=0.52,angle=0]{fig551.pdf}}\\ \vspace{-3mm}
\centerline {\bf\it
图~\!5.5.1 ~瑞利分布由~M-H 算法产生的马氏链部分路径图}\label{fig551}%\vspace{5mm}
\end{figure}
}



\frame{\frametitle{续例~5.5.2}
\vspace{-2mm}

\begin{itemize}

\item 本例中我们的目的是说明~M-H 抽样方法的应用，
对瑞利分布，有更高效率的产生随机数方法. 比如瑞利分布的分位数
可以表示为
$$x_q=F^{-1}(q)=\sigma[-2\ln(1-q)]^{1/2}, 0<q<1.$$
因此可以使用逆变换方法生成随机数.\medskip\pause

\item 下面比较~$\sigma=4$ 的瑞利分布的理论分位数和生成链的分位数拟合程度的~QQ 图，
和以~$\sigma=4$ 的瑞利分布
为目标分布的链生成样本的直方图.\medskip\pause

\item 图~5.5.2 的左边是样本的直方图和密度函数曲线图，右边的~QQ
图是判断由链产生的样本分位数和目标分布的理论分位数拟合好坏的一种方法.
由于~QQ
图上的点基本上集中在一条直线附近，这显示样本分位数和理论分位数是高度近似一致的.
\end{itemize}
}

\frame{\frametitle{图~5.5.2}
\begin{figure}[h,p,t]\vspace{-7mm}
\centering{
\includegraphics[scale=0.42,angle=0]{fig552.pdf}}\\ \vspace{0mm}
\centerline %\caption
{\bf\it 图~\!5.5.2 ~M-H 链样本的直方图和~M-H 链分位数的~QQ
图}\label{fig552}
\end{figure}
}



\subsection{5.5.4 ~M-H 抽样方法的几个变种}
\frame{\frametitle{1. Metropolis 抽样方法}
\begin{itemize}

\item 根据提议分布~$g$ 的不同选择，M-H 抽样方法衍生出了几个不同的变种：Metropolis
抽样方法，随机游动~Metropolis 抽样方法，独立抽样方法和逐分量~M-H
抽样方法等. 下面将分别介绍这些方法.\medskip\pause

\item {\bf Metropolis 抽样方法} ~在~Metropolis 抽样方法中, 提议分布是对称的.
即~$g(\cdot|X_n)$ 满足\vspace{-2mm}
$$
g(X|Y)=g(Y|X),\vspace{-2mm}
$$
因此接受概率为\vspace{-2mm}
$$\alpha(X_t,Y)=\min\left\{1,\frac{f(Y)}{f(X_t)}\right\}.\vspace{-2mm}
$$
\end{itemize}
}

\frame{\frametitle{ 2. 随机游动~Metropolis方法}
\begin{itemize}

\item {\bf 随机游动~Metropolis方法} ~随机游动~Metropolis
抽样方法是~Metropolis 方法的一个特例. 假设候选点~$Y$
从一个对称的提议分布~$g(Y|X_t)=g(|X_t-Y|)$ 中产生的.
则在每一次迭代中， 从~$g(\cdot)$ 中产生一个随机增量~$Z,$
然后~$Y=X_t+Z.$ 比如随机增量~$Z$ 可以从均值为~$0$ 的正态分布中产生，
此时候选点~$Y|X_t\sim N(X_t,\sigma^2), ~\!\sigma^2>0.$\medskip\pause

\item 随机游动~Metropolis 算法下得到的链，
其收敛性常常对刻度参数的选择比较敏感. 当增量的方差太大时，
大部分的候选点会被拒绝，此时算法的效率很低. 如果增量的方差太小，
则候选点就几乎都被接受，因此此时随机游动~Metropolis 算法
下得到的链就几乎是随机游动了，效率也较低.  Robert et al. (1996)
提出一种选择刻度参数的方法是监视接受率，拒绝率应该在区间~[0.15,0.5]
之内才可以保证得到的链有较好的性质.
\end{itemize}
}

\frame{\frametitle{随机游动-Metropolis 算法的例子}
\begin{itemize}

\item {\bf 例~5.5.3} ~(随机游动~Metropolis)
~~使用提议分布~$N(X_t,\sigma^2)$ 和随机游动~Metropolis 算法
产生自由度为~$\nu=4$ 的~$t$ 分布随机数, 并 对方差~$\sigma^2$
的几个个不同值重复此过程，从中并选出 ~$\sigma^2$
的一个最佳值，使链的收敛性最好. 丢弃链的前~500 个迭代~(预烧期) 值,
比较~4 个生成链 观测值的十分位数和自由度为~$\nu$ 的~$t$
分布理论上的十分位数
的拟合状况，并找出与理论分布十分位数拟合最好的那个链.\medskip\pause


\item $t_\nu$ 的密度正比于~$(1+x^2/\nu)^{-(\nu+1)/2},$ 因此
$$\alpha(x_t,y)=\min\left\{1,\frac{f(y)}{f(x_t)}\right\}
=\min\left\{1,\frac{(1+y^2/\nu)^{-(\nu+1)/2}}{(1+x_t^2/\nu)^{-(\nu+1)/2}}\right\}.$$
下面我们仍然使用~$dt$ 来计算~$t$ 密度在给定点处的值.
\end{itemize}
}

\frame{\frametitle{续例~5.5.3}
\begin{itemize}

\item 上述四种方差的选择, 4 个拒绝率为~$ 0.0035, ~0.1000,
~0.4390, $ $~0.9010~.$ 只有第三个链的拒绝率在区间~[0.15,0.5] 之间.
我们可以在不同的提议分布方差下, 检查所得链的收敛性.\medskip\pause


\item 由图~5.5.3 可以看出: $\sigma=0.05$ 时, 增量太小,
几乎每个候选点都被接受了, 链在~2000 次迭代后还没有收敛.
$\sigma^2=0.5$ 时, 链的收敛较慢. $\sigma=2$ 时, 链很快收敛.
而当~$\sigma=16$ 时, 接受的概率太小, 链虽然收敛了, 但是效率很低.
\end{itemize}
}

\frame{\frametitle{续例~5.5.3}
\begin{figure}[h,p,t]\vspace{-2mm}
\centering{
\includegraphics[scale=0.52,angle=0]{fig553.pdf}}\\ \vspace{2mm}
\centerline %\caption
{\bf\it 图~\!5.5.3 ~由不同方差的提议分布生成随机游动~M-H
链}\label{fig553}
\end{figure}
}

\frame{\frametitle{续例~5.5.3}
\begin{itemize}

\item 一般说来在~MCMC 问题中，我们不能获得目标分布的理论分位数用于比较的目的.\medskip\pause

\item 但在本例中理论分位数可以通过从具有一定自由度的~$t$ 分布中直接获得,
故可用于比较.\medskip\pause

\item 忽略预烧期的~500 个值，通过应用函数获得模拟分位数的值.
目标分布的理论分位数 和~4 个链~$rw1,~rw2,~rw3$ 和~$rw4$
的样本分位数见表~5.5.1~\!.由表可见~$rw3$
的样本分位数与理论分位数比较接近.
\end{itemize}
}

\frame{\frametitle{续例~5.5.3}
\begin{table}[ht]
\centering %\caption
{\bf\it 表~5.5.1 ~目标分布的理论分位数和~4 个链的样本分位
}\label{tab6.5.1}\vspace{1mm}
\begin{tabular}{cccccc}
\hline
 {~~}&{$Q$}~~&$rw1$~~&$rw2$~~&$rw3$~~&$rw4$\\
\hline

5\% ~~& -2.132~~& 23.589~~& -1.474~~& -2.053~~& -2.929\\

10\%~~& -1.533~~& 23.707~~& -1.108~~& -1.525~~& -2.045\\

20\%~~& -0.941~~& 23.939~~& -0.402~~& -1.077~~& -0.927\\

30\%~~& -0.569~~& 24.105~~&  0.043~~& -0.625~~& -0.500\\

40\%~~& -0.271~~& 24.211~~&  0.513~~& -0.333~~& -0.187\\

50\%~~&  0.000~~& 24.459~~&  0.936~~& -0.045~~& -0.025\\

60\%~~&  0.271~~& 24.704~~&  1.766~~&  0.261~~&  0.300\\

70\%~~&  0.569~~& 24.877~~&  4.270~~&  0.533~~&  0.622\\

80\%~~&  0.941~~& 24.992~~& 13.652~~&  0.936~~&  1.051\\

90\%~~&  1.533~~& 25.129~~& 16.868~~&  1.482~~&  1.892\\

95\%~~&  2.132~~& 25.210~~& 19.250~~&  1.860~~&  2.495\\
\hline
\end{tabular}
\end{table}
}



\frame{\frametitle{一个简单的投资模型}

\noindent {\bf 例~5.5.4 }~(贝叶斯推断：一个简单的投资模型) ~一般，
\!不同的投资所得的回报是不独立的. 为了减少风险，
因此使用投资组合以保证有价证券的回报是负相关的.
这里不讨论回报的相关性， 而是对每天组合里的每个证券的收益进行排序.
假设有~5 种股票被跟踪记录了~250 交易日每天的表现, 在每一个交易日，
收益最大的股票被标记出来. ~\!用~$X_i$ 表示股票~$i$ 在~250
个交易日中胜出天数，则记录到得的频数~$(x_1,\cdots,x_5)$
为~r.v.~~\!$(X_1,\cdots,X_5)$\!~ 的观测值. 基于历史数据，假设这~5
种股票在任何给定的一个交易日能胜出的先验机会比率为\vspace{-2mm}
$$1:(1-\beta):(1-2\beta):2\beta:\beta,\vspace{-2mm}
$$
这里~$\beta\in(0,0.5)$ 是一个未知的参数.  在有了当前这~250
个交易日的数据后，使用贝叶斯方法对此比例进行更新.
}

\frame{\frametitle{续例~5.5.4}
\begin{itemize}

\item {\bf 解} ~由前所述, $\bm{X}=(X_1,\cdots,X_5)$ 在给定~$\beta$ 的条件下服从多项分布,
概率向量为\vspace{-2mm}
$$
\bm{p}=\bm{p}(\beta)=\left(\frac13,\frac{1-\beta}{3},\frac{1-2\beta}{3},\frac{2\beta}{3},\frac{\beta}{3}\right).\vspace{-3mm}
$$\pause

\item 令~$\beta$ 的先验分布~$\pi(\beta)$ 为~(0,~\!0.5) 上的均匀分布，$\bm{p}$ 的似然函数
~(即样本$\bm{x}$的分布) 为~$l(\bm{p}|\bm{x}),$
因此后验分布为\vspace{-3mm}
$$
\pi(\beta|x_1,\cdots,x_5)\!\propto\!
l\big(\bm{p}(\beta)|\bm{x}\big)\cdot\pi(\beta)\!\propto\!\frac{250!}{x_1!\cdots
x_5!}p_1^{x_1}p_2^{x_2}p_3^{x_3}p_4^{x_4}p_5^{x_5}.\vspace{-4mm}
$$\pause

\item 我们不能直接从此后验分布中产生随机数.
一种估计~$\beta$ 的方法是用随机游动~Metropolis 算法产生一个马氏链，
使其平稳分布为此后验分布，提议分布为对称的均匀分布，然后从此链中产生目标分布的随机数来估计~$\beta.$
\end{itemize}
}

\frame{\frametitle{续例~5.5.4}
\begin{itemize}

\item 此时，
接受的概率为\vspace{-2mm}
$$
\alpha(X_t,Y)=\min\left\{1,\frac{f(Y)}{f(X_t)}\right\}.\vspace{-2mm}
$$
其中~$f(\cdot)=\pi(\cdot|x_1,\cdots,x_5)$
为目标分布，且\vspace{-2mm}
$$
\frac{f(Y)}{f(X_t)}\!=\!\frac{(1/3)^{x_1}((1\!-\!Y)/3)^{x_2}((1\!-\!2Y)/3)^{x_3}((2Y)/3)^{x_4}(Y/3)^{x_5}}
{(1/3)^{x_1}((1\!-\!X_t)/3)^{x_2}((1\!-\!2X_t)/3)^{x_3}((2X_t)/3)^{x_4}(X_t/3)^{x_5}}.\vspace{-2mm}
$$
\pause

%$$
%\frac{f(Y)}{f(X_t)}=\frac{(1/3)^{x_1}((1-Y)/3)^{x_2}((1-2Y)/3)^{x_3}((2Y)/3)^{x_4}(Y/3)^{x_5}}
%{(1/3)^{x_1}((1-X_t)/3)^{x_2}((1-2X_t)/3)^{x_3}((2X_t)/3)^{x_4}(X_t/3)^{x_5}}.
%$$\pause

%\item 我们首先生成观测数据:\medskip\pause

\item 下列的频数分别是~5 种股票模拟胜出的天数:  $84,~64,~48,$ $~39,
~15~.$ 基于当年观测数据的胜出分布，我们希望估计$\beta.$
\medskip\pause

\item 下面使用随机游动~Metropolis 算法生产随机数. 这里需要两个均匀分布的随机变量，其中一个对称均匀分布用于产生提议分布，而另一个均匀分布用于
决定接受还是拒绝候选点. %下列的函数“prob”用于计算目标分布密度函数
%(忽略了密度函数的正则化常数).


\end{itemize}
}

\frame{\frametitle{图6.5.4}
\begin{figure}[h,p,t]\vspace{-7mm}
\centering{
\includegraphics[scale=0.52,angle=0]{fig554.pdf}}\\ \vspace{-1mm}
\centerline %\caption
{\bf \it 图~\!5.5.4 ~随机游动~Metropolis 算法生成~$\beta$
的链和直方图}\label{fig554}\vspace{-2mm}
\end{figure}

}

\frame{\frametitle{续例~5.5.4}
\begin{itemize}

\item 图~5.5.4 的左图表明：链的路径图显示链已经近似地收敛到目标分布.
去掉链的预烧期，生成的链可以用于估计~$\beta$. 由图~5.5.4
的右边的直方图似乎可以看出~$\beta$ 的均值接近 0.2~.\medskip\pause

\item 由链的样本产生的~5~种股票的胜出的天数、胜出频率、以及~MCMC方 法估计多项分布先验概率比的列表如下：
\smallskip\pause

胜出天数：~84, ~64, ~48, ~39, ~15\smallskip\pause

胜出频率：~0.336, ~0.256, ~0.192, ~0.156, ~0.060\smallskip\pause

先验机会比：~0.333, ~0.267, ~0.200, ~0.133, ~0.067\smallskip\pause

生成链的均值：~0.2164363\smallskip\pause

生成链的标准差：~0.02319191\smallskip\pause
\end{itemize}
}

\frame{\frametitle{3. ~独立抽样方法}
\begin{itemize}

\item Metropolis-Hastings 抽样方法的另一个特殊情形是独立抽样.\medskip\pause

\item 独立抽样中的提议分布不依赖于链的前一步状态值. 因此
~$g(Y|X_n)=g(Y),$ 接受概率为\vspace{-3mm}
$$
\alpha(X_n,Y)=\min\left\{1,\frac{f(Y)g(X_n)}{f(X_n)g(Y)}\right\}.\vspace{-3mm}
$$\pause

\item 独立抽样方法容易实施，
而且在提议分布和目标分布很接近时也趋于表现很好，
但是当提议分布和目标分布差别很大时，其表现就较差.\medskip\pause

\item Robert (1996) 讨论了独立抽样的收敛性，
并且说道:“独立抽样算法作为单独的算法很少是有用的”.
但是不管怎么样，我们仍然用下例来说明这种方法的应用.
因为独立抽样方法在混合的~MCMC 方法中是比较有用的.
\end{itemize}
}

\frame{\frametitle{独立抽样方法的例子}
\begin{itemize}

\item  {\bf 例~5.5.5}~(独立抽样) ~假设从一个正态混合分布\vspace{-2mm}
$$
pN(\mu_1,\sigma_1^2)+(1-p)N(\mu_2,\sigma_2^2)\vspace{-2mm}
$$
中观测到一个样本~$\bm{z}=(z_1,\cdots,z_n),$ 求~$p$
的估计.\medskip\pause

\item {\bf 解} ~显然, 混合正态的密度为\vspace{-2mm}
$$
f(z|p)=pf_1(z|p)+(1-p)f_2(z|p),\vspace{-2mm}
$$
其中~$f_1,f_2$ 分别为两个正态的密度.\medskip\pause

\item 设~$p$ 的先验分布~$\pi(p)$ 为
~$(0,1)$ 上的均匀分布~$U(0,1),$ 则~$p$ 的后验分布\vspace{-3mm}
$$\pi(p|\bm z)\propto f(\bm z|p)\pi(p)=\prod\limits_{j=1}^n\big[ pf_1(z _j|p)+(1-p)f_2(z_j|p)\big]\vspace{-2mm}
$$

\item 提议分布的支撑应和~$p$ 的取值范围~$(0,1)$ 相同.
%$p$ 的先验分布
%最明显的选择就是贝塔分布.
在没有先验信息的情况下，这里使用贝塔分布~$Be(1,1)$ 作为提议分布.
%($Be(1,1)$就是均匀分布$U(0.1)$).

\end{itemize}
}

\frame{\frametitle{续例~5.5.5}
\begin{itemize}

\item 候选点~$Y$ 被接受的概率为\vspace{-3mm}
$$
\alpha(X_n,Y)=\min\left\{1,\frac{\pi(Y|\bm{z})g(X_n)}{\pi(X_n|\bm{z})g(Y)}\right\}.\vspace{-2mm}
$$
其中~$g$ 提议分布, 其密度函数~$g(y)\propto
y^{a-1}(1-y)^{b-1},~\pi(\cdot|\bm{z})$ 为目标分布. 此处\vspace{-3mm}
$$
\frac{\pi(y|\bm{z})g(x_t)}{\pi(x_t|\bm{z})g(y)}\!=\!\frac{x_t^{a-1}(1-x_t)^{b-1}\prod\limits_{j=1}^n[yf_1(z_j|y)\!+\!(1-y)f_2(z_j|y)]}
{y^{a-1}(1-y)^{b-1}\prod\limits_{j=1}^n[x_tf_1(z_j
|x_t)\!+\!(1-x_t)f_2(z_j|x_t)]}.\vspace{-3mm}
$$\pause

%$$
%\frac{\pi(y|\bm{z})g(x_n)}{\pi(x_n|\bm{z})g(y)}\!=\!\frac{x_n^{a-1}(1-x_n)^{b-1}\prod_{j=1}^n[yf_1(z_j|y)\!+\!(1-y)f_2(z_j|y)]}
%{y^{a-1}(1-y)^{b-1}\prod_{j=1}^n[x_nf_1(z_j
%|x_n)\!+\!(1-x_n)f_2(z_j|x_n)]}.\vspace{-2mm}
%$$\pause

\item 下面我们进行模拟, 提议分布取为~$U(0,1)=Be(1,1)$. 观测数据从下述正态混合中产生\vspace{-2mm}
$$0.2N(0,1)+0.8N(5,1)\vspace{-5mm}
$$
\end{itemize}
}

\frame{\frametitle{图~5.5.5}

 链的路径图和丢掉~100~个预烧期样本后的直方图如下:\vspace{-8mm}


\begin{figure}[h,p,t]\vspace{-7mm}
\centering{
\includegraphics[scale=0.52,angle=0]{fig555.pdf}}\\ \vspace{-1mm}
%\caption
\centerline{\bf\it 图~\!5.5.5 ~提议分布为~$Be(1,1)$ 时独立抽样生成
的链和直方图}\label{fig555}\vspace{-3mm}
\end{figure}

}

\frame{\frametitle{续例~5.5.5}
\begin{itemize}

\item 丢掉~100 个预烧期样本后，由图~5.5.5 左图可见
链的时间状态图显示链混合的很好, 很快收敛到平稳分布. 图~5.5.5
右图是生成链的直方图，由图可见保留样本的样本均值是~0.2553,
它就是~$p$ 的估计值. %链的路径图和直方图的代码如下：
\end{itemize}
}

\frame{\frametitle{4. ~逐分量的~M-H 抽样方法}
\begin{itemize}

\item 当状态空间为多维时, 不整体更新~$\mathbf{X}_n,$
而是对其分量进行逐个更新, 即称为逐分量的~M-H 抽样方法.
这样做更方便和更有效率.\medskip\pause

\item 记\vspace{-3mm}
\begin{eqnarray*}
&& \mathbf{X}_n=(X_{n,1},\ldots,X_{n,k}),\\
&&
\mathbf{X}_{n,-i}=(X_{n,1},\ldots,X_{n,i-1},X_{n,i+1},\ldots,X_{n,k}).
\end{eqnarray*}\\\vspace{-2mm}
分别表示在第~$n$ 步链的状态， 及在第~$n$ 步除第~$i$
个分量外其他分量的状态. $f(\mathbf{x})=f(x_1,\cdots,x_k)$
为目标分布，\vspace{-2mm}
$$
f(x_i|\mathbf{x}_{-i})=f(\bm x)\Big/\int
f(x_1,\cdots,x_k)dx_i\vspace{-2mm}
$$
表示~$X_i$ 对其他分量的条件密度.
\end{itemize}
}

\frame{\frametitle{逐分量的~M-H 抽样方法}
\begin{itemize}

\item 逐分量的~M-H 抽样方法由~$k$ 步构成：
令~$X_{n,i}$ 表示在第~$n$ 次迭代后~$\mathbf{X}_n$ 第~$i$
个分量的状态， 则在第~$n+1$ 步迭代的第~$i$ 步中， 使用~M-H
算法更新~$X_{n,i}.$ 做法如下：\medskip\pause

\item 对~$i=1,\cdots,k,$
从第~$i$ 个提议分布~$q_i(\cdot|X_{n,i},\mathbf{X}^*_{n,-i})$
中产生~$Y_i,$ 这里\vspace{-3mm}
$$\mathbf{X}^*_{n,-i}=(X_{n+1,1},\ldots,X_{n+1,i-1},X_{n,i+1},\ldots,X_{n,k}).\vspace{-4mm}
$$\pause

\item 然后以概率\vspace{-2mm}
$$
\alpha(\bm{X}^*_{n,-i},X_{n,i},Y_i)=\min\left\{1,\frac{f(Y_i|\mathbf{X}^*_{n,-i})
q_i(X_{n,i}|Y_i,\mathbf{X}^*_{n,-i})}
{f(X_{n,i}|\mathbf{X}^*_{n,-i})q_i(Y_i|X_{n,i},\mathbf{X}^*_{n,-i})}\right\}\vspace{-2mm}
$$
若~$Y_i$ 被接受， 则令~$X_{n+1,i}=Y_i$；
否则令~$X_{n+1,i}=X_{n,i}~\!.$
\end{itemize}
}

\frame{\frametitle{逐分量的~M-H 抽样方法的例子}
\begin{itemize}

\item {\bf 例~5.5.6} ~考虑~54 位老年人的智力测试成绩~(Wechsler Adult Intelligence Scale, WAIS, 0-20分).
研究的兴趣在于发现老年痴呆症与~WAIS 成绩的关系.\medskip\pause

\item {\bf 解} ~我们采用如下简单的~logistic 回归模型:\vspace{-2mm}
$$
Y_i\sim B(1,\pi_i),\quad
\ln\frac{\pi_i}{1-\pi_i}=\beta_0+x_i\beta_1,~i=1,\ldots,54.\vspace{-2mm}
$$
其中~$Y_i=1$ 表示第~i 个人老年痴呆，$x_i$ 表示第~i 个人的~WAIS
成绩.\medskip\pause

\item 似然函数为\vspace{-2mm}
\begin{eqnarray*}
f({\bf y}|\beta_0,\beta_1)\!\!\!\!&&=
\prod_{i=1}^n\left(\frac{e^{\beta_0+x_i\beta_1}}{1+e^{\beta_0+x_i\beta_1}}\right)^{y_i}
\left(\frac{1}{1+e^{\beta_0+x_i\beta_1}}\right)^{1-y_i}\\
&&= \exp\bigg\{n\bar{y}\beta_0+\beta_1\sum_{i=1}^n
x_iy_i-\sum_{i=1}^n \ln\Big(1+e^{\beta_0+x_i\beta_1}\Big)\bigg\}.
\end{eqnarray*}
\end{itemize}
}

\frame{\frametitle{续例~5.5.6}
\begin{itemize}

\item 考虑~$\beta_0,~\beta_1$ 的先验分布~$\pi$ 为独立的
正态分布, 即~$\pi(\beta_0,\beta_1)=$
$\pi_1(\beta_0)~\!\pi_2(\beta_1),$ 而\vspace{-2mm}
$$\beta_j\sim N(\mu_{\D \beta_j},\sigma^2_j), \quad j=0,1.\vspace{-2mm}
$$
其中当~$\mu_{\D \beta_j}=0,~\sigma^2_j$ 很大时,
上述先验分布接近无信息先验.\medskip\pause

\item 因此后验分布为\vspace{-2mm}
\begin{eqnarray*}
\pi(\beta_0,\beta_1|{\bf y})\!\!\!\!&&\propto f({\bf y}|\beta_0,\beta_1)\pi(\beta_0,\beta_1)\\
&&\propto\exp\bigg\{\sum_{i=1}^n\Big[\big(\beta_0+\beta_1x_i\big)y_i-\ln\big(1+e^{\beta_0+x_i\beta_1}\big)\Big]\\
&&~~~-\frac{(\beta_0-\mu_{\beta_0})^2}{2\sigma^2_0}-\frac{(\beta_1-\mu_{\beta_1})^2}{2\sigma^2_1}\bigg\}.
\end{eqnarray*}
\end{itemize}
}

\frame{\frametitle{续例~5.5.6}
\begin{itemize}

\item 下面我们用逐分量的~M-H 算法从此分布中产生随机数.
提议分布取为\vspace{-2mm}
\begin{eqnarray*}
\bm{\beta}'\!\sim\!
N_2\big(\bm\beta,diag\{\bar{s}^2_{\beta_0},\bar{s}^2_{\beta_1}\}\big)\Leftrightarrow
\beta'_0 \!\sim\! N(\beta_0,\bar{s}^2_{\beta_0}),~\beta'_1 \!\sim\!
N(\beta_1,\bar{s}^2_{\beta_1}).
\end{eqnarray*}\\\vspace{-2mm}
此处~$\bm{\beta}'=(\beta'_0,\beta'_1)^\tau,$
而~$\bm\beta=(\beta_0,\beta_1)^\tau.$\medskip\pause

\item 在~M-H 算法中, 按分量进行逐个更新, 其
优势在于应用方便, 不需要考虑调节参数. 算法如下:\\\medskip
对~$t=1,\cdots,T$:\medskip\pause

\item (1) ~令~$\bm\beta=\big(\beta_0^{(t-1)},\beta_1^{(t-1)}\big)^\tau$.\medskip\pause

\item (2)
~从提议分布~$N(\beta_0,\bar{s}^2_{\beta_0})$ 产生候选点~$\beta'_0.$
\end{itemize}
}

\frame{\frametitle{续例~5.5.6}
\vspace{-2mm}

\begin{itemize}

\item (3) ~令~$\bm\beta'=\big(\beta'_0,\beta_1^{(t-1)}\big)^\tau$\!,
    接受概率~$\alpha_0(\bm\beta,\bm{\beta}')=\min\left\{1,A\right\},$
    其中\vspace{-1mm}
    $$ A= \frac{\pi(\beta_0',\beta_1|\bm{y})}
    {\pi(\beta_0,\beta_1|\bm{y})}=\frac{f(\textbf{y}|\beta'_0,\beta_1)\pi(\beta'_0,\beta_1)}
    {f(\textbf{y}|\beta_0,\beta_1)\pi(\beta_0,\beta_1)}.\vspace{-3mm}
    $$\pause

\item (4) ~以概率~$\alpha_0(\bm\beta,\bm{\beta}')
$ 接受~$\bm\beta=\bm{\beta}',$ 即将~$\bm\beta$
的第一个分量~$\beta_0$ 更新为~$\beta'_0;$ 否则~$\bm{\beta}$
保持不变.\smallskip\pause

\item (5) ~从提议分布~$N(\beta_1,\bar{s}^2_{\beta_1})$ 产生候选点~$\beta'_1$.\smallskip\pause

\item (6) ~令~\!$\bm{\beta}'=(\beta_0,\beta'_1)^\tau$\!, 接受概率
~\!$\alpha_1(\bm\beta,\bm{\beta}')=\min\left\{1,B\right\},$\!
其中\vspace{-2mm}
$$
B=\frac{\pi(\beta_0,\beta'_1|\bm{y})}
    {\pi(\beta_0,\beta_1|\bm{y})}=\frac{f(\textbf{y}|\beta_0,\beta'_1)\pi(\beta_0,\beta'_1)}
{f(\textbf{y}|\beta_0,\beta_1)\pi(\beta_0,\beta_1)}.\vspace{-4mm}
$$\pause

\item (7)
~以概率~$\alpha_1(\bm\beta,\bm{\beta}')$
接受~$\bm\beta=\bm{\beta}',$ 即将~$\bm\beta$ 的 第二个分量~$\beta_1$
更新为~$\beta'_1;$ 否则~$\bm\beta$ 保持不变.\smallskip\pause

\item (8) ~令~$\bm{\beta}^{(t)}=\bm{\beta}$.
\end{itemize}
}

\frame{\frametitle{图~5.5.6}\vspace{-2mm}

\begin{figure}[h,p,t]\vspace{-7mm}
\centering{
\includegraphics[scale=0.52,angle=0]{fig556.pdf}}\\ \vspace{2mm}
%\caption
\centerline{\bf\it 图~\!5.5.6 ~逐分量~M-H 算法~logistic
回归参数~$\beta_0$ 和~$\beta_1$ 的诊断图}\label{fig556}
\end{figure}\vspace{-2mm}
}

\frame{\frametitle{续例~5.5.6}
\begin{itemize}

\item 逐分量随机游动~M-H
算法~logistic 回归参数的诊断图，以及去掉预烧期后由链获得的
~$\bm\beta$ 的两个分量~$\beta_0$ 和~$\beta_1$
的样本均值和样本标准差如下：\medskip\pause

\item 图~5.5.6 最上面两个小图，即~$\beta_0$ 和~$\beta_1$ 的时间状态图显示链的混合较好.
中间两个小图， 即~$\beta_0$ 和~$\beta_1$ 的遍历均值图，去掉前~2000
次迭代~(即预烧期) 后趋于平稳，也表明链的收敛性较好.
最下面两个小图，即链的自相关函数图中显示马氏链的抽样步长~$L=20$
次迭代，也就是每间隔~20
个样本抽取一个，获得的样本的自相关性就很低了，这几种诊断方法都显示链的收敛性较好.
最后，$\beta_0$ 和~$\beta_1$
的链的样本均值分别是~$2.8793214,~-0.3737996.$
~样本标准差分别是~$1.3479687,$ $0.1287645.$
\end{itemize}
}

\subsection{5.5.5 ~Gibbs 抽样方法}

\frame{\frametitle{引言}
\begin{itemize}

\item Gibbs 抽样最早是由~Geman 和~Geman (1984) 提出
并被用于~Gibbs 格子点分布，因而由此得名的. 它是~M-H
抽样的另一种特殊情形.\medskip\pause

\item Gibbs 抽样方法特别适合于目标分布是多元的场合，其最令人感兴趣的方面是为了产生不可约的、
非周期、并以高维空间的目标分布作为其平稳分布的马氏链，
只需要从一些一元分布中进行抽样就可以了.
即将从多元目标分布中抽样转化为从一元目标分布抽样，这是~Gibbs
抽样的重要性所在.
\end{itemize}
}

\frame{\frametitle{ ~Gibbs 抽样方法}
\begin{itemize}

\item 令~$\bm X=(X_1,\cdots,X_k)$ 为~$\mathbb{R}^k$ 中的随机变量， 其联合分布~$f(\bm
x)$ $=f(x_1,\cdots,x_k)$ 为目标抽样分布.  定义~$k-1$
维的随机变量\vspace{-2mm}
$$
\mathbf{X}_{-j}=(X_1,\cdots,X_{j-1},X_{j+1},\cdots,X_k),\vspace{-5mm}
$$\pause

\item 记~$X_j|\mathbf{X}_{-j}$ 的{\bf 满条件分布}密度为
~$f\big(x_j~\!|~\!\mathbf{x}_{-j}\big),~j=1,\cdots,k.$ 则~Gibbs
抽样方法是从这~$k$ 个条件分布中产生候选点，以解决直接从~$f$
中进行抽样的困难. 算法如下:\medskip\pause

\item (1) ~在~$t=0$ 时, 初始化~$\mathbf{X}(0)$;

\end{itemize}
}

\frame{\frametitle{Gibbs 抽样方法}

(2) ~对~$t=1,2,\cdots,J$ 重复下列步骤\medskip\pause

~~~~~(a) 令~$(x_1,x_2,\cdots,x_k)=\mathbf{X}(t-1)$.\medskip\pause

~~~~~(b) 对每个分量~$j=1,\ldots,k$,\medskip\pause

\indent\indent\indent (i) 从~$f(X_j|\mathbf{x}_{-j})$
中产生候选点~$X_j^*(t)$.\medskip\pause

\indent\indent\indent (ii) 更新~$x_j=X_j^*(t)$.\medskip\pause

~~~~~(c) 令~$\mathbf{X}(t)=(X_1^*(t),\ldots,X_k^*(t))$
(每个候选点都被接受)\medskip\pause

~~~~~(d) 增加~$t$ 返回到~(a).


}

\frame{\frametitle{Gibbs 抽样方法}
\begin{itemize}

\item 注意在上述算法~(b) 步抽样中，各个分量依次被更新：\vspace{-2mm}
\begin{eqnarray*}
&& x_1(t)\sim f\big(x_1~\!|~\!x_2(t-1),x_3(t-1),\cdots,x_k(t-1)\big);\\
&& x_2(t)\sim f\big(x_2~\!|~\!x_1(t),x_3(t-1),\cdots,x_k(t-1)\big)\\
&&~~\vdots\\
&& x_k(t)\sim f\big(x_k~\!|~\!x_1(t),x_2(t),\cdots,x_{k-1}(t)\big).
\end{eqnarray*}\\\vspace{-2mm}\pause

\item 从一元分布~\vspace{-2mm}
$$f\big(x_j~\!|~\!x_1(t),x_2(t),\cdots,x_{j-1}(t),x_{j+1}(t-1),\cdots,x_k(t-1)\big)\vspace{-2mm}
$$
中抽样是比较容易的，
因为~$f\big(x_j~\!|~\!\mathbf{x}_{-j}\big)\propto f(\mathbf{x}),$
其中除了变量~$x_j$ 外, 其他变量都是常数.\medskip\pause
\end{itemize}
}

\frame{\frametitle{Gibbs 抽样方法}
\begin{itemize}

\item Gibbs 抽样方法的合理性并不需要通过验证其是~M-H 抽样
方法的一种特例来证明. Gibbs
抽样方法的一个特别之处就是满条件分布可以唯一确定联合分布.
这就是著名的~{\bf Hammersley-Clifford 定理}，详见韦和张~(2013)
定理~5.5.2.\medskip\pause

\item 前面所述的~Gibbs 抽样方法可直接在贝叶斯框架下实现.
令~$\bm{\theta}=(\theta_1,\cdots,\theta_k),$ 而\vspace{-3mm}
$$\bm{\theta}_{-j}=(\theta_1,\cdots,\theta_{j-1},\theta_{j+1},\cdots,\theta_k).\vspace{-5mm}
$$\pause

\item 将~$\bm {x}$ 用感兴趣的参数~$\bm\theta$
代替，将目标分布~$f(\mathbf{x})$ 用
后验分布~$\pi(\bm\theta|\mathbf{x})$ 代替, 将满条件分布
~$f\big(x_j~\!|~\!\mathbf{x}_{-j}\big),~j=1,\cdots,k$
用~$\pi\big(\theta_j~\!|~\!\bm{\theta}_{-j},\mathbf{x}\big),~j=1,\cdots,k$
代替.%\medskip
\end{itemize}
}


\frame{\frametitle{Gibbs 抽样方法}
\begin{itemize}

\item 因此，在贝叶斯分析中，Gibbs 抽样方法产生马氏链~$\big\{\bm{\theta}^{(i)},$ $i=1,2,\cdots\big\}
$ 的方法概括如下：\medskip\pause

\item (1) ~选择参数的一个初始值~$\bm{\theta}^{(0)};$\medskip\pause

\item (2) ~对~$t=1,2,\cdots,J$ 重复下列步骤\medskip\pause

~~~~~(a)
~令~$(\theta_1,\theta_2,\cdots,\theta_k)=\bm{\theta}^{(t-1)}$\smallskip\pause

~~~~~(b) ~对每个分量~$j=1,\ldots,k$\smallskip\pause

~~~~~~~~~~~~\indent\indent\indent (i)
从~$\pi\big(\theta_j|\bm{\theta}_{-j},\mathbf{x}\big)$
中产生候选点~$\theta_j^{(t)}$.\smallskip\pause

~~~~~~~~~~~~\indent\indent\indent (ii)
更新~$\theta_j=\theta_j^{(t)}$.\smallskip\pause

~~~~~(c)
令~$\bm{\theta}^{(t)}=(\theta_1^{(t)},\ldots,\theta_k^{(t)}).$\smallskip\pause

~~~~~(d) 增加~$t$ 返回到~(a).\medskip\pause
\end{itemize}
}

\frame{\frametitle{Gibbs 抽样方法}
\begin{itemize}

\item 注意在上述算法~(b) 步抽样中，各个分量依次被更新：\medskip

~~~~~1)
~抽取~$\theta_1^{(t)}\sim\pi\big(\theta_1~\!\big|~\!\theta_2^{(t-1)},\theta_3^{(t-1)},
\cdots,\theta_k^{(t-1)},\mathbf{x}\big);$\medskip\pause

~~~~~2)
~抽取~$\theta_2^{(t)}\sim\pi\big(\theta_2~\!\big|~\!\theta_1^{(t)},\theta_3^{(t-1)},
\cdots,\theta_k^{(t-1)},\mathbf{x}\big);$\medskip\pause

~~~~~~$\vdots$\medskip

~~~~~k)
~抽取~$\theta_k^{(t)}\sim\pi\big(\theta_k~\!\big|~\!\theta_1^{(t)},\theta_2^{(t)},
\cdots,\theta_{k-1}^{(t)},\mathbf{x}\big).$\medskip\pause

\item 下例说明如何利用~Gibbs 抽样生成二元正态分布的样本.%\medskip\pause
\end{itemize}
}

\frame{\frametitle{Gibbs 抽样生成二元正态分布随机数的例子}
\begin{itemize}

\item {\bf 例~5.5.7} ~~使用~Gibbs 抽样产生二元正态分布\vspace{-2mm}
$$N\big(\mu_1,\mu_2,\sigma_1^2,\sigma_2^2,\rho\big)\vspace{-2mm}
$$
的随机数.\medskip\pause

\item {\bf 解} ~在二元正态场合，$X_1|X_2$ 以及~$X_2|X_1$
仍然服从正态分布，且易知\vspace{-3mm}
\begin{eqnarray*}
&&E[X_1|X_2=x_2]=\mu_1+\rho\frac{\sigma_1}{\sigma_2}(x_2-\mu_2),\\
&&Var[X_1|X_2=x_2]=(1-\rho^2)\sigma_1^2.
\end{eqnarray*}\\\vspace{-2mm}\pause

\item 类似可得~$X_2|X_1$ 的分布. 因此\vspace{-2mm}
\begin{eqnarray*}
&&f(x_1|x_2)\sim N\Big(\mu_1+\rho\frac{\sigma_1}{\sigma_2}(x_2-\mu_2),(1-\rho^2)\sigma_1^2\Big),\\
&&f(x_2|x_1)\sim
N\Big(\mu_2+\rho\frac{\sigma_2}{\sigma_1}(x_1-\mu_1),(1-\rho^2)\sigma_2^2\Big).
\end{eqnarray*}
\end{itemize}
}

\frame{\frametitle{续例~5.5.7}

从而使用~Gibbs 算法如下\medskip\pause

~~~~(1) ~令~$(x_1,x_2)=\mathbf{x}(t-1)$;\medskip\pause

~~~~(2) ~从~$f(x_1|x_2)$ 中产生候选点~$X_1^*(t)$.\medskip\pause

~~~~(3) ~更新~$x_1=X_1^*(t)$.\medskip\pause

~~~~(4) ~从~$f(x_2|x_1)$ 中产生~$X_2^*(t)$.\medskip\pause

~~~~(5) ~更新~$x_2=X_2^*(t)$.\medskip\pause

~~~~(6) ~令~$X(t)=(X_1^*(t),X_2^*(t))$.\medskip

生成链的散点图如下： }

\frame{\frametitle{图~5.5.7}

\begin{figure}[h,p,t]\vspace{-7mm}
\centering{
\includegraphics[scale=0.42,angle=0]{fig557.pdf}}\\ \vspace{0mm}
%\caption
\centerline{\bf\it 图~\!5.5.7 ~Gibbs
抽样生成的二元正态链的散点图}\label{fig557}\vspace{-3mm}
\end{figure}
}

\frame{\frametitle{续例~5.5.7}
\begin{itemize}

\item 产生的链开始的~1000 个观测值被丢弃掉，剩下的观测值存在~$x$ 中，
由此观察值获得样本均值、样本协方差矩阵和样本相关系数矩阵如下所示.
各参数的样本估计值离真值很近，散点图~5.5.7 也显示出
二元正态所具有的椭圆对称性和负相关性特征.\medskip\pause

\item 模拟样本均值为 $-0.0735532, ~2.0297997$\medskip\pause

\item 样本协方差矩阵

~~~~$~0.9609837 ~-0.3506705$

~~~~$-0.3506705 ~~0.2327129$\medskip\pause


\item 样本相关系数矩阵

~~~~$~1.0000000 ~-0.7415338$

~~~~$-0.7415338 ~~1.0000000$
\end{itemize}
}

\frame{\frametitle{Gibbs 抽样算法的另一个例子}
\begin{itemize}

\item {\bf 例~5.5.8} ~(身体温度数据) 考虑~Mackowiak et al.~(1992) 的数据，该数据记录了
~130 个人的身体温度~(华氏)，性别和每分钟的心跳.
实验的目的是检验~Carl Wunderlich 的观点--
健康成年人的体温平均为~37$^oC$~(~\!=~\!98.6$^o~\!F$).\medskip\pause

\item 记温度为~$Y_i,~i=1,\cdots,n$，并假设正态模型\vspace{-2mm}
$$Y_i\sim N(\mu,\sigma^2).\vspace{-5mm}
$$\pause

\item 令~$\mu$ 和~$\sigma^2$ 的先验独立，分别为正态和逆伽玛分布，即\vspace{-2mm}
$$\mu\sim N(\mu_0,\sigma_0^2),\quad \sigma^2\sim IG(a_0,b_0),\vspace{-5mm}
$$\pause

\item 此时我们的目标分布为~$\mu,\sigma^2$ 的后验分布\vspace{-2mm}
$$\pi(\mu,\sigma^2|\bm{y})\propto f(|\bm{y}|\mu,\sigma^2)\pi_1(\mu)\pi_2(\sigma^2).
$$
\end{itemize}
}

\frame{\frametitle{续例~5.5.8}
\begin{itemize}

\item 为使用~Gibbs 抽样算法，我们必须首先计算两个满条件分布
~$\pi_1(\mu|\sigma^2,\bm{y})$ 与~$\pi_2(\sigma^2|\mu,\bm{y}).$
经过计算得到\vspace{-2mm}
\begin{eqnarray*}
\mu\big|(\sigma^2,\bm{y})\!\!\!\!&&\sim
N\big(\omega\bar{y}+(1-\omega)\mu_0,~\omega\sigma^2/n\big), \quad
\omega=\frac{\sigma_0^2}{\sigma^2/n+\sigma^2_0}.\\
\sigma^2\big|(\mu,\bm{y})\!\!\!\!&&\sim
\Gamma^{-1}\Big(a_0+\frac{n}{2},~b_0+\frac12\sum_{i=1}^n(y_i-\mu)^2\Big).
\end{eqnarray*}\\\vspace{-2mm}\pause

\item 使用上述结果，Gibbs 抽样算法如下：\medskip\pause

对~$t=1,\cdots,J:$\medskip

\item (1) ~令~$\mu=\mu^{(t-1)},\sigma=\sigma^{(t-1)}$.%\medskip\pause
\end{itemize}
}

\frame{\frametitle{续例~5.5.8}
\begin{itemize}

\item (2) ~令~$\omega=\sigma_0^2\big/[\sigma^2/n+\sigma^2_0]$, 计算\vspace{-3mm}
$$m=\omega\bar{y}+(1-\omega)\mu_0,~~s^2=\omega~\!\frac{\sigma^2}{n}.\vspace{-5mm}
$$.
\pause

\item (3) ~从~$N(m,s^2)$ 中产生~$\mu$.\medskip\pause

\item (4) ~令~$\mu^{(t)}=\mu$.\medskip\pause

\item (5) ~计算~$a=a_0+\frac{n}{2}, b=b_0+\sum_{i=1}^n(y_i-\mu)^2\big/2$.\medskip\pause

\item (6) ~从~$G(a,b)$ 中产生~$\tau$.\medskip\pause

\item (7) ~令~$\sigma^2=1/\tau$ 以及~$\sigma^{(t)}=\sigma.$\medskip\pause

算法收敛性的诊断图如下：
\end{itemize}
}

\frame{\frametitle{图~5.5.8}

\begin{figure}[h,p,t]\vspace{-5mm}
\centering{
\includegraphics[scale=0.46,angle=0]{fig558.pdf}}\\ \vspace{0mm}
\centerline %\caption
{\bf\it 图~\!5.5.8 ~利用~Gibbs 抽样对正态模型的~MCMC
诊断图}\label{fig558}
\end{figure}
}

\frame{\frametitle{续例~5.5.8}
\begin{itemize}

\item 图~5.5.8 算法收敛性的诊断图，包括~$\mu$ 和~$\sigma$ 的轨迹图、遍历均值图和自相关函数图.
\medskip\pause

\item 这些图都显示算法已经达到收敛.
这是因为轨迹图没有不规则之处，遍历均值图很快就
达到平稳正态，自相关函数图中只有样本的第一步自相关性较高，其余都较低.
这些都表示 由链生成的样本的表现都很好.\medskip\pause

\item 模拟样本均值为~$98.2492694, ~0.7370387$\medskip\pause

\item 模拟样本标准差为~$ 0.06522787, ~0.04698446$
\end{itemize}
}


\section{5.6 ~R 与~WinBUGS 软件简介}
\subsection{5.6.1 ~WinBUGS 软件简介}
\frame{\frametitle{WinBUGS 软件简介}
\begin{itemize}

\item R~(~S 和~Splus 的开源版本) 是一套完整的数据处理、计算和制图软件系统.
功能包括：数据存储和处理系统；
数组运算工具~(其向量、矩阵运算功能尤其强大)；完整连贯的统计分析工具；
优秀的统计制图功能；简便而强大的编程语言.\medskip\pause

\item 对贝叶斯计算来说，可以直接使用~Gibbs 抽样和~Metroplis 算法
进行编程，前面都是如此做的. 对于高强度的计算问题，可以使用~Fortran
或~C 语言然后从~R 中通过接口调用.\medskip\pause

\item BUGS~(Bayesian inference Using Gibbs Sampling)
是一个免费软件，其允许用户指定复杂的多层模型并使用~MCMC
算法来估计模型中的未知参数. WinBUGS 是~BUGS 在~Windows
平台下的版本，而~openBUGS 是~BUGS 的开源版本.
还有一些其它的版本，可以访问~$http://www.mrc-bsu.cam.ac.$ $uk/bugs/$
查看详细介绍.
\end{itemize}
}

\frame{\frametitle{WinBUGS 软件简介}
\begin{itemize}
\item WinBUGS，它与使用其他软件得到的分析结果是一致的，即便是它们使用不同的MCMC算法.
使用~WinBUGS 软件进行贝叶斯分析的详细介绍也可以参考~Ntzoufras
(2009).\medskip\pause

\item WinBUGS 软件最新版本是~WinBUGS 1.4.3.
下载地址是如前所述的网站，即~http://www.mrc-bsu.cam.ac.uk/bugs/\medskip\pause

\item 双击~WinBUGS 版本的安装文件即可完成~WinBUGS 的安装.
~WinBUGS
软件图标位于安装目录下，为操作方便，建议用快捷方式将其拖至桌面.
双击该图标将其打开， 阅读自动弹出的许可证协议窗口后将其关闭.
单击左上角主菜单的~File 项，在下拉菜单中选择~New,
打开一个新的文件窗口，在此区域编写程序.

\end{itemize}
}
\subsection{5.6.2 ~使用~R 调用~WinBUGS}
\frame{\frametitle{使用~R 调用~WinBUGS}
\begin{itemize}

\item 另外一种方式. 使用~R 里的包~R2WinBUGS
来调用~WinBUGS. 安装好~WinBUGS 和此包后， 就可以在~R 里使用~WinBUGS
进行贝叶斯分析了. 在~R 里使用~WinBUGS 进行贝叶斯分析的主要函数为
~bugs, 在运行此函数后，WinBUGS
会自动打开并运行相应的过程，然后自动关闭.\medskip\pause

\item 在~R 里还可以调用包~coda 来作更多的分析. 如绘制自相关图,
诊断收敛性等. 请参考~coda 的说明文档.

\end{itemize}
}


\end{CJK*}
\end{document}

\frame{\frametitle{}
\begin{itemize}

\item

\item

\item

\end{itemize}
}

\frame{\frametitle{}
\begin{itemize}

\item

\item

\item

\end{itemize}
}
