вычисления функции риски в R для стандартного нормального распределения

Question

на основе формул, приведенных в базе данных Mathematica UUPDE http://blog.wolfram.com/2013/02/01/the-ultimate-univariate-probability-distribution-explorer/ Я нанесенная функции риски стандартного нормальное распределение в R.

кажется правильными в определенном диапазоне , числовые проблемы возникают для больших значений, см. прилагаемый рисунок. Ниже полный код R.

Любые комментарии были бы очень оценены. смотри рисунок

PDF = function(x) { 1/(sqrt(2*pi))*exp(-x^2/2) } 
erf <- function(x) 2 * pnorm(x * sqrt(2)) - 1 
erfc <- function(x) 2 * pnorm(x * sqrt(2), lower = FALSE) 
CDF = function(x) { 1/2 * (1 + erf(x/(sqrt(2)))) } 
HF = function(x) { sqrt(2/pi)/(exp(x^2/2)*(2-erfc(-x/sqrt(2)))) } 
SF = function(x) { 1 - 1/2 *erfc(-x/sqrt(2)) } 

par(mar=c(3,3,1.5,0.5), oma=c(0,0,0,0), mgp=c(2,1,0)) 
par(mfrow = c(2, 2)) 

x = seq(from = -4,to = 10,by = .001) 

##### PDF 
a = PDF(x) 
plot(x,a,'l',main='',ylab="PDF",xlab="x") 
grid(nx = NULL,ny = NULL,col = "grey",lty = "dotted",lwd = par("lwd"),equilogs = TRUE) 

##### CDF 
a = CDF(x) 
plot(x,a,'l',main='',ylab="CDF",xlab="x") 
grid(nx = NULL,ny = NULL,col = "grey",lty = "dotted",lwd = par("lwd"),equilogs = TRUE) 

##### HF 
a = HF(x) 
plot(x,a,'l',main='',ylab="HF",xlab="x") 
grid(nx = NULL,ny = NULL,col = "grey",lty = "dotted",lwd = par("lwd"),equilogs = TRUE) 

##### SF 
a = SF(x) 
plot(x,a,'l',main='',ylab="SF",xlab="x") 
grid(nx = NULL,ny = NULL,col = "grey",lty = "dotted",lwd = par("lwd"),equilogs = TRUE) 

Answer 1

Функция опасности является функция плотности, деленной на функции оставшегося в живых. Проблема с вашим кодом заключается в том, что вы принимаете это определение по номинальной стоимости и выполняете операцию простого деления; когда и числитель, и знаменатель очень маленькие (порядка 1е-300), что происходит в хвосте распределения, эта операция становится численно неустойчивой. Для такого рода проблем более подходящим решением является вычисление логарифмов числителя и знаменателя (которые являются отрицательными числами умеренного размера, а не крошечными числами), вычитают логарифмический знаменатель из логарифмического числителя и затем экспонентуют.

R содержит все части, необходимые для расчета. Вы можете получить функцию оставшегося в живых через pnorm(x,lower=FALSE); вы можете получить плотность и функции выживших на шкале журнала, используя log=TRUE и log.p=TRUE в dnorm() и pnorm() соответственно. Итак:

HF <- function(x) { 
    exp(dnorm(x,log=TRUE)-pnorm(x,lower=FALSE,log.p=TRUE)) 
} 
curve(HF,from=-4,to=10) 

Эта стратегия может быть обобщена для вычисления функции риски для любого распределения при условии, что лог-плотность и срубы выживших функций доступны (в общем случае для распределения foo R обеспечивает функцию плотности dfoo и CDF pfoo, которые могут быть заменены выше).