# Computing an integral for expectation over a normal distribution

**URL:** <https://discourse.julialang.org/t/computing-an-integral-for-expectation-over-a-normal-distribution/53121>\
**Category:** New to Julia\
**Tags:** numerics\
**Created:** [January 10, 2021, 4:58pm UTC](https://discourse.julialang.org/t/computing-an-integral-for-expectation-over-a-normal-distribution/53121 "2021-01-10T16:58:07Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![samerb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/samerb/32/9242_2.png) [@samerb](https://discourse.julialang.org/u/samerb)\
**Post date:** [January 10, 2021, 4:58pm UTC](https://discourse.julialang.org/t/computing-an-integral-for-expectation-over-a-normal-distribution/53121/1 "2021-01-10T16:58:07Z")

</div>

Solved: The main problem was that Normal accepts Normal(μ,σ). I thought it was Normal(μ,σ^2). The other discrepancies are still a little mysterious but this explains why my Gauss Hermite was so off.

I want to compute the integral of a logit function with respect to a normal distribution. Specifically I want to compute: \int\_{-\infty}^{\infty} \frac{exp(5\beta)}{1+\exp(5\beta)} \cdot \exp( \frac{-(\beta - \mu)^2}{2\sigma^2}) d\beta with \mu=0.5, \sigma = 2.

I try three different approaches to computing the integral: 1) [GaussGK](https://github.com/JuliaMath/QuadGK.jl), 2) Monte Carlo, 3) [Gauss-Hermite](https://en.wikipedia.org/wiki/Gauss%E2%80%93Hermite_quadrature) with [FastGaussQuadrature package](https://github.com/JuliaApproximation/FastGaussQuadrature.jl#the-algorithm-for-gauss-hermite).

For the first two approaches I get the same answer up to 3 digits but the third approach (Gauss-Hermite) is not close to the other two. **Question: why not**? I also would’ve thought the first two approaches would have been a little closer, but that’s not my main concern.

First I write a function for the logit part (to deal with big numbers problems) using softmax from StatsFuns package and check that it’s correct.

```julia-auto
function f(β,x)
   exp(β'*x)/(1 + exp(β'*x)) 
end

function logi(β,x=5)
    softmax([β'*x,0])[1]
end

# check to make sure logi == f
x = 5*rand(10)
β = 5*rand(10)

f(β,x) == logi(β,x) # -> true

```

## 1) QuadGK package:

```julia-auto
using StatsFuns
using Distributions
d = Normal(0.5,2)
function bilog(β, x=5, dist=Normal(0.5,2))
   return logi(β,x) * pdf(dist, β)
end

using QuadGK
tol = 10^-14
quant_tol_l, quant_tol_u = (quantile(d,tol), quantile(d,1-tol) ) 
# bounds chosen so that not much is left over 

integral, err = quadgk(bilog, quant_tol_l, quant_tol_u, rtol=tol) 
# -> (0.5971674, 5.3732205e-15)

```

## 2) Monte Carlo approach

```julia-auto
nodes = quantile(d,rand(1000000))
mean(logi.(nodes, 5)) # -> 0.5975626 
# would have thought 0.5975626 would match 0.597167 from 
# above a little more closely but not too bad
# I have run this a few to times since it is random so this is not a "special" draw.

```

## 3) Gauss-Hermite approach.

Here’s where the real problem is. I try to use Gauss-Hermite quadrature which is described [here](https://en.wikipedia.org/wiki/Gauss%E2%80%93Hermite_quadrature). I follow the formula for the change of variables as described there but I am getting a significantly different answer from above. Any idea why?

```julia-auto
using FastGaussQuadrature
nodes, weights = gausshermite( 1000000 )
sum(weights) # -> 1.772 should this be 1?
vals = (β -> logi(β*sqrt(2)*sqrt(2) + 0.5)).(nodes) 
(1/sqrt(π))*vals'*weights # -> 0.63406 not equal .597

```

---

<div class="post-metadata">

**Author:** ![jlperla](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jlperla/32/34332_2.png) [@jlperla](https://discourse.julialang.org/u/jlperla)\
**Post date:** [January 10, 2021, 5:03pm UTC](https://discourse.julialang.org/t/computing-an-integral-for-expectation-over-a-normal-distribution/53121/2 "2021-01-10T17:03:12Z")

</div>

See if [GitHub - QuantEcon/Expectations.jl: Expectation operators for Distributions.jl objects](https://github.com/QuantEcon/Expectations.jl) works. It hopefully does the change of variables correctly, so you could use it or copy the logic

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 10, 2021, 5:34pm UTC](https://discourse.julialang.org/t/computing-an-integral-for-expectation-over-a-normal-distribution/53121/3 "2021-01-10T17:34:27Z")

</div>

> [@samerb](#):
>
> `quadgk(bilog, quant_tol_l, quant_tol_u, rtol=tol) `

(Note that you can just integrate from `-Inf` to `+Inf`, and it will give you about the same accuracy and efficiency.)

---

<div class="post-metadata">

**Author:** ![tbeason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tbeason/32/15898_2.png) [@tbeason](https://discourse.julialang.org/u/tbeason)\
**Post date:** [January 10, 2021, 5:40pm UTC](https://discourse.julialang.org/t/computing-an-integral-for-expectation-over-a-normal-distribution/53121/4 "2021-01-10T17:40:52Z")

</div>

The weights sum to 1.772 because they are not normalized – that is what the `1/sqrt(pi)` is doing.

Not sure why the result is different.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 10, 2021, 5:41pm UTC](https://discourse.julialang.org/t/computing-an-integral-for-expectation-over-a-normal-distribution/53121/5 "2021-01-10T17:41:06Z")

</div>

> [@samerb](#):
>
> `vals = (β -> logi(β*sqrt(2)*sqrt(2) + 0.5)).(nodes) `

Looks like you forgot the Jacobian factor from the change of variables?

---

<div class="post-metadata">

**Author:** ![samerb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/samerb/32/9242_2.png) [@samerb](https://discourse.julialang.org/u/samerb)\
**Post date:** [January 10, 2021, 5:50pm UTC](https://discourse.julialang.org/t/computing-an-integral-for-expectation-over-a-normal-distribution/53121/6 "2021-01-10T17:50:21Z")

</div>

The result of QuantEcon Expectations is:

```julia
using Expectations
dist = Normal(0.5, 2)
E = expectation(dist)
E(β -> logi(β))

```

which gives 0.58975 as opposed to 0.5971674 which is what I got from the first approach. confirms there something wrong with my Gauss-Hermite.

---

<div class="post-metadata">

**Author:** ![samerb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/samerb/32/9242_2.png) [@samerb](https://discourse.julialang.org/u/samerb)\
**Post date:** [January 10, 2021, 5:57pm UTC](https://discourse.julialang.org/t/computing-an-integral-for-expectation-over-a-normal-distribution/53121/7 "2021-01-10T17:57:50Z")

</div>

Ok I see (at least) the main problem! Normal accepts Normal(μ,σ). I thought it was Normal(μ,σ^2). IMHO that’s confusing since N(μ,σ^2) is the standard in my experience but what I know. Other differences seem to be from different number of nodes used (though monte carlo still seems to be a bit off).

---

<div class="post-metadata">

**Author:** ![jlperla](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jlperla/32/34332_2.png) [@jlperla](https://discourse.julialang.org/u/jlperla)\
**Post date:** [January 10, 2021, 6:17pm UTC](https://discourse.julialang.org/t/computing-an-integral-for-expectation-over-a-normal-distribution/53121/8 "2021-01-10T18:17:22Z")

</div>

Do you use the same number of nodes to compare?

---

<div class="post-metadata">

**Author:** ![tbeason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tbeason/32/15898_2.png) [@tbeason](https://discourse.julialang.org/u/tbeason)\
**Post date:** [January 10, 2021, 6:20pm UTC](https://discourse.julialang.org/t/computing-an-integral-for-expectation-over-a-normal-distribution/53121/9 "2021-01-10T18:20:50Z")

</div>

No, you can see he used 1m nodes for bare GH and the default for Expectations, which is like 30?