GMM normal distribution with R

### blog: http://reakkt.com
# based on: http://cran.r-project.org/web/packages/gmm/vignettes/gmm_with_R.pdf

library(gmm)

n <- 100

x.mean  <- runif(1,-10,10)
x.sd    <- runif(1,1,10)

x <- rnorm(n,x.mean,x.sd)

plot(density(x))

### MLE

objFun <- function(p) {

  err <- 0

  if (p[2]>0)

    err <- -sum(dnorm(x,p[1],p[2],log=TRUE))

  return(err)

}

### GMM normal distribution

g2 <- function(tet,x) {

  m1 <- (tet[1]-x)
  m2 <- (tet[2]^2 - (x - tet[1])^2)

  f <- cbind(m1,m2)

  return(f)

}

g3 <- function(tet,x) {

  m1 <- (tet[1]-x)
  m2 <- (tet[2]^2 - (x - tet[1])^2)
  m3 <- x^3-tet[1]^3-3*tet[1]*tet[2]^2 # x^3-tet[1]*(tet[1]^2+3*tet[2]^2)

  f <- cbind(m1,m2,m3)

  return(f)

}

g4 <- function(tet,x) {
# normal distribution moments: http://en.wikipedia.org/wiki/Normal_distribution#Moments

  m1 <- (tet[1]-x)
  m2 <- (tet[2]^2 - (x - tet[1])^2)
  m3 <- x^3-(tet[1]^3+3*tet[1]*tet[2]^2)
  m4 <- x^4-(tet[1]^4+6*tet[1]^2*tet[2]^2+3*tet[2]^4)

  f <- cbind(m1,m2,m3,m4)

  return(f)

}

###

x.mean; mean(x)
x.sd; sd(x)

o <- optim(c(0,1),objFun) #method="SANN",control=list(maxit=50000))

o$par

t0 <- c(mean(x),sd(x))

x.gmm <- gmm(g2,x,t0) # 2 momenty
coef(x.gmm)

x.gmm <- gmm(g3,x,t0) # 3 momenty
coef(x.gmm)

x.gmm <- gmm(g4,x,t0) # 3 momenty
coef(x.gmm)