R（或任何语言）中偏态正态分布的非线性最小二乘回归
Nonlinear least squares regression of skewed normal distribution in R (or any language)
第一次海报。如有使用不当的礼仪或词汇，请提前致歉。
我有来自 USGS 河流调查的化学浓度 (y) 与时间 (x) 的时间序列数据。它表现出我想通过非线性最小二乘回归建模的偏斜正态分布。我能够将正态分布曲线拟合到数据，但似乎无法将 "skewness" 合并到模型中。
我从 Whuber 在这里给出的答案得出了我的正态分布拟合...线性回归最佳多项式（或更好的使用方法）？
我的数据和代码...
y <- c(0.532431978850729, 0.609737363640599, 0.651964078008195, 0.657368066358271,
0.741496240155044, 0.565435828629966, 0.703655525439792, 0.718855614453251,
0.838983191559565, 0.743767469276213, 0.860155614137561, 0.81923941209205,
1.07899884812998, 0.950877380129941, 1.01284743983765, 1.11717867112622,
1.08452873942528, 1.14640319037414, 1.35601176845714, 1.55587090166098,
1.81936731953165, 1.79952819117948, 2.27965075864338, 2.92158756334143,
3.28092981974249, 1.09884083379528, 4.52126319475028, 5.50589160306292,
6.48951979830975, 7.61196542128105, 9.56700470248019, 11.0814901164772,
13.3072954022821, 13.8519364143597, 11.4108376964234, 8.72143939873907,
5.12221325838613, 2.58106436004881, 1.0642701141608, 0.44945378376047,
0.474569233285229, 0.128299654944011, 0.432876244482592, 0.445456125461339,
0.435530646939433, 0.337503495863836, 0.456525976632425, 0.35851011819921,
0.525854215793115, 0.381206935673774, 0.548351975353343, 0.365384673834335,
0.418990479166088, 0.50039125911365, 0.490696977485334, 0.376809405620949,
0.484559448760701, 0.569111550743562, 0.439671715276438, 0.353621820313257,
0.444241243031233, 0.415197754444015, 0.474852839357701, 0.462144150397257,
0.535339727332139, 0.480714031175711)
#creating an arbitrary vector to represent time
x <- seq(1,length(y), by=1)
#model of normal distribution
f <- function(x, theta) {
m <- theta[1]; s <- theta[2]; a <- theta[3]; b <- theta[4];
a*exp(-0.5*((x-m)/s)^2) + b
}
# Estimate some starting values.
m.0 <- x[which.max(y)]; s.0 <- (max(x)-min(x))/4; b.0 <- min(y); a.0 <- (max(y)-min(y))
# Do the fit. (It takes no time at all.)
fit <- nls(y ~ f(x,c(m,s,a,b)), data.frame(x,y), start=list(m=m.0, s=s.0, a=a.0, b=b.0))
# Display the estimated location of the peak and its SE.
summary(fit)$parameters["m", 1:2]
par(mfrow=c(1,1))
plot(c(x,0),c(y,f(coef(fit)["m"],coef(fit))), main="Data", type="n",
xlab="Time", ylab="Concentration")
curve(f(x, coef(fit)), add=TRUE, col="Red", lwd=2)
points(x,y, pch=19)
那么，关于如何调整模型以适应偏度，有什么建议吗？
干杯，
杰米
数据是河流水样中某些化合物的浓度与时间的关系，不是吗？如果我绘制 y vs x，假设样本是定期采集的，我会看到一个浓度峰值，因此时间依赖性似乎是某种物理 and/or 化学现象，可以建模为 y = f( b, x) + e，其中 f 是 chemical/physical 现象的参数 b 的函数，x 代表时间。 e 项是随机误差，在化学中通常样品是独立测量的，因此 e ~ N(0, s^2)。然后用 nls 拟合 f(b, x)。
可以使用广义加性模型 (GAM) 吗？ GAM 功能强大且灵活，但难以解释模型系数。所以决定将取决于你的目的。如果目的是评估趋势，或者目的是预测浓度（在已知时间范围内），那么GAM可能是一个不错的选择。
library(mgcv)
library(ggplot2)
dat <- data.frame(x = 1:length(y), y = y)
fit_gam <- gam(y ~ s(x, k = 20), data = dat)
ggplot(dat, aes(x = x, y = y)) +
geom_point() +
geom_line(data = data.frame(x = x, y = fit_gam$fitted.values),
color = "red") +
ggtitle("Data") +
xlab("Cocentration") +
ylab("Time") +
theme_bw() +
theme(panel.grid = element_blank())
下面是应用stat_smooth的另一个选项来拟合同一个GAM模型。
ggplot(dat, aes(x = x, y = y)) +
geom_point() +
stat_smooth(method = "gam", formula = y ~ s(x, bs = "tp", k = 20)) +
ggtitle("Data") +
xlab("Cocentration") +
ylab("Time") +
theme_bw() +
theme(panel.grid = element_blank())
我和 python 的一个天才朋友谈过，他帮助我构建了右偏正态分布方程。我已经在下面发布了 R 脚本。
我想做的是用右偏分布模型替换正态分布模型。吸引我的不是脚本编写，而是我为右偏分布编写一般方程的能力（我的伙伴也是数学天才）。
我对 www 竖起大拇指，因为出于所有深入的目的，他们回答了我的问题。我喜欢他们也采用了使用 GAM 的不同方法，尽管我对模型产生的系数很感兴趣。
我的下一个计划是整合模型曲线下的面积，以及置信区间曲线下的面积。
初次使用 Whosebug 的体验很好。谢谢你们。
f <- function(x, theta) {
m <- theta[1]; s <- theta[2]; a <- theta[3]; b <- theta[4]; k <- theta[5]
a*exp(k*((x - m))/s - sqrt(((x - m))/s*((x - m))/s+1)) + b
}
# Estimate some starting values.
m.0 <- x[which.max(y)]; s.0 <- (max(x)-min(x))/4; b.0 <- min(y); a.0 <- (max(y)-min(y)); k.0 <- -0.5
# Do the fit. (It takes no time at all.)
fit <- nls(y ~ f(x,c(m,s,a,b, k)), data.frame(x,y), start=list(m=m.0, s=s.0, a=a.0, b=b.0, k=k.0))
# Display the estimated location of the peak and its SE.
summary(fit)$parameters["m", 1:2]
par(mfrow=c(1,1))
plot(c(x,0),c(y,f(coef(fit)["m"],coef(fit))), main="Data", type="n",
xlab="Time", ylab="Concentration")
curve(f(x, coef(fit)), add=TRUE, col="Red", lwd=2)
points(x,y, pch=19)
第一次海报。如有使用不当的礼仪或词汇，请提前致歉。
我有来自 USGS 河流调查的化学浓度 (y) 与时间 (x) 的时间序列数据。它表现出我想通过非线性最小二乘回归建模的偏斜正态分布。我能够将正态分布曲线拟合到数据，但似乎无法将 "skewness" 合并到模型中。
我从 Whuber 在这里给出的答案得出了我的正态分布拟合...线性回归最佳多项式（或更好的使用方法）？
我的数据和代码...
y <- c(0.532431978850729, 0.609737363640599, 0.651964078008195, 0.657368066358271,
0.741496240155044, 0.565435828629966, 0.703655525439792, 0.718855614453251,
0.838983191559565, 0.743767469276213, 0.860155614137561, 0.81923941209205,
1.07899884812998, 0.950877380129941, 1.01284743983765, 1.11717867112622,
1.08452873942528, 1.14640319037414, 1.35601176845714, 1.55587090166098,
1.81936731953165, 1.79952819117948, 2.27965075864338, 2.92158756334143,
3.28092981974249, 1.09884083379528, 4.52126319475028, 5.50589160306292,
6.48951979830975, 7.61196542128105, 9.56700470248019, 11.0814901164772,
13.3072954022821, 13.8519364143597, 11.4108376964234, 8.72143939873907,
5.12221325838613, 2.58106436004881, 1.0642701141608, 0.44945378376047,
0.474569233285229, 0.128299654944011, 0.432876244482592, 0.445456125461339,
0.435530646939433, 0.337503495863836, 0.456525976632425, 0.35851011819921,
0.525854215793115, 0.381206935673774, 0.548351975353343, 0.365384673834335,
0.418990479166088, 0.50039125911365, 0.490696977485334, 0.376809405620949,
0.484559448760701, 0.569111550743562, 0.439671715276438, 0.353621820313257,
0.444241243031233, 0.415197754444015, 0.474852839357701, 0.462144150397257,
0.535339727332139, 0.480714031175711)
#creating an arbitrary vector to represent time
x <- seq(1,length(y), by=1)
#model of normal distribution
f <- function(x, theta) {
m <- theta[1]; s <- theta[2]; a <- theta[3]; b <- theta[4];
a*exp(-0.5*((x-m)/s)^2) + b
}
# Estimate some starting values.
m.0 <- x[which.max(y)]; s.0 <- (max(x)-min(x))/4; b.0 <- min(y); a.0 <- (max(y)-min(y))
# Do the fit. (It takes no time at all.)
fit <- nls(y ~ f(x,c(m,s,a,b)), data.frame(x,y), start=list(m=m.0, s=s.0, a=a.0, b=b.0))
# Display the estimated location of the peak and its SE.
summary(fit)$parameters["m", 1:2]
par(mfrow=c(1,1))
plot(c(x,0),c(y,f(coef(fit)["m"],coef(fit))), main="Data", type="n",
xlab="Time", ylab="Concentration")
curve(f(x, coef(fit)), add=TRUE, col="Red", lwd=2)
points(x,y, pch=19)
那么，关于如何调整模型以适应偏度，有什么建议吗？
干杯， 杰米
数据是河流水样中某些化合物的浓度与时间的关系，不是吗？如果我绘制 y vs x，假设样本是定期采集的，我会看到一个浓度峰值，因此时间依赖性似乎是某种物理 and/or 化学现象，可以建模为 y = f( b, x) + e，其中 f 是 chemical/physical 现象的参数 b 的函数，x 代表时间。 e 项是随机误差，在化学中通常样品是独立测量的，因此 e ~ N(0, s^2)。然后用 nls 拟合 f(b, x)。
可以使用广义加性模型 (GAM) 吗？ GAM 功能强大且灵活，但难以解释模型系数。所以决定将取决于你的目的。如果目的是评估趋势，或者目的是预测浓度（在已知时间范围内），那么GAM可能是一个不错的选择。
library(mgcv)
library(ggplot2)
dat <- data.frame(x = 1:length(y), y = y)
fit_gam <- gam(y ~ s(x, k = 20), data = dat)
ggplot(dat, aes(x = x, y = y)) +
geom_point() +
geom_line(data = data.frame(x = x, y = fit_gam$fitted.values),
color = "red") +
ggtitle("Data") +
xlab("Cocentration") +
ylab("Time") +
theme_bw() +
theme(panel.grid = element_blank())
下面是应用stat_smooth的另一个选项来拟合同一个GAM模型。
ggplot(dat, aes(x = x, y = y)) +
geom_point() +
stat_smooth(method = "gam", formula = y ~ s(x, bs = "tp", k = 20)) +
ggtitle("Data") +
xlab("Cocentration") +
ylab("Time") +
theme_bw() +
theme(panel.grid = element_blank())
我和 python 的一个天才朋友谈过，他帮助我构建了右偏正态分布方程。我已经在下面发布了 R 脚本。
我想做的是用右偏分布模型替换正态分布模型。吸引我的不是脚本编写，而是我为右偏分布编写一般方程的能力（我的伙伴也是数学天才）。
我对 www 竖起大拇指，因为出于所有深入的目的，他们回答了我的问题。我喜欢他们也采用了使用 GAM 的不同方法，尽管我对模型产生的系数很感兴趣。
我的下一个计划是整合模型曲线下的面积，以及置信区间曲线下的面积。
初次使用 Whosebug 的体验很好。谢谢你们。
f <- function(x, theta) {
m <- theta[1]; s <- theta[2]; a <- theta[3]; b <- theta[4]; k <- theta[5]
a*exp(k*((x - m))/s - sqrt(((x - m))/s*((x - m))/s+1)) + b
}
# Estimate some starting values.
m.0 <- x[which.max(y)]; s.0 <- (max(x)-min(x))/4; b.0 <- min(y); a.0 <- (max(y)-min(y)); k.0 <- -0.5
# Do the fit. (It takes no time at all.)
fit <- nls(y ~ f(x,c(m,s,a,b, k)), data.frame(x,y), start=list(m=m.0, s=s.0, a=a.0, b=b.0, k=k.0))
# Display the estimated location of the peak and its SE.
summary(fit)$parameters["m", 1:2]
par(mfrow=c(1,1))
plot(c(x,0),c(y,f(coef(fit)["m"],coef(fit))), main="Data", type="n",
xlab="Time", ylab="Concentration")
curve(f(x, coef(fit)), add=TRUE, col="Red", lwd=2)
points(x,y, pch=19)