ShihanRan's Blog Life is like a Markov Chain.

Statistical Computing


This is my course note taken in Chinese.

Monte Carlo Methods for Hypothesis Testing

  $H_0$为真 $H_1$为真
拒绝$H_0$ Type I error $\le \alpha$  
不拒绝$H_0$   Type II error

关于一型错误率的估计

  • 由于一型错误率是指在零假设正确的时候,拒绝它的概率。因此我们可以在零假设成立的条件下,重复m次检验过程,观察到的拒绝率则为一型错误率

  • 对于估计的错误率的不确定度的估计,一般都是用方差来做。

    • 在这里我们求到的错误率是一个bernoulli分布,因为每次都是判断p-value和$\alpha$的大小,然后决策结果为0和1,所以方差就应该为bernoulli的方差p(1-p)。
    • 但由于我们使用的均值估计量,所以方差的估计应该为$\frac{\hat p (1-\hat p)}{m}$
  • 在构造置信区间的时候,我们利用了$\chi^2$分布和t分布,所以我们==需要知道的一个证明==为:

    (1) $\bar X \sim N(\mu,\sigma)$ then $\frac{\sum(X_i - \bar X)^2}{\sigma^2}\sim \chi^2_{n-1}$

    (2) $\bar X$ and $\sum(X_i - \bar X)^2$ are independent, and therefore $\frac{\sqrt{n}(\bar X - \mu)}{s} = \frac{\bar X - \mu}{SE(\bar X)}\sim t_{n-1}$

Power of a Test

  • 首先关注“功效power”的定义

    • 不管$\theta$是在哪个参数空间,始终是拒绝的$H_0$

  • 较低的二型错误率对应较高的功效,因为二型错误率 ($H_1$为真,但却没有拒绝$H_0$)为$1-\pi(\theta_1)$

  • 由于一般$\Theta_1$都是一个比较大的空间,而不是像$\Theta_0$一样一般是一个点,所以在用蒙特卡诺的时候,记得每次先选择一个特定的$\Theta_1$,再在这个对立假设空间中开始蒙特卡诺。

  • 注意:$H_1$假设下,T服从的是非中心的t分布,而不是对称的t分布。

Power Comparisons

  • 可以用Monte Carlo来比较不同检验的功效,以此找到错误率最低的检验方法。

The Bootstrap

