Using the mills ratio result, let X∼N(μ,σ2), then
However, when calculating it in R. I don't obtain the correct results as
> mu <- 1 > sigma <- 2 > a <- 3 > x <- rnorm(1000000, mu, sigma) > x <- x[x < a] > mean(x)  0.4254786 > > mu - sigma * dnorm(a, mu, sigma) / pnorm(a, mu, sigma)  0.7124
What am I doing wrong?
Your formula implementation is wrong because,
As you can see, we have an extra σ in the denominator of fX,μ,σ(x), which yields:
dnorm method gives you fX,μ,σ(x), where you need to multiply it with σ to obtain ϕ. Since your σ=2, this can be practically done via subtracting the second term again, which is 1−0.7124=0.2876:
which is close to your estimate.