Mengapa nls () memberi saya kesalahan "matriks gradien singular pada estimasi parameter awal"?

21

Saya memiliki beberapa data dasar tentang pengurangan emisi dan biaya per mobil:

q24 <- read.table(text = "reductions  cost.per.car
    50  45
    55  55
    60  62
    65  70
    70  80
    75  90
    80  100
    85  200
    90  375
    95  600
    ",header = TRUE, sep = "")

Saya tahu bahwa ini adalah fungsi eksponensial, jadi saya berharap dapat menemukan model yang cocok dengan:

    model <- nls(cost.per.car ~ a * exp(b * reductions) + c, 
         data = q24, 
         start = list(a=1, b=1, c=0))

tapi saya mendapatkan kesalahan:

Error in nlsModel(formula, mf, start, wts) : 
  singular gradient matrix at initial parameter estimates

Saya telah membaca banyak pertanyaan tentang kesalahan yang saya lihat dan saya temui bahwa masalahnya mungkin karena saya membutuhkan nilai yang lebih baik / berbeda start( initial parameter estimatesmasuk akal sedikit) tetapi saya tidak yakin, mengingat data yang saya miliki, bagaimana saya akan memperkirakan parameter yang lebih baik.

Amanda
sumber
Saya akan menyarankan memulai penguraian Anda dengan mencari situs kami untuk pesan kesalahan .
whuber
3
Sebenarnya, saya melakukan itu dan pencarian saya untuk kesalahan penuh muncul pertanyaan setengah dipanggang dengan tiga titik data dan tidak ada jawaban. Tetapi pencarian Anda yang lebih spesifik memang mendapatkan beberapa hasil. Mungkin karena Anda memiliki lebih banyak pengalaman di sini dan tahu istilah mana yang paling relevan.
Amanda
Satu hal yang saya temukan tentang kesalahan perangkat lunak adalah bahwa pencarian untuk pesan kesalahan tertentu (biasanya dalam tanda kutip) adalah cara paling pasti untuk mengetahui apakah telah dibahas sebelumnya. (Ini berlaku di seluruh Internet, tidak hanya di situs SE.) Seperti yang dikatakan oleh pesan "ditahan" kami, jika penelitian tambahan Anda tidak menyelesaikan masalah Anda, silakan kembali dan dorong kami sedikit: pertanyaan ini ada di persimpangan statistik dan komputasi dan mungkin memaparkan beberapa masalah yang sangat menarik di sini.
whuber
1
Kesesuaian untuk nilai awal Anda sangat jauh dari data; bandingkan exp(50)dan exp(95)dengan nilai-y pada x = 50 dan x = 95. Jika Anda menetapkan c=0dan mengambil log y (membuat hubungan linier), Anda dapat menggunakan regresi untuk mendapatkan taksiran awal untuk log ( ) dan yang akan mencukupi untuk data Anda (atau jika Anda mencocokkan sebuah baris dengan sumbernya, Anda dapat meninggalkan at 1 dan cukup gunakan estimasi untuk ; itu juga cukup untuk data Anda). Jika jauh di luar interval yang cukup sempit di sekitar kedua nilai tersebut, Anda akan mengalami beberapa masalah. [Atau coba algoritma lain]b a b bababb
Glen_b -Reinstate Monica
1
Terima kasih @Glen_b. Saya berharap saya bisa menggunakan R sebagai pengganti kalkulator grafik untuk bekerja melalui buku teks statistik intro (dan melompati kursus itu sendiri) jadi saya mulai dengan hanya wawasan statistik paling sederhana, tetapi banyak pengalaman melakukan pengirisan dan pencobaan lain dalam R .
Amanda

Jawaban:

38

Secara otomatis menemukan nilai awal yang baik untuk model nonlinear adalah seni. (Ini relatif mudah untuk dataset satu kali ketika Anda hanya bisa memplot data dan membuat beberapa tebakan yang bagus secara visual.) Salah satu pendekatan adalah membuat linierisasi model dan menggunakan estimasi kuadrat terkecil.

Dalam hal ini, model memiliki bentuk

E(Y)=aexp(bx)+c

untuk parameter yang tidak diketahui . Kehadiran eksponensial mendorong kita untuk menggunakan logaritma - tetapi penambahan membuatnya sulit untuk melakukan itu. Perhatikan, meskipun, bahwa jika positif maka akan kurang dari nilai yang diharapkan terkecil dari --dan karena itu mungkin sedikit kurang dari yang terkecil diamati nilai . (Jika bisa negatif Anda juga harus mempertimbangkan nilai yang sedikit lebih besar dari nilai diamati terbesar .)c a c Y Y a c Ya,b,ccacYYacY

Mari kita jaga dengan menggunakan estimasi awal kira-kira setengah dari minimum pengamatan . Model sekarang dapat ditulis ulang tanpa istilah aditif berduri sebagaic 0 y icc0yi

E(Y)c0aexp(bx).

Bahwa kita dapat mengambil log dari:

log(E(Y)c0)log(a)+bx.

Itu adalah pendekatan linier ke model. Baik danlog(a)b

Berikut adalah kode yang direvisi:

c.0 <- min(q24$cost.per.car) * 0.5
model.0 <- lm(log(cost.per.car - c.0) ~ reductions, data=q24)
start <- list(a=exp(coef(model.0)[1]), b=coef(model.0)[2], c=c.0)
model <- nls(cost.per.car ~ a * exp(b * reductions) + c, data = q24, start = start)

Outputnya (untuk data contoh) adalah

Nonlinear regression model
  model: cost.per.car ~ a * exp(b * reductions) + c
   data: q24
        a         b         c 
 0.003289  0.126805 48.487386 
 residual sum-of-squares: 2243

Number of iterations to convergence: 38 
Achieved convergence tolerance: 1.374e-06

Konvergensi terlihat bagus. Mari kita plot itu:

plot(q24)
p <- coef(model)
curve(p["a"] * exp(p["b"] * x) + p["c"], lwd=2, col="Red", add=TRUE)

Angka

Itu bekerja dengan baik!

ya<0


Metode lain untuk memperkirakan nilai awal bergantung pada pemahaman apa artinya, yang dapat didasarkan pada pengalaman, teori fisik, dll. Contoh tambahan dari kesesuaian nonlinier (cukup sulit) yang nilai awalnya dapat ditentukan dengan cara ini dijelaskan dalam jawaban saya di /stats//a/15769 .

Analisis visual sebar (untuk menentukan perkiraan parameter awal) dijelaskan dan diilustrasikan di /stats//a/32832 .

Dalam beberapa keadaan, urutan cocok nonlinear dibuat di mana Anda dapat mengharapkan solusi berubah perlahan. Dalam hal ini, sering kali nyaman (dan cepat) untuk menggunakan solusi sebelumnya sebagai perkiraan awal untuk yang berikutnya . Saya ingat menggunakan teknik ini (tanpa komentar) di /stats//a/63169 .

whuber
sumber
2

Pustaka ini dapat menyelesaikan masalah saya dengan nls singular gradient: http://www.r-bloggers.com/a-better-nls/ Contoh:

library(minpack.lm)
nlsLM(function, start=list(variable=2,variable2=12))
joshring
sumber
Fungsi itu sepertinya dipanggil nls.lmsekarang.
Matt
-1

Jadi ... Saya rasa saya salah membaca ini sebagai fungsi eksponensial. Yang saya butuhkan adalahpoly()

model <- lm(cost.per.car ~ poly(reductions, 3), data=q24)
new.data <- data.frame(reductions = c(91,92,93,94))
predict(model, new.data)

plot(q24)
lines(q24$reductions, predict(model, list(reductions = q24$reductions)))

Atau, menggunakan lattice:

xyplot(cost.per.car ~ reductions, data = q24,
       panel = function(x, y) {
         panel.xyplot(x, y)
         panel.lines(x, predict(model,list(reductions = x) ))
       }, 
       xlab = "Reductions", 
       ylab = "Cost per car")
Amanda
sumber
2
Ini tidak menjawab pertanyaan yang Anda ajukan - itu mengubahnya menjadi sesuatu yang berbeda (dan agak kurang menarik, IMHO).
whuber
6
Meskipun, itu mungkin memecahkan masalah menyesuaikan fungsi untuk mewakili data, jawaban Anda yang diterima bukan harapan untuk pertanyaan Anda. Mr. @whuber memberikan Anda penjelasan yang bagus dan pantas menerima jawaban yang diterima.
Lourenco