生成指定均值、标准差且限区间的对数正态随机数问题求解

你遇到的核心问题是：直接使用原始对数正态分布的meanlog和sdlog参数生成截断随机数时，截断操作会改变分布的矩（均值、标准差），导致最终结果不符合预期。要解决这个问题，需要先根据目标均值、标准差和截断范围，反推截断后对数正态分布对应的meanlog和sdlog参数，再用这些参数生成随机数。

具体解决方案 #

方法一：矩匹配优化法（精准匹配目标矩） #

通过优化算法求解满足目标均值、标准差的截断对数正态分布参数，步骤如下：

1. 编写函数计算截断对数正态分布的均值和标准差

trunc_lnorm_moments <- function(params, L, U) {
  mu <- params[1]
  sigma <- params[2]
  
  # 转换为标准正态分布的分位点
  z_low <- (log(L) - mu) / sigma
  z_high <- (log(U) - mu) / sigma
  
  # 计算标准正态的密度和累积分布值
  phi_low <- dnorm(z_low)
  phi_high <- dnorm(z_high)
  Phi_low <- pnorm(z_low)
  Phi_high <- pnorm(z_high)
  
  # 截断后的均值
  mean_trunc <- exp(mu + sigma^2/2) * (pnorm(z_high - sigma) - pnorm(z_low - sigma)) / (Phi_high - Phi_low)
  # 截断后的标准差
  var_trunc <- exp(2*mu + sigma^2) * (exp(sigma^2)*(pnorm(z_high - 2*sigma) - pnorm(z_low - 2*sigma)) - (pnorm(z_high - sigma) - pnorm(z_low - sigma))^2) / (Phi_high - Phi_low)^2
  sd_trunc <- sqrt(var_trunc)
  
  return(c(mean_trunc, sd_trunc))
}

2. 定义优化目标函数（最小化计算矩与目标矩的误差）

objective <- function(params, target_mean, target_sd, L, U) {
  calculated <- trunc_lnorm_moments(params, L, U)
  (calculated[1] - target_mean)^2 + (calculated[2] - target_sd)^2
}

3. 求解最优参数并生成随机数

library(EnvStats)
set.seed(123)

# 替换为你的目标均值、标准差，以及截断范围
target_mean <- 150000  # 示例目标均值
target_sd <- 80000     # 示例目标标准差
L <- exp(10)
U <- exp(13)

# 生成初始参数猜测（基于原始对数正态矩公式）
initial_mu <- log(target_mean) - 0.5 * log(1 + (target_sd/target_mean)^2)
initial_sigma <- sqrt(log(1 + (target_sd/target_mean)^2))

# 优化求解最优参数
opt_result <- optim(c(initial_mu, initial_sigma), objective, 
                    target_mean=target_mean, target_sd=target_sd, L=L, U=U)
optimal_mu <- opt_result$par[1]
optimal_sigma <- opt_result$par[2]

# 生成符合要求的截断随机数
random_nums <- rlnormTrunc(100000, meanlog=optimal_mu, sdlog=optimal_sigma, min=L, max=U)

# 验证结果
mean(random_nums)
sd(random_nums)

方法二：接受-拒绝抽样法（直观但效率较低） #

如果截断范围较宽，可先生成原始对数正态随机数，再过滤掉超出范围的数值。注意：这种方法得到的截断分布矩会偏离原始分布的目标矩，仅适用于对矩要求不严格的场景。

set.seed(123)
n <- 100000
L <- exp(10)
U <- exp(13)

# 基于目标均值、标准差计算原始对数正态参数
target_mean <- 150000
target_sd <- 80000
mu <- log(target_mean) - 0.5 * log(1 + (target_sd/target_mean)^2)
sigma <- sqrt(log(1 + (target_sd/target_mean)^2))

# 抽样并过滤
random_nums <- numeric(0)
while(length(random_nums) < n) {
  candidates <- rlnorm(n * 1.2, meanlog=mu, sdlog=sigma)  # 多生成20%提高效率
  valid <- candidates >= L & candidates <= U
  random_nums <- c(random_nums, candidates[valid])
}
random_nums <- random_nums[1:n]

# 验证结果
mean(random_nums)
sd(random_nums)

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