在R中使用LogNormal分布进行插值的实现问题

1. 先修正代码中的基础问题 #

你的代码存在两处语法错误：meanlog <- mean(log(y2))和SDlog <- sd(log(y2))都缺少闭合括号。更关键的是，这两行逻辑错误——对数正态分布的参数应该基于**土地规模（x变量）**的对数估计，而非持有面积（y2）。

2. 核心逻辑梳理 #

我们的目标是：基于原有分组数据拟合对数正态分布，利用该分布的累积概率计算目标区间的持有面积占比，再乘以总持有量得到最终结果。

3. 完整实现步骤 #

步骤1：整理原始数据并拟合对数正态分布

首先将原始分类转换为区间边界，计算区间中点（用持有面积加权）来拟合分布参数：

library(assertthat)

# 整理原始分组数据
original_data <- data.frame(
  lower = c(0, 0.8, 1.6, 2.4, 3.2, 4.0, 6.0),
  upper = c(0.8, 1.6, 2.4, 3.2, 4.0, 6.0, Inf),
  area = c(18012, 66155, 80224, 61555, 47754, 56234, 38257)
)

# 计算区间中点，处理无限区间（将6+的中点设为8，可按需调整）
original_data$midpoint <- (original_data$lower + original_data$upper) / 2
original_data$midpoint[original_data$upper == Inf] <- 8

# 加权拟合对数正态分布的参数
weighted_meanlog <- weighted.mean(log(original_data$midpoint), w = original_data$area)
weighted_sdlog <- sqrt(weighted.mean((log(original_data$midpoint) - weighted_meanlog)^2, w = original_data$area))

步骤2：编写对数正态插值函数

定义函数计算目标区间对应的持有面积：

lnorm_interpolate_area <- function(target_lower, target_upper, total_area, meanlog, sdlog) {
  # 参数合法性校验
  assert_that(is.numeric(target_lower), is.numeric(target_upper),
              is.scalar(total_area), is.scalar(meanlog), is.scalar(sdlog),
              target_lower < target_upper)
  
  # 计算目标区间的累积概率占比
  p_lower <- ifelse(target_lower == 0, 0, plnorm(target_lower, meanlog = meanlog, sdlog = sdlog))
  p_upper <- ifelse(is.infinite(target_upper), 1, plnorm(target_upper, meanlog = meanlog, sdlog = sdlog))
  area_proportion <- p_upper - p_lower
  
  # 返回估算的持有面积
  return(total_area * area_proportion)
}

步骤3：应用函数到目标分类

# 定义目标分类区间
target_classes <- data.frame(
  lower = c(0, 1, 2, 3, 4, 5),
  upper = c(1, 2, 3, 4, 5, 10)
)

# 计算每个目标区间的持有面积
target_classes$estimated_area <- mapply(lnorm_interpolate_area,
                                        target_lower = target_classes$lower,
                                        target_upper = target_classes$upper,
                                        MoreArgs = list(total_area = 368191,
                                                        meanlog = weighted_meanlog,
                                                        sdlog = weighted_sdlog))

# 查看结果
print(target_classes)

4. 注意事项 #

• 对于6.0及以上这类无限区间，中点值的选择会影响拟合结果，可根据数据的实际分布调整（比如用行业默认的尾部中点值）。
• 如果原始数据的分组更精细，拟合的分布会更准确；当前用中点近似是分组数据下的常用方法。

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