使用 R 拟合“多峰”对数正态分布
Fit “multimodal” lognormal distributions using R
我的问题与 类似,但我想在 R 中完成。数据框是
x<-c(0.35,0.46,0.60,0.78,1.02,1.34,1.76,2.35,3.17,4.28,5.77,7.79,10.50,14.20,19.10,25.80)
y<-c(32.40,43.00,37.20,26.10,17.40,14.00,19.90,36.90,48.60,55.30,64.60,70.20,63.90,47.60,22.70,10.30)
df<-data.frame(x,y)
plot(df,log='xy')
这里绘制的是数据的样子。在 x 尺度的单位中,有一种模式约为 0.5,另一种模式约为 8。
如何将“多峰”对数正态分布拟合到此类数据(在本例中有 2 条曲线)?这是我尝试过的。非常感谢任何解决问题的帮助或指导。
ggplot(data=df, aes(x=x, y=y)) +
geom_point() +
stat_smooth(method="nls",
formula=y ~ a*dlnorm(x, meanlog=8, sdlog=2.7),
method.args = list(start=c(a=2e6)),
se=FALSE,color = "red", linetype = 2)+
scale_x_log10()+
scale_y_log10()
我假设你想要 nls
。您可以通过在方程中定义两个参数来考虑两种模式,例如 a
和 b
。定义两个 start=
ing 值。 (请注意,我此时只是猜测了所有值。)
fit <- nls(y ~ a*dlnorm(x, meanlog=.5, sdlog=.5) + b*dlnorm(x, meanlog=8, sdlog=2.7),
data=df1, start=list(a=1, b=1))
summary(fit)
# Formula: y ~ a * dlnorm(x, meanlog = 0.5, sdlog = 0.5) + b * dlnorm(x,
# meanlog = 8, sdlog = 2.7)
#
# Parameters:
# Estimate Std. Error t value Pr(>|t|)
# a -81.97 16.61 -4.934 0.00022 ***
# b 30695.42 2417.90 12.695 4.53e-09 ***
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Residual standard error: 11.92 on 14 degrees of freedom
#
# Number of iterations to convergence: 1
# Achieved convergence tolerance: 4.507e-07
fitted()
已经根据数据框的 x
值给出了 y
的拟合值。
fitted(fit)
# [1] 45.56775 44.59130 38.46212 27.34071 15.94205 12.76579 21.31640
# [8] 36.51385 48.68786 53.60069 53.56958 51.40254 48.41267 44.95541
# [15] 41.29045 37.41424
# attr(,"label")
# [1] "Fitted values"
您也可以为此使用 predict()
。
stopifnot(all.equal(predict(fit), as.numeric(fitted(fit))))
但是,为了获得更平滑的线条,您需要 predict
ion(即 y
值)沿 x
轴的一组更精细的 x
值。
plot(df1, log='xy')
x.seq <- seq(0, max(df$x), .1)
lines(x=x.seq, y=predict(fit, newdata=data.frame(x=x.seq)), col=2)
旁注: 即使这很常见,通过命名你的数据框 df
你使用的是相同的名称对于 F 分布的密度函数 df()
,这可能会导致混淆!为此,我使用了 df1
。
数据:
df1 <- structure(list(x = c(0.35, 0.46, 0.6, 0.78, 1.02, 1.34, 1.76,
2.35, 3.17, 4.28, 5.77, 7.79, 10.5, 14.2, 19.1, 25.8), y = c(32.4,
43, 37.2, 26.1, 17.4, 14, 19.9, 36.9, 48.6, 55.3, 64.6, 70.2,
63.9, 47.6, 22.7, 10.3)), class = "data.frame", row.names = c(NA,
-16L))
我的问题与
x<-c(0.35,0.46,0.60,0.78,1.02,1.34,1.76,2.35,3.17,4.28,5.77,7.79,10.50,14.20,19.10,25.80)
y<-c(32.40,43.00,37.20,26.10,17.40,14.00,19.90,36.90,48.60,55.30,64.60,70.20,63.90,47.60,22.70,10.30)
df<-data.frame(x,y)
plot(df,log='xy')
这里绘制的是数据的样子。在 x 尺度的单位中,有一种模式约为 0.5,另一种模式约为 8。
如何将“多峰”对数正态分布拟合到此类数据(在本例中有 2 条曲线)?这是我尝试过的。非常感谢任何解决问题的帮助或指导。
ggplot(data=df, aes(x=x, y=y)) +
geom_point() +
stat_smooth(method="nls",
formula=y ~ a*dlnorm(x, meanlog=8, sdlog=2.7),
method.args = list(start=c(a=2e6)),
se=FALSE,color = "red", linetype = 2)+
scale_x_log10()+
scale_y_log10()
我假设你想要 nls
。您可以通过在方程中定义两个参数来考虑两种模式,例如 a
和 b
。定义两个 start=
ing 值。 (请注意,我此时只是猜测了所有值。)
fit <- nls(y ~ a*dlnorm(x, meanlog=.5, sdlog=.5) + b*dlnorm(x, meanlog=8, sdlog=2.7),
data=df1, start=list(a=1, b=1))
summary(fit)
# Formula: y ~ a * dlnorm(x, meanlog = 0.5, sdlog = 0.5) + b * dlnorm(x,
# meanlog = 8, sdlog = 2.7)
#
# Parameters:
# Estimate Std. Error t value Pr(>|t|)
# a -81.97 16.61 -4.934 0.00022 ***
# b 30695.42 2417.90 12.695 4.53e-09 ***
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Residual standard error: 11.92 on 14 degrees of freedom
#
# Number of iterations to convergence: 1
# Achieved convergence tolerance: 4.507e-07
fitted()
已经根据数据框的 x
值给出了 y
的拟合值。
fitted(fit)
# [1] 45.56775 44.59130 38.46212 27.34071 15.94205 12.76579 21.31640
# [8] 36.51385 48.68786 53.60069 53.56958 51.40254 48.41267 44.95541
# [15] 41.29045 37.41424
# attr(,"label")
# [1] "Fitted values"
您也可以为此使用 predict()
。
stopifnot(all.equal(predict(fit), as.numeric(fitted(fit))))
但是,为了获得更平滑的线条,您需要 predict
ion(即 y
值)沿 x
轴的一组更精细的 x
值。
plot(df1, log='xy')
x.seq <- seq(0, max(df$x), .1)
lines(x=x.seq, y=predict(fit, newdata=data.frame(x=x.seq)), col=2)
旁注: 即使这很常见,通过命名你的数据框 df
你使用的是相同的名称对于 F 分布的密度函数 df()
,这可能会导致混淆!为此,我使用了 df1
。
数据:
df1 <- structure(list(x = c(0.35, 0.46, 0.6, 0.78, 1.02, 1.34, 1.76,
2.35, 3.17, 4.28, 5.77, 7.79, 10.5, 14.2, 19.1, 25.8), y = c(32.4,
43, 37.2, 26.1, 17.4, 14, 19.9, 36.9, 48.6, 55.3, 64.6, 70.2,
63.9, 47.6, 22.7, 10.3)), class = "data.frame", row.names = c(NA,
-16L))