Replace approximate M values with analyically obtained ones

[oguzhanogreden]

I'm following the notation of Woods & Lin (2009). M matrix is defined in page 105.

Right now this matrix is filled with approximate values. I think can obtain the exact values using moment generating functions and a bit of manual labor.

[philchalmers]

These values look fine to me, they should be whole numbers under the Z ~ N(0,1) case. In what sense are they approximations?

[oguzhanogreden]

I sampled a vector obs where each element was a realization of Z, then obtained mean(obs^R). I repeated this to estimate E[Z^R].

[philchalmers]

That's one way to do it without the MGF. I just verified your moments using integrate(), like so:

> f <- function(Z) Z * dnorm(Z, mean=0, sd=1)
> integrate(f, -Inf, Inf)
0 with absolute error < 0
> 
> f <- function(Z) Z^2 * dnorm(Z, mean=0, sd=1)
> integrate(f, -Inf, Inf)
1 with absolute error < 1.2e-07
> 
> f <- function(Z) Z^3 * dnorm(Z, mean=0, sd=1)
> integrate(f, -Inf, Inf)
0 with absolute error < 0
> 
> f <- function(Z) Z^5 * dnorm(Z, mean=0, sd=1)
> integrate(f, -Inf, Inf)
0 with absolute error < 0
> 
> f <- function(Z) Z^6 * dnorm(Z, mean=0, sd=1)
> integrate(f, -Inf, Inf)
15 with absolute error < 7.9e-05
> 
> f <- function(Z) Z^8 * dnorm(Z, mean=0, sd=1)
> integrate(f, -Inf, Inf)
105 with absolute error < 0.0014
> 
> f <- function(Z) Z^9 * dnorm(Z, mean=0, sd=1)
> integrate(f, -Inf, Inf)
0 with absolute error < 0

This is what I was referring to before when computing the expected values beforehand for Z, which would be fairly cheap given the outer product form E(U * t(U) ), but computing the moments analytically would be more efficient.

[oguzhanogreden]

A-ha, simple but effective. Didn't think of this myself. I'll work the analytical forms out this week.

[oguzhanogreden]

Turns out I got most correct. Only E[X^18] and E[X^20] needed adjustment.