Jak zminimalizować resztkową sumę kwadratów dopasowania wykładniczego?

Mam następujące dane i chciałbym dopasować do niego model ujemnego wzrostu wykładniczego:

Days <- c( 1,5,12,16,22,27,36,43)
Emissions <- c( 936.76, 1458.68, 1787.23, 1840.04, 1928.97, 1963.63, 1965.37, 1985.71)
plot(Days, Emissions)
fit <- nls(Emissions ~ a* (1-exp(-b*Days)), start = list(a = 2000, b = 0.55))
curve((y = 1882 * (1 - exp(-0.5108*x))), from = 0, to =45, add = T, col = "green", lwd = 4)

Kod działa i drukowana jest linia dopasowania. Jednak dopasowanie nie jest wizualnie idealne, a resztkowa suma kwadratów wydaje się być dość duża (147073).

Jak możemy poprawić nasze dopasowanie? Czy dane w ogóle pozwalają na lepsze dopasowanie?

Nie mogliśmy znaleźć rozwiązania tego problemu w sieci. Każda bezpośrednia pomoc lub link do innych stron / postów jest bardzo mile widziana.

r nonlinear-regression fitting nls

— Strohmi
źródło

W takim przypadku, jeśli weźmie się pod uwagę model regresji

, gdzie

, to otrzymujesz podobne estymatory. Wykreślając obszary ufności, można zaobserwować, w jaki sposób wartości te są zawarte w regionach zaufania. Nie możesz oczekiwać idealnego dopasowania, chyba że interpolujesz punkty lub zastosujesz bardziej elastyczny model nieliniowy.

{Emissions}_{i} = f ({Days}_{i}, a, b) + ϵ_{i}

$\text{Emissions}_i=f(\text{Days}_i,a,b)+\epsilon_i$

ϵ_{i} \sim N (0, σ)

$\epsilon_i\sim N(0,\sigma)$

Zmieniłem tytuł, ponieważ „ujemny model wykładniczy” oznacza coś innego niż opisano w pytaniu.

— whuber

Dziękujemy za wyjaśnienie pytania (@whuber) i dziękuję za odpowiedź (@Prastrastinator). Jak obliczyć i wykreślić regiony ufności. A jaki byłby bardziej elastyczny model nieliniowy?

— Strohmi,

Potrzebujesz dodatkowego parametru. Zobacz, co się stanie

fit <- nls(Emissions ~ a* (1- u*exp(-b*Days)), start = list(a = 2000, b = 0.1, u=.5));  beta <- coefficients(fit); curve((y = beta["a"] * (1 - beta["u"] * exp(-beta["b"]*x))), add = T)

— whuber

@ whuber - może powinieneś zamieścić to jako odpowiedź?

— łucznik

Prawo (ujemne) wykładnicze ma postać . Jeśli jednak zezwalasz na zmiany jednostek w wartościach i , powiedzmy, że i , wówczas prawo zostanie wyrażone jako $y=-\exp(-x)$ $x$ $y$ $y = \alpha y' + \beta$ $x = \gamma x' + \delta$

α y^{'} + β = y = - \exp (- x) = - \exp (- γ x^{'} - δ),

$\alpha y' + \beta = y = -\exp(-x) = -\exp(-\gamma x' - \delta),$

który algebraicznie jest równoważny

