• API
• FAQ
• Tools
• Archive
SHARE
TWEET

GMM normal distribution with R

mjaniec Aug 28th, 2012 169 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
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.

Top