如何用逆变换法模拟伯努利分布最大值函数及在R中生成M变量？

好的，我来拆解这两个问题，一步步给你讲清楚怎么用逆变换法实现，包括R语言的代码。

首先明确：这里的“最大值函数”应该是指n个独立伯努利变量的最大值——给定$B_1 \sim \text{Bernoulli}(p_1), B_2 \sim \text{Bernoulli}(p_2), ..., B_n \sim \text{Bernoulli}(p_n)$且相互独立，求$M_n = \max{B_1,B_2,...,B_n}$的模拟方法。

伯努利变量的取值只有0和1，所以最大值$M_n$的取值也只有0或1：

• $M_n=0$当且仅当所有$B_i=0$，概率为$P(M_n=0)=\prod_{i=1}^n (1-p_i)$
• $M_n=1$当且仅当至少有一个$B_i=1$，概率为$P(M_n=1)=1 - \prod_{i=1}^n (1-p_i)$

逆变换法的核心是利用均匀分布的逆CDF生成目标变量。对于$M_n$，它的CDF是：

• $F(x)=0$，当$x<0$
• $F(x)=\prod_{i=1}^n (1-p_i)$，当$0 \leq x <1$
• $F(x)=1$，当$x \geq1$

具体步骤：

1. 生成一个服从$\text{Uniform}(0,1)$的随机数$U$
2. 如果$U \leq \prod_{i=1}^n (1-p_i)$，则$M_n=0$；否则$M_n=1$

对应的R代码实现：

# 输入：p_vec是包含每个伯努利变量参数p的向量
generate_max_bernoulli <- function(p_vec) {
  # 计算所有变量取0的概率
  prob_all_zero <- prod(1 - p_vec)
  # 生成均匀分布随机数
  U <- runif(1)
  # 根据逆变换规则返回结果
  if (U <= prob_all_zero) {
    return(0)
  } else {
    return(1)
  }
}

# 示例：3个伯努利变量，参数分别为0.3,0.5,0.2
set.seed(456)
generate_max_bernoulli(c(0.3, 0.5, 0.2))

这个变量$M$的含义是：在无限序列$B_0,B_1,B_2,...$中，所有取值为1的$B_k$对应的$k$的最大值；如果所有$B_k$都为0，那么$M=-1$（这个定义避免歧义）。每个$B_k \sim \text{Bernoulli}(p_k)$，且序列相互独立。

第一步：推导M的CDF和逆变换规则 #

首先明确$M$的概率分布：

• $P(M=-1) = $所有$B_k$都为0的概率$ = \prod_{k=0}^\infty (1-p_k)$（记为$Q_0$，要求这个无穷乘积收敛）
• 对于非负整数$m$，$P(M=m) = P(B_m=1 \text{ 且 所有}k>m\text{的}B_k=0) = p_m \times \prod_{k=m+1}^\infty (1-p_k)$（记$\prod_{k=m+1}^\infty (1-p_k)$为$Q_{m+1}$）

$M$的CDF $F(t)$为：

• $F(-1) = Q_0$
• 对于$m\geq0$，$F(m) = P(M\leq m) = $所有$k>m$的$B_k=0$的概率$ = Q_{m+1}$

注意$Q_m$是递增序列（因为$Q_{m+1}=Q_m/(1-p_m)$，而$1-p_m \leq1$，所以$Q_{m+1}\geq Q_m$），且$Q_m$最终趋近于1（当$p_k$足够快地趋近于0时）。

逆变换法步骤：

1. 生成$U \sim \text{Uniform}(0,1)$
2. 如果$U \leq Q_0$，说明所有$B_k$都是0，返回$M=-1$
3. 否则，从$m=0$开始递推计算$Q_{m+1}$，找到最小的$m$使得$F(m)=Q_{m+1}\geq U$，此时$M=m$（因为$F(m-1)=Q_m < U$，说明$M>m-1$）

第二步：R语言实现 #

因为涉及无穷乘积，我们需要用递推的方式计算，直到乘积收敛（变化足够小）。下面是完整的实现代码：

# 辅助函数：计算Q₀=∏_{k=0}^∞(1-p(k))，直到乘积收敛
compute_Q0 <- function(p, tol = 1e-15) {
  Q <- 1
  k <- 0
  while(TRUE) {
    term <- 1 - p(k)
    # 如果当前项几乎是1，再乘下去不会改变Q，停止
    if (term >= 1 - tol) break
    new_Q <- Q * term
    # 如果Q的变化小于阈值，停止
    if (abs(new_Q - Q) < tol) break
    Q <- new_Q
    k <- k + 1
  }
  return(Q)
}

# 生成M的主函数
generate_M <- function(p, tol = 1e-15) {
  # 生成均匀分布随机数
  U <- runif(1)
  # 计算所有B_k为0的概率Q₀
  Q0 <- compute_Q0(p, tol)
  
  # 情况1：所有B_k都是0
  if (U <= Q0) {
    return(-1)
  }
  
  # 情况2：存在至少一个B_k=1，递推找M
  Q_m <- Q0  # 初始Q_m=Q₀=∏_{k=0}^∞(1-p(k))
  m <- 0
  cum_prob <- Q0  # 累积概率初始为P(M=-1)
  
  while(TRUE) {
    current_p <- p(m)
    # 处理p(m)接近1的情况（避免除以0）
    if (1 - current_p < tol) {
      return(m)
    }
    
    # 计算Q_{m+1}=∏_{k=m+1}^∞(1-p(k))
    Q_m_plus_1 <- Q_m / (1 - current_p)
    # 计算P(M=m)
    prob_m <- current_p * Q_m_plus_1
    # 更新累积概率
    cum_prob <- cum_prob + prob_m
    
    # 如果U落在当前累积概率范围内，返回m
    if (U <= cum_prob) {
      return(m)
    }
    
    # 如果累积概率已经接近1，剩下的概率可以忽略，返回当前m
    if (cum_prob >= 1 - tol) {
      return(m)
    }
    
    # 更新Q_m为Q_{m+1}，m加1继续循环
    Q_m <- Q_m_plus_1
    m <- m + 1
  }
}

# 示例：定义p(k)=0.5^(k+1)，即p0=0.5,p1=0.25,p2=0.125,...
p_example <- function(k) 0.5^(k+1)
set.seed(789)
# 生成10个M的样本
replicate(10, generate_M(p_example))

代码说明 #

• compute_Q0函数用来计算无穷乘积$Q_0$，通过循环直到乘积的变化小于阈值，避免无限循环
• generate_M函数先判断是否所有$B_k$都是0，否则递推计算每个$m$对应的概率，直到找到$U$所在的区间
• 加入了对$p(m)$接近1的特殊处理，避免除以0的错误
• 当累积概率接近1时，停止循环，因为剩下的概率可以忽略不计

内容的提问来源于stack exchange，提问作者Narsis