y^{'} = \frac{- 1}{α} \exp (- γ x^{'} - δ) - β = a (1 - u \exp (- b x^{'}))

$y' = \frac{-1}{\alpha} \exp(-\gamma x' - \delta) - \beta = a\left(1 - u\exp(-b x')\right)$

na podstawie trzech parametrów , , oraz . Można rozpoznać jako parametr skali dla , jako parametr skali dla i , jak pochodzący z lokalizacji parametr . $a = -\beta/\alpha$ $u = 1/(\beta\exp(\delta))$ $b = \gamma$ $a$ $y$ $b$ $x$ $u$ $x$

Zasadniczo parametry te można zidentyfikować na podstawie wykresu :

Parametr jest wartością poziomej asymptoty, nieco mniejszą niż . $a$ $2000$
Parametr jest względną wartością, jaką krzywa podnosi od początku do poziomej asymptoty. Tutaj wzrost jest zatem nieco mniejszy niż ; względnie, to około asymptoty. $u$ $2000 - 937$ $0.55$
Ponieważ , gdy równa się trzykrotności wartości krzywa powinna wzrosnąć do około lub jego sumy. wzrostu z do prawie stawia nas około ; skanowanie całego wykresu wskazuje, że zajęło to od do dni. Wezwanie Chodźmy go dla uproszczenia, skąd $\exp(-3) \approx 0.05$ $x$ $1/b$ $1-0.05$ $95\%$ $95\%$ $937$ $2000$ $1950$ $20$ $25$ $24$ . (Ta metoda gałki ocznej skali wykładniczej jest standardem w niektórych dziedzinach, w których często stosuje się wykresy wykładnicze). $b \approx 3/24 = 0.125$ $95\%$

Zobaczmy, jak to wygląda:

plot(Days, Emissions)
curve((y = 2000 * (1 - 0.56 * exp(-0.125*x))), add = T)

Eyeball fit

Nieźle jak na początek! (Nawet pomimo pisania 0.56zamiast 0.55, co i tak było dość przybliżone.) Możemy go wypolerować za pomocą nls:

fit <- nls(Emissions ~ a * (1- u * exp(-b*Days)), start=list(a=2000, b=1/8, u=0.55))
beta <- coefficients(fit)
plot(Days, Emissions)
curve((y = beta["a"] * (1 - beta["u"] * exp(-beta["b"]*x))), add = T, col="Green", lwd=2)

NLS fit

Dane wyjściowe nlszawierają obszerne informacje na temat niepewności parametrów. Np. Prosty summarypodaje standardowe błędy szacunków:

> summary(fit)

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
a 1.969e+03  1.317e+01  149.51 2.54e-10 ***
b 1.603e-01  1.022e-02   15.69 1.91e-05 ***
u 6.091e-01  1.613e-02   37.75 2.46e-07 ***

Możemy czytać i pracować z całą macierzą kowariancji oszacowań, co jest przydatne do szacowania równoczesnych przedziałów ufności (przynajmniej dla dużych zestawów danych):

> vcov(fit)
             a             b             u
a 173.38613624 -8.720531e-02 -2.602935e-02
b  -0.08720531  1.044004e-04  9.442374e-05
u  -0.02602935  9.442374e-05  2.603217e-04

nls obsługuje wykresy profili dla parametrów, podając bardziej szczegółowe informacje o ich niepewności:

> plot(profile(fit))

$a$

Profile plot

$2$ $1945$ $1995$

— Whuber
źródło

res <- residuals(fit); res %*% res

u

$u$

2724

$2724$

147073

$147073$

Wszystko dobrze i dobry whuber. Ale może OP miał jakiś powód, aby wybrać model wykładniczy (a może tylko dlatego, że jest dobrze znany). Myślę, że najpierw należy przyjrzeć się resztom dla modelu wykładniczego. Wykreśl je względem potencjalnych zmiennych towarzyszących, aby zobaczyć, czy jest tam struktura, a nie tylko duży przypadkowy szum. Zanim przejdziesz do bardziej wyrafinowanych modeli, spróbuj sprawdzić, czy bardziej wymyślny model może pomóc.

— Michael R. Chernick

Dlaczego nie spojrzysz na oryginalną fabułę, Michael? Dzięki temu stanie się oczywiste, dlaczego potrzebny jest co najmniej jeden dodatkowy parametr. Należy również pamiętać, że w komentarzu do pytania OP zadał: „jaki byłby bardziej elastyczny model nieliniowy?” Jednym z implikacji wstępnej analizy przedstawionej w tej odpowiedzi jest to, że należy uznać za niezwykłe dopasowanie wykładnika o mniej niż trzech parametrach: w takich przypadkach musi istnieć pewne nieodłączne ograniczenie (takie jak wewnętrznie określona jednostka miary lub wewnętrzna lokalizacja dla

x

$x$

Nie krytykowałem twojej odpowiedzi! Nie widziałem żadnych fabuł. Wszystko, co zasugerowałem, to to, że wykresy reszt względem potencjalnych zmiennych towarzyszących powinny być pierwszym krokiem do znalezienia lepszego modelu. Gdybym pomyślał, że mam tam odpowiedź, dałbym odpowiedź, a nie podniosłem swój punkt widzenia jako stały. Myślałem, że dałeś świetną odpowiedź i byłem jednym z tych, którzy dali ci +1.

— Michael R. Chernick