计算标准正态分布的风险函数

Question

基于Mathematica UUPDE database中给出的公式我已经在 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

风险函数是密度函数除以幸存者函数。您的代码的问题在于，您是从表面上接受这个定义并进行简单的除法运算；当分子和分母都是非常小的值（大约 1e-300）时，这种情况发生在分布的尾部，此操作在数值上变得不稳定。对于此类问题，更合适的解决方案是计算分子和分母（中等大小的负数而不是小数）的对数，从中减去对数分母对数分子，然后求幂。

R 提供了进行此计算所需的所有部分。您可以通过 pnorm(x,lower=FALSE) 获取幸存者功能；您可以分别在 dnorm() 和 pnorm() 中使用 log=TRUE 和 log.p=TRUE 来获得对数尺度上的密度和幸存者函数。所以：

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 可以在上面替换）。

计算标准正态分布的风险函数

calculating hazard function for the standard normal distribution

r

normal-distribution

hazard