如何在R中为偏移对数正态分布数据拟合曲线并获取方程

数据准备 #

你提供的数据集可通过以下R代码重建：

df = data.frame(
  x = c(400.0000, 410.3448, 420.6897, 431.0345, 441.3793, 451.7241, 462.0690, 472.4138, 482.7586, 493.1034, 503.4483, 513.7931, 524.1379, 534.4828, 544.8276, 555.1724, 565.5172, 575.8621, 586.2069, 596.5517, 606.8966, 617.2414, 627.5862, 637.9310, 648.2759, 658.6207, 668.9655, 679.3103, 689.6552, 700.0000),
  y = c(0.0158, 0.0373, 0.0379, 0.0303, 0.0250, 0.0212, 0.0234, 0.0295, 0.0329, 0.0364, 0.0469, 0.0678, 0.0851, 0.0693, 0.0508, 0.0482, 0.0617, 0.2510, 0.6570, 0.7360, 0.6690, 0.5810, 0.5060, 0.4390, 0.3640, 0.2980, 0.2390, 0.1480, 0.0496, 0.0122)
)

拟合偏移对数正态分布 #

偏移对数正态分布的核心拟合公式为（适配你的数据量级，加入缩放系数）：
$$y = \frac{A}{(x - \xi)\sigma\sqrt{2\pi}} e^{-\frac{(\ln(x - \xi) - \mu)^2}{2\sigma2}}$$
参数说明：

• $A$：缩放系数，匹配你的y值范围
• $\xi$：偏移参数，控制分布的位置偏移
• $\mu$：对数均值，对应分布峰值的对数位置
• $\sigma$：对数标准差，控制分布的宽窄

步骤1：安装并加载工具包 #

使用nls2自动搜索初始参数（避免手动猜测的麻烦）：

# 首次运行需安装包
install.packages("nls2")
library(nls2)

步骤2：定义模型函数 #

shifted_lognormal <- function(x, A, xi, mu, sigma) {
  term1 <- A / ((x - xi) * sigma * sqrt(2 * pi))
  term2 <- exp(-(log(x - xi) - mu)^2 / (2 * sigma^2))
  term1 * term2
}

步骤3：拟合模型 #

先通过网格搜索获取初始参数，再用非线性最小二乘法优化：

# 设置参数搜索范围
param_grid <- expand.grid(
  A = seq(10, 100, by = 10),
  xi = seq(350, 400, by = 10),
  mu = seq(5, 7, by = 0.5),
  sigma = seq(0.1, 1, by = 0.1)
)

# 网格搜索找初始值
nls_start <- nls2(y ~ shifted_lognormal(x, A, xi, mu, sigma),
                  data = df,
                  start = param_grid,
                  algorithm = "brute-force")

# 拟合最终模型
final_model <- nls(y ~ shifted_lognormal(x, A, xi, mu, sigma),
                   data = df,
                   start = as.list(coef(nls_start)))

# 查看拟合参数与统计结果
summary(final_model)

步骤4：验证拟合效果 #

将拟合曲线与原始数据可视化，检查匹配度：

plot(df$x, df$y, pch = 16, col = "blue", main = "偏移对数正态分布拟合", xlab = "x", ylab = "y")
x_seq <- seq(400, 700, by = 1)
y_pred <- predict(final_model, newdata = data.frame(x = x_seq))
lines(x_seq, y_pred, col = "red", lwd = 2)
legend("topright", legend = c("原始数据", "拟合曲线"), col = c("blue", "red"), pch = c(16, NA), lty = c(NA, 1))

步骤5：估算任意x值的y值 #

直接用拟合模型预测400-700区间内的任意x值，例如估算x=550的y值：

predict(final_model, newdata = data.frame(x = 550))

替代方案：无nls2时的手动初始参数 #

如果nls2安装有问题，可根据数据特征手动设置初始参数，再拟合：

final_model <- nls(y ~ shifted_lognormal(x, A, xi, mu, sigma),
                   data = df,
                   start = list(A=60, xi=390, mu=5.3, sigma=0.2))

后续步骤与上述一致。

内容的提问来源于stack exchange，提问作者dark-walrus