# GMM normal distribution with R

Aug 28th, 2012
373
0
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)