Bagaimana cara meminimalkan jumlah residu kuadrat dari fit eksponensial?

14

Saya memiliki data berikut dan ingin mencocokkan model pertumbuhan eksponensial negatif dengan itu:

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)

Kode berfungsi dan garis pas diplot. Namun, kecocokan secara visual tidak ideal, dan jumlah kuadrat residu tampaknya cukup besar (147073).

Bagaimana kita dapat meningkatkan kecocokan kita? Apakah data memungkinkan kesesuaian yang lebih baik?

Kami tidak dapat menemukan solusi untuk tantangan ini di internet. Bantuan atau tautan langsung ke situs web / posting lain sangat dihargai.

Strohmi
sumber
1
Dalam hal ini, jika Anda mempertimbangkan model regresi , di mana ϵ iN ( 0 , σ ) , maka Anda mendapatkan penduga yang serupa. Dengan memplot wilayah kepercayaan, seseorang dapat mengamati bagaimana nilai-nilai ini terkandung dalam wilayah kepercayaan. Anda tidak dapat mengharapkan pasangan yang sempurna kecuali jika Anda menginterpolasi poin atau menggunakan model nonlinear yang lebih fleksibel. Emisisaya=f(Berhari-harisaya,Sebuah,b)+ϵsayaϵsayaN(0,σ)
Saya mengubah judul karena "model eksponensial negatif" berarti sesuatu yang berbeda dari yang dijelaskan dalam pertanyaan.
whuber
Terima kasih telah membuat pertanyaan menjadi lebih jelas (@whuber) dan terima kasih atas jawaban Anda (@Procrastinator). Bagaimana saya bisa menghitung dan memplot wilayah kepercayaan. Dan, apa yang akan menjadi model nonlinear yang lebih fleksibel?
Strohmi
4
Anda memerlukan parameter tambahan. Lihat apa yang terjadi 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 - mungkin Anda harus memposting itu sebagai jawaban?
jbowman

Jawaban:

16

Hukum eksponensial (negatif) berbentuk . Ketika Anda mengizinkan perubahan unit dalam nilai x dan y , katakanlah ke y = α y + β dan x = γ x + δ , maka hukum akan dinyatakan sebagaiy=-exp(-x)xyy=αy+βx=γx+δ

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

yang secara aljabar setara dengan

y=-1αexp(-γx-δ)-β=Sebuah(1-kamuexp(-bx))

menggunakan tiga parameter , u = 1 / ( β exp ( δ ) ) , dan b = γ . Kita dapat mengenali a sebagai parameter skala untuk y , b sebagai parameter skala untuk x , dan u berasal dari parameter lokasi untuk x .Sebuah=-β/αkamu=1/(βexp(δ))b=γSebuahybxkamux

Sebagai patokan, parameter ini dapat diidentifikasi sekilas dari plot :

  • Parameter adalah nilai asymptote horizontal, sedikit kurang dari 2000 .Sebuah2000

  • Parameter adalah jumlah relatif kurva naik dari asal ke asimtot horizontal. Di sini, kenaikan itu sedikit kurang dari 2000 - 937 ; relatif, itu sekitar 0,55 asimtot.kamu2000-9370,55

  • Karena , ketika x sama dengan tiga kali nilai 1 / b kurva seharusnya naik menjadi sekitar 1 - 0,05 atau 95 % dari totalnya. 95 % kenaikan dari 937 menjadi hampir 2000 menempatkan kita sekitar 1950 ; pemindaian di plot menunjukkan ini membutuhkan waktu 20 hingga 25 hari. Mari kita menyebutnya 24 untuk kesederhanaan, mana b 3 / 24exp(-3)0,05x1/b10,0595%95%93720001950202524 . (Metode 95 % untuk mengamati skala eksponensial ini merupakan standar di beberapa bidang yang sering menggunakan plot eksponensial.)b3/24=0,12595%

Mari kita lihat seperti apa ini:

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

Bola mata cocok

Tidak buruk untuk permulaan! (Bahkan meskipun mengetik 0.56di tempat 0.55., Yang merupakan perkiraan kasar pula) Kami bisa memolesnya dengan 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 cocok

Output nlsberisi informasi luas tentang ketidakpastian parameter. Misalnya , sederhana summarymemberikan kesalahan estimasi standar:

> 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 ***

Kita dapat membaca dan bekerja dengan seluruh matriks kovarian estimasi, yang berguna untuk memperkirakan interval kepercayaan simultan (setidaknya untuk dataset besar):

> 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 mendukung plot profil untuk parameter, memberikan informasi lebih rinci tentang ketidakpastiannya:

> plot(profile(fit))

a

Plot profil

219451995

whuber
sumber
res <- residuals(fit); res %*% resu2724147073
Semuanya baik dan baik. Tapi mungkin OP punya alasan untuk memilih model eksponensial (atau mungkin hanya karena itu sudah terkenal). Saya pikir pertama residual harus dilihat untuk model eksponensial. Plotkan mereka terhadap kovariat potensial untuk melihat apakah ada struktur di sana dan bukan hanya noise acak yang besar. Sebelum beralih ke model yang lebih canggih, cobalah untuk melihat apakah model yang lebih bagus mungkin bisa membantu.
Michael R. Chernick
3
x
2
Saya tidak mengkritik jawaban Anda! Saya tidak melihat plot residual. Yang saya sarankan adalah bahwa plot residu vs kovariat potensial harus menjadi langkah pertama dalam menemukan model yang lebih baik. Jika saya pikir saya punya jawaban untuk diletakkan di sana saya akan memberikan jawaban daripada mengangkat poin saya sebagai konstan. Saya pikir Anda memberikan respons yang bagus dan saya termasuk di antara mereka yang memberi Anda +1.
Michael R. Chernick