一些理解

  • 非参数Monte Carlo方法,总体分布完全未知,样本是唯一已有的信息,将观测到的样本视为一个有限总体,从中进行随机抽样来估计总体的特征以及对抽样总体作出统计推断。

  • (当目标总体有指定的参数信息,比如正态中已知均值和方差则总体分布完全已知了,此时我们可以利用Monte Carlo方法从总体中抽样。

  • 以样本为总体进行的抽样,可以据此估计统计量的分布,而统计量的一些性质,比如偏差和标准差等等也可以通过再抽样来进行估计。

  • 关系可以表述为:

    • $F \rightarrow X \rightarrow F_n$, 真实的分布生成了我们观测到的X样本,我们根据这个样本估计出一个经验分布函数$F_n$。

    • $F_n \rightarrow X^* \rightarrow F_n^* $,我们利用Bootstrap方法,从X中再抽样,等价于从$F_n$中产生随机样本$X^*$。

    • 需要注意的是,如果$F_n$并不能逼近真实分布F,则重复抽样下的分布$F_n^*$也不会靠近F。

Bootstrap Estimation of Standard Error

  • B大概到50就足够了。
  • 在R中有现成的package可以调用:bootbootstrap,只需要传入数据和统计量函数

Bootstrap Estimation of Bias

  • 对于估计bias,需要更大的B,大概2000.

Bootstrap Confidence Intervals for $\theta$

The Standard Normal Bootstrap Confidence Interval

  • 也即统计学I学的Wald test。
  • 若$\hat\theta$为$\theta$的无偏估计, 那么$\theta$的一个渐近的95%标准正态 Bootstrap 置信区间为$\hat\theta\pm 1.96\hat{se}_B(\hat\theta)$。
  • 缺点:假设多(要为无偏估计,要有$\theta$分布是正态的假设或CLT成立)
  • 优点:好计算。

The Percentile Bootstrap Confidence Interval

  • 为Bootstrap重复样本$\hat \theta^*$的百分位数
  • 与标准正态置信区间的差异主要在于样本相关系数的分布是不是靠近正态分布,当很靠近正态的时候,两个区间就会趋向一致。

The Basic Bootstrap Confidence Interval

  • 与百分位数置信区间的差异在于,需要多算一个$\hat\theta$

The Bootstrap t Interval

  • 步骤是在bootstrap里嵌套bootstrap:

  • 并没有使用t分布作为推断分布,而是使用“再抽样方法”得到一个“t类型”的统计量的分布。

The Jackknife

  • 思想:
    • 类似于“leave-one-out”的交叉验证方法,定义第i个Jackknife样本为丢掉第i个样本后的剩余样本。

bias的Jackknife估计

  • ==这里注意证明==:为什么前面要有一个(n-1)

se的Jackknife估计

  • ==这里注意证明==:为什么前面要有一个(n-1)/n

失效情况:

  • 当估计量$\hat \theta$不够平滑(也即:当数据的微小改变就能引起估计量的巨大波动时),Jackknife方法就可能会失效。
  • 中位数就是其中一个不平滑统计量的例子。

Introduction to Markov Chain Monte Carlo Methods

  • 重点:构造合适的马氏链
  • 概括:先构造一个马氏链,在运行此链足够长时间后,该马氏链就会达到 平稳状态,再利用构造出来的分布去用MC方法求积分。
  • 不可约+非周期+正常返 ⟹ 唯一的平稳分布。

The Metropolis-Hastings Algorithm

  • 目的:对于马氏链,在给定一个$X_t$所处的状态下,产生下一步的状态$X_{t+1}$。
  • 算法步骤:
    • 构造一个合适的Proposal distribution $g(\cdot \mid X_t)$
    • 从$g(\cdot \mid X_t)$中产生Y
    • 若Y被接受,则$X_{t+1}=Y$,否则$X_{t+1}=X_t$
  • 关于Proposal distribution应该怎么选?
    • 要使得马氏链最后的平稳分布即为expected distribution f
    • 需要满足三大条件:不可约,正常返,非周期。
    • (往往是:和目标分布具有相同的支撑集的一些分布会满足这些条件。

Metropolis-Hastings Sampler

  • 算法步骤:

  • 注意:

    • 接受概率为:
      • $\alpha(X_t,Y)=min\left(1,\frac{f(Y)g(X_t\mid Y)}{f(X_t)g(Y\mid X_t)}\right)$
    • 条件:目标分布f需要正比于某个常数,即使这个常数不知道(才能上下约掉
    • 不管接受还是拒绝,马氏链在时间点上都会前进一步(如果是拒绝的情况,则是短的水平平移线
  • ==这里注意证明==:在MH算法下得到的马氏链的平稳分布为f

The Metropolis sampler

  • MH的一个特例,但与MH算法的不同在于,Metropolis算法中的Proposal distribution必须要是对称的。

Random Walk Metropolis

  • 是Metropolis算法的一个特例(所以Proposal distribution也要是对称的 $g(Y\mid X_t)=g(\mid X_t - Y\mid)$
  • 多了一个增量项Z,有$Y=X_t+Z$
  • 拒绝率最好是在$[0.15,0.5]$之间

The Independence Sampler

  • MH算法的一个特例,与MH算法不同的是,Independence算法中的Proposal distribution不依赖于链的前一步状态

  • 失效:当Proposal distribution与目标分布差别很大时表现较差。

The Gibbs Sampler

  • MH的一个特例,与MH算法不同的是,它主要用于目标分布是多元分布的场合。

  • 一般爱考的就是二元正态。

EM Algorithm

  • 目的:利用观测的数据,不断迭代估计未知参数。
  • 思想:replace a difficult maximization with a sequence of easier maximization whose limit is the answer to the original problem
  • 步骤:
    • E步:计算 $Q(\theta,\theta^*)$,对完全似然进行一个期望
    • M步:根据参数最大化$Q(\theta,\theta^{(t)})$(即:对参数求导),对于完全的数据进行一个ML估计。
  • 例子:
    • Exponential families
    • Finite normal mixtures
  • 收敛性证明:
    • EM每一次迭代都会增加对数似然

Monitoring Convergence

The Gelman-Rubin Method

  • 前提:从不同的初始值出发,链达到平稳后,应该表现的是一样的。更准确的说,Gelman 和 Rubin 指出链内的方差和链之间的方差应该是相同的。

Comments

Content