Differentiating pnorm() in R to show that the PDF of a continuous variable is the derivative of the-CodePudding

I need to demonstrate that the probability density function is the derivative of the CDF. Any distribution will do, but I have been trying with the normal. I have got as far as:

set.seed(53)
b <- rnorm(500)
db <- density(b)
plot(db)

Then I can calculate the cumulative probabilities using pnorm(b), but then I don't know how to differentiate, because D() requires an expression rather than pnorm(). Could anyone help, please?

CodePudding user response：

Here's the console scrape from where I demonstrated the near equality (to 5 or 7 decimal places) of the integral of dnorm to pnorm from -Inf to selected values of "x": The Fundamental Theorem of Calculus says that if the integral of a function f(x) is g(x) then f(x) is the derivative of g(x). (Or words to that effect.)

> sapply(c(0,Inf), function(x) integrate(dnorm, lower=-Inf, upper=x))
             [,1]         [,2]        
value        0.5          1           
abs.error    4.680562e-05 9.361124e-05
subdivisions 3            3           
message      "OK"         "OK"        
call         expression   expression  
> sapply(c(0,Inf), function(x) integrate(dnorm, lower=-Inf, upper=x)$value)
[1] 0.5 1.0
> sapply(seq(-3,3, by=0.5), function(x) integrate(dnorm, lower=-Inf, upper=x)$value)
 [1] 0.001349899 0.006209665 0.022750132 0.066807201 0.158655254 0.308537539
 [7] 0.500000000 0.691462461 0.841344751 0.933192799 0.977249868 0.993790335
[13] 0.998650102
> pnorm(seq(-3,3, by=0.5)
  )
 [1] 0.001349898 0.006209665 0.022750132 0.066807201 0.158655254 0.308537539
 [7] 0.500000000 0.691462461 0.841344746 0.933192799 0.977249868 0.993790335
[13] 0.998650102

I wasn't sure that the D() was "smart" enough to to the symbolic differentiation, but I shouldn't have been so skeptical. This bit of console interaction was done following the examples on the ?deriv help page:

> D(quote(pnorm(x)), "x")
dnorm(x)

Also ... here's something you can get with deriv:

> norm.expr <- expression(pnorm(x))
> deriv(norm.expr, "x")
expression({
    .value <- pnorm(x)
    .grad <- array(0, c(length(.value), 1L), list(NULL, c("x")))
    .grad[, "x"] <- dnorm(x)
    attr(.value, "gradient") <- .grad
    .value
})