Misalkan saya memiliki fungsi seperti:
f <- function(x){
exp(x) / (1 + exp(x))
}
itu seharusnya bekerja untuk setiap nilai nyata x, tetapi sebenarnya mengembalikan NaN ketika x adalah 710 atau lebih besar. Saya bertanya-tanya apa cara yang tepat untuk menangani masalah ini. Saya menyadari mudah untuk membuatnya hanya mengembalikan 1, tetapi mungkin itu bukan perilaku yang baik dari sudut pandang seorang ahli statistik. Apakah ada yang punya komentar atau saran?
Jawaban:
Dalam hal ini
NaN
(bukan angka) dikembalikan karena perhitungan luapan eksponensial dalam aritmatika presisi ganda.Ekspresi yang setara secara aljabar, diperluas dalam seri MacLaurin sekitar , adalah0
Karena ini adalah seri bergantian, kesalahan yang dibuat dalam menjatuhkan istilah apa pun tidak lebih besar dari ukuran istilah berikutnya. Jadi ketika , kesalahannya tidak lebih besar dari relatif terhadap nilai sebenarnya. Itu jauh lebih tepat daripada perhitungan statistik apa pun yang perlu dilakukan, jadi Anda boleh mengganti nilai pengembalian dengan dalam situasi ini.x > 710 exp( - 710 ) ≈10- 308≈2- 1024 1
Menariknya,
R
tidak akan menghasilkanNaN
ketika eksponensial mengalir . Dengan demikian Anda bisa memilih versi perhitungan yang lebih andal, tergantung pada tandax
, seperti padaMasalah ini muncul di hampir semua platform komputasi (saya belum melihat pengecualian) dan mereka akan bervariasi dalam bagaimana mereka menangani overflow dan underflow. Eksponensial terkenal karena menciptakan masalah seperti ini, tetapi mereka tidak sendirian. Oleh karena itu, tidak cukup hanya dengan memiliki solusi
R
: seorang ahli statistik yang baik memahami prinsip-prinsip aritmatika komputer dan tahu bagaimana menggunakannya untuk mendeteksi dan mengatasi kekhasan lingkungan komputernya.sumber
Yang lain sudah membahas masalah komputasi, jadi saya akan menyerahkannya kepada mereka. Karena saya menganggap Anda bekerja dengan R, saya pikir saya akan menunjukkan paket boot dilengkapi dengan fungsi logit terbalik sendiri untuk Anda gunakan yang cukup stabil secara komputasi:
require(boot) inv.logit(710)
tampaknya mengevaluasi ke 1 seperti yang diinginkan.
sumber
plogis(710)
dapatkan hasil yang sama. (Memanginv.logit
hanya alias untukplogis
.)