mjaniec

GMM normal distribution with R

Aug 28th, 2012
252
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. ### blog: http://reakkt.com
  2. # based on: http://cran.r-project.org/web/packages/gmm/vignettes/gmm_with_R.pdf
  3.  
  4. library(gmm)
  5.  
  6. n <- 100
  7.  
  8. x.mean  <- runif(1,-10,10)
  9. x.sd    <- runif(1,1,10)
  10.  
  11. x <- rnorm(n,x.mean,x.sd)
  12.  
  13. plot(density(x))
  14.  
  15. ### MLE
  16.  
  17. objFun <- function(p) {
  18.  
  19.   err <- 0
  20.  
  21.   if (p[2]>0)
  22.    
  23.     err <- -sum(dnorm(x,p[1],p[2],log=TRUE))
  24.  
  25.   return(err)
  26.  
  27. }
  28.  
  29. ### GMM normal distribution
  30.  
  31. g2 <- function(tet,x) {
  32.  
  33.   m1 <- (tet[1]-x)
  34.   m2 <- (tet[2]^2 - (x - tet[1])^2)
  35.  
  36.   f <- cbind(m1,m2)
  37.  
  38.   return(f)
  39.  
  40. }
  41.  
  42. g3 <- function(tet,x) {
  43.  
  44.   m1 <- (tet[1]-x)
  45.   m2 <- (tet[2]^2 - (x - tet[1])^2)
  46.   m3 <- x^3-tet[1]^3-3*tet[1]*tet[2]^2 # x^3-tet[1]*(tet[1]^2+3*tet[2]^2)
  47.  
  48.   f <- cbind(m1,m2,m3)
  49.  
  50.   return(f)
  51.  
  52. }
  53.  
  54. g4 <- function(tet,x) {
  55. # normal distribution moments: http://en.wikipedia.org/wiki/Normal_distribution#Moments
  56.  
  57.   m1 <- (tet[1]-x)
  58.   m2 <- (tet[2]^2 - (x - tet[1])^2)
  59.   m3 <- x^3-(tet[1]^3+3*tet[1]*tet[2]^2)
  60.   m4 <- x^4-(tet[1]^4+6*tet[1]^2*tet[2]^2+3*tet[2]^4)
  61.  
  62.   f <- cbind(m1,m2,m3,m4)
  63.  
  64.   return(f)
  65.  
  66. }
  67.  
  68. ###
  69.  
  70. x.mean; mean(x)
  71. x.sd; sd(x)
  72.  
  73. o <- optim(c(0,1),objFun) #method="SANN",control=list(maxit=50000))
  74.  
  75. o$par
  76.  
  77. t0 <- c(mean(x),sd(x))
  78.  
  79. x.gmm <- gmm(g2,x,t0) # 2 momenty
  80. coef(x.gmm)
  81.  
  82. x.gmm <- gmm(g3,x,t0) # 3 momenty
  83. coef(x.gmm)
  84.  
  85. x.gmm <- gmm(g4,x,t0) # 3 momenty
  86. coef(x.gmm)
RAW Paste Data