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

Question

Первый раз постер. Заранее извиняюсь, если я использую неправильный этикет или словарь.

У меня есть данные временного ряда химической концентрации (y) и времени (x) из речного обследования USGS. Это демонстрирует косое нормальное распределение, которое я хотел бы смоделировать с помощью нелинейной регрессии наименьших квадратов. Я могу приспособить нормальную кривую распределения к данным, но, похоже, не могу включить "асимметрию" в модель.

Я пришел к своему нормальному распределению соответствия из ответа, данного здесь Уубером. Линейная регрессия - лучший полином (или лучший подход к использованию)?

мои данные и код ...

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)

Итак, есть ли какие-либо предложения о том, как настроить модель с учетом асимметрии?

Ура, джем ie

Answer 1

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

Answer 2

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

Answer 3

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

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

Я дал большие пальцы до www because для всех интенсивных целей, они ответили на мой вопрос. Мне нравится, что они также использовали другой подход, используя 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)