Calculation of hazard function in R for standard normal distribution

Based on the formulas given in the UUPDE Mathematica database http://blog.wolfram.com/2013/02/01/the-ultimate-univariate-probability-distribution-explorer/ I built a standard normal distribution of the hazard function in R.

In a certain range it seems correct, numerical problems arise with large values, see the attached figure. Below is the full R code.

Any comments would be greatly appreciated. See FigurePDF, CDF, HF and SF plots for standard normal distribution

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)
+4
source share
1 answer

- , . , ; , ( 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)

enter image description here

, ( foo R dfoo CDF pfoo, ).

+9

Source: https://habr.com/ru/post/1654689/


All Articles