Hola Freddy el resultado no le da es porque usted esta utilizando, la distribución uniforme U ~ [0,1] y para intentar calcula una integral que se mueve entre -infinito e infinito como es el valor esperado de la distribución normal, por lo tanto no use la uniforme, utilice distribuciones continuas que tengan ese soporte como la normal estándar o la t-student.
Por ejemplo si quiero calcular la siguiente integral, de ∫ 5/(1+x^2) con x ϵ (-∞, ∞), utilizando el método de montecarlo. Entonces sea h(x) = 5/(1+x^2) y mi propuesta f(x) = N(0,1) la dist. normal estándar.
Integral_impropia <- function(h){
h.aste <- function(x){h(x)*exp(0.5*x^2)}
n <- 2500000
set.seed(94)
u <- rnorm(n)
res <- sqrt(2*pi)*mean(h.aste(u))
return(res)
}
h <- function(x){5/(1+x^2)}
integrate(h,-Inf,Inf)
15.70796 with absolute error < 2.6e-09
Integral_impropia(h)
[1] 15.61579
Ahora si quiere calcular el valor esperado de la distribución normal solo cambie h en el código.
h2 <- function(x){1/sqrt(2*pi)*x*exp(-x^2/2)}
integrate(h2,-Inf,Inf)
0 with absolute error < 0
Integral_impropia(h2)
[1] -0.0004724279