Logit Beta分布及其R语言随机模拟算法

Logit Beta分布是一个在广义线性模型中时常遇到的分布，通常是作为模型算法的一个中间部分。因其自身特点，在采样的时候，很容易得到取值非常大的样本点，从而导致算法发散。所以本篇简单介绍一下这个分布，并给出对其进行采样的策略。

---

Logit Beta分布

Diaconis-Ylvisaker分布是Natural Exponential Family的共轭先验分布，假设$X\sim DY(\alpha ,\kappa ;\psi )$，则
$$f(x)\propto \exp (\alpha x-\kappa \psi (x))$$

在Data Model为二项分布的广义线性模型（主要是Logit模型及其变型）中，
$$\psi (x)=\log (1+e^x)$$

从下面的定理可以看出这个分布与beta分布的关系，而这种关系经常出现在分层logit模型中，所以这个分布也被称为logit beta分布。

定理 假设$w\sim \text{Beta}(\alpha ,\kappa -\alpha )$，则
$$X=\log \frac{w}{1-w}\sim DY(\alpha ,\kappa ;\psi )$$

证明 $X$是随机变量$w$的连续可微且单调递增的变换，其逆变换为$$w=\frac{e^X}{1+e^X}$$

所以$X$的概率密度为
f ( x ) = f w ( w ( x ) ) w ′ ( x ) = Γ ( κ ) Γ ( α ) Γ ( κ − α ) ( e x 1 + e x ) α − 1 ( 1 1 + e x ) κ − α − 1 e x ( 1 + e x ) 2 = Γ ( κ ) Γ ( α ) Γ ( κ − α ) e α x ( 1 + e x ) κ = Γ ( κ ) Γ ( α ) Γ ( κ − α ) exp ⁡ ( α x − κ log ⁡ ( 1 + e x ) )

f(x)​=fw​(w(x))w′(x)=Γ(α)Γ(κ−α)Γ(κ)​(1+exex​)α−1(1+ex1​)κ−α−1(1+ex)2ex​=Γ(α)Γ(κ−α)Γ(κ)​(1+ex)κeαx​=Γ(α)Γ(κ−α)Γ(κ)​exp(αx−κlog(1+ex))​

Logit Beta分布的采样算法

Logit Beta分布的直接采样法可能遇到的问题
根据上文的定理，要得到Logit Beta分布的随机样本，只需要生成一个beta分布的样本$w\sim \text{beta}(\alpha ,\kappa -\alpha )$，然后根据变换$x=\log \frac{w}{1-w}$就可以得到logit beta分布的样本。但beta分布与这个变换的特征综合起来就会带来一些问题。

考虑$\text{beta}(\alpha ,\kappa -\alpha )$，它的核为
$$g(x)=x^{\alpha -1}(1-x{)}^{\kappa -\alpha -1},x\in (0,1)$$

如果$\alpha <1$，则当$x\to 0$时，$g(x)\to +\infty$；如果$\kappa -\alpha <1$，则当$x\to 1$时，$g(x)\to +\infty$；也就是说如果$\alpha <1$或者$\kappa -\alpha <1$，beta分布的密度都会集中到$0$或者$1$附近。

而关于变换$$h(x)=\log \frac{x}{1-x},x\in (0,1)$$

它是一个单调递增的函数，当$x\to 0$时，$h(x)\to -\infty$；当$x\to 1$时，$h(x)\to +\infty$；综合变换与beta分布的特征，如果$\alpha ,\kappa -\alpha$均小于1，直接采样法就会得到很多绝对值非常大的logit beta样本，所以用直接采样法作为算法的中间步骤极大概率导致算法无法收敛。

对直接采样法的改进
根据上文的分析，我们不希望beta分布的样本太靠近0和1，所以引入一个threshold $\delta$，记生成的beta随机样本为$w$，按下面的思路处理：

1. 如果$w\in (\delta ,1-\delta )$，就用$x=\log \frac{w}{1-w}$作为logit beta样本；
2. 如果$w\in (0,\delta ]$，用替代方案$x=-\frac{1}{\alpha }A-\log (B+e^{-\frac{A}{\alpha }})$，其中$A\sim \text{Gamma}(1,1),B\sim \text{Gamma}(\kappa -\alpha ,1)$
3. 如果$w\in [1-\delta ,1)$，用替代方案$x=\frac{1}{\kappa -\alpha }A+\log (B+e^{-\frac{A}{\kappa -\alpha }})$，其中$A\sim \text{Gamma}(1,1),B\sim \text{Gamma}(\alpha ,1)$

根据这个替代方法可以写出logit beta的采样代码（Jonathan Bradley 2018）

#' Simulate logitBeta
#'
#' This code simulates logit beta random variables
#' @param alpha A n-dimensional vector of shape parameters
#' @param kappa A second n-dimensional vector of shape parameters. kappa must be greater than alpha.
#' @return W An n-dimensional vector of simulated values
#' @export

logitbetasim<-function(alpha,kappa){

  temp1=rbeta(length(alpha),alpha,kappa-alpha)

  indics = temp1<=1e-12
  indics2 = (1-temp1)<=1e-12

  W = log(temp1/(1-temp1))

  bet2 = kappa-alpha

  #small shape and scale stuff
  if (sum(indics)>0){
    X1 = -(1/alpha[indics==1])*rgamma(sum(indics),matrix(1,sum(indics),1),matrix(1,sum(indics),1))
    Y1 = rgamma(sum(indics),bet2[indics==1],matrix(1,sum(indics),1))
    W[indics==1] = X1 - log(exp(X1)+Y1);
  }

  if (sum(indics2)>0){
    X1 = -(1/bet2[indics2==1])*rgamma(sum(indics2),matrix(1,sum(indics2),1),matrix(1,sum(indics2),1))
    Y1 = rgamma(sum(indics2),alpha[indics2==1],matrix(1,sum(indics2),1))
    W[indics2==1] = -X1 + log(exp(X1)+Y1);
  }

  return(W)
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34