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

14

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.

Strohmi
źródło
1
W takim przypadku, jeśli weźmie się pod uwagę model regresji , gdzie ϵ iN ( 0 , σ ) , 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. Emissionsi=f(Daysi,a,b)+ϵiϵiN(0,σ)
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,
4
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
1
@ whuber - może powinieneś zamieścić to jako odpowiedź?
łucznik

Odpowiedzi:

16

Prawo (ujemne) wykładnicze ma postać . Jeśli jednak zezwalasz na zmiany jednostek w wartościach x i y , powiedzmy, że y = α y + β i x = γ x + δ , wówczas prawo zostanie wyrażone jakoy=exp(x)xyy=αy+βx=γx+δ

αy+β=y=exp(x)=exp(γxδ),

który algebraicznie jest równoważny

y=1αexp(γxδ)β=a(1uexp(bx))

na podstawie trzech parametrów = - β / α , u = 1 / ( β exp ( δ ) ) , oraz b = γ . Można rozpoznać jako parametr skali dla Y , B jako parametr skali dla X i U , jak pochodzący z lokalizacji parametr x .a=β/αu=1/(βexp(δ))b=γaybxux

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

  • Parametr jest wartością poziomej asymptoty, nieco mniejszą niż 2000 .a2000

  • Parametr jest względną wartością, jaką krzywa podnosi od początku do poziomej asymptoty. Tutaj wzrost jest zatem nieco mniejszy niż 2000 - 937 ; względnie, to około 0,55 asymptoty.u20009370.55

  • Ponieważ , gdy x równa się trzykrotności wartości 1 / b, krzywa powinna wzrosnąć do około 1 - 0,05 lub 95 % jego sumy. 95 % wzrostu z 937 do prawie 2000 stawia nas około 1950 r . ; skanowanie całego wykresu wskazuje, że zajęło to od 20 do 25 dni. Wezwanie Chodźmy go 24 dla uproszczenia, skąd b 3 / 24exp(3)0.05x1/b10.0595%95%93720001950202524 . (Ta 95 % metoda gałki ocznej skali wykładniczej jest standardem w niektórych dziedzinach, w których często stosuje się wykresy wykładnicze).b3/24=0.12595%

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

219451995

Whuber
źródło
res <- residuals(fit); res %*% resu2724147073
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
3
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 dlax
2
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