Нелинейная регрессия по методу наименьших квадратов асимметричного нормального распределения в R (или любом другом языке)

Можете ли вы использовать обобщенную аддитивную модель (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())

Данные представляют собой концентрацию определенного соединения в зависимости от времени в пробах воды из реки, не так ли? Если я построю зависимость y от x, предполагая, что пробы отбирались через равные промежутки времени, я увижу пик концентрации, поэтому зависимость от времени представляется каким-то физическим и/или химическим явлением, которое можно смоделировать как y = f(b, x) + e, где f — функция параметров b химических/физических явлений, а x — время. Член e - это случайная ошибка, в химии обычно образцы измеряются независимо, таким образом, e ~ N (0, s ^ ​​2). Затем вы подгоняете f(b, x) к nls.

Я поговорил с приятелем, который хорошо разбирается в питоне, и он помог мне составить правильное уравнение нормального распределения с перекосом. Я разместил сценарий R ниже.

Что я хотел сделать, так это заменить модель нормального распределения моделью распределения с правым перекосом. Меня увлекло не столько написание сценария, сколько моя способность написать общее уравнение для распределения с правильным перекосом (мой приятель тоже гений математики).

Я дал большой палец вверх www, потому что для всех интенсивных целей они ответили на мой вопрос. Мне нравится, что они также использовали другой подход, используя GAM, хотя меня интересуют коэффициенты, создаваемые моделью.

Мой следующий план состоит в том, чтобы интегрировать площадь под модельной кривой, а также площадь под кривыми доверительного интервала.

Первый опыт работы со stackoverflow был удачным. Спасибо вам всем.

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)