I join below a piece of code I have written for GMM in the special case of
spherical noise. Note that the end condition is a fixed number of iterations.
you may want to replace it by some criterion on the change in loglikelihood
or something like that.
You may want to test it by typing:
# generate some bivariate data
t <- seq(.15,3.05,length=100)
X <- cbind(t,t+1.25*sin(2*t)) + matrix(rnorm(200,0,.1),100,2)
# GMM with 20 classes (or "points")
gm <- gmm(X,npts=20)
You may also want to free the variance parameters (var.equal=FALSE) and the
prior probabilities (prior.probs=TRUE), though it is not recommended to free
both of them.
Hope this helps,
--
Yvonnick Noel
Universite de Lille 3
UFR de psychologie
#####################################################################################
gmm <- function(X,npts,tmax=40,var.equal=TRUE,prior.prob=FALSE,plot.true=T)
{
# Needed packages
require(mva)
require(scatterplot3d)
#require(MASS)
X <- as.matrix(X)
n <- nrow(X)
p <- ncol(X)
priors <- matrix(1,n,1) %*% matrix(1/npts,1,npts)
# For drawing distribution width in the 2D case
theta <- seq(0,2*pi,length=360)
circle <- cbind(cos(theta),sin(theta))
# Starting model: sample of observed data
model <- X[sample(npts),]
# Initialize inverse variance parameter
d2 <- dist2(X,model)
beta <- rep((n*p*npts)/sum(d2),npts)
for(t in 1:tmax)
{
# Probabilities
R <- exp(-.5 * (d2 %*% diag(beta))) * priors
# Responsibilities
R <- R / apply(R,1,"sum")
G <- apply(R,2,"sum")
model <- (t(R)/G) %*% X
# Plotting (each 5 iterations)
if((plot.true)&&(((t%/%5)-(t/5))==0))
{
if(p==2)
{
plot(X,pch=19,cex=.5,col="red",xlab="",ylab="",main="")
points(model,col="blue",pch=19,cex=.7)
for(l in 1:npts)
{
icirc <- circle * 1.96 * sqrt(1/beta[l])
lines((matrix(1,360,1) %*% model[l,]) +
icirc,col="darkgreen")
}
}
else if(p==3)
{
gph <-
scatterplot3d(X[,c(1,3,2)],highlight.3d=TRUE,pch=19,xlab="",ylab="",zlab="",main="")
gph$points3d(model[,c(1,3,2)],col="blue",pch=19)
}
else if(p>3)
{
pc<-princomp(model)
gph <-
scatterplot3d(predict(pc,X)[,c(1,3,2)],highlight.3d=TRUE,pch=19,xlab="",ylab="",zlab="",main="")
gph$points3d(pc$scores[,c(1,3,2)],col="blue",pch=19)
}
}
d2 <- dist2(X,model)
# Constrain to equal variances
if(var.equal) { beta <- rep((n*p)/sum(R*d2),npts) }
# In case of unequal variances
else { beta <- (p*G)/apply(R*d2,2,"sum") }
# Take prior probabilities into account...
if(prior.prob) { priors <- matrix(1,n,1) %*% (G / n) }
}
list(model=model,beta=beta,d2=d2,priors=priors)
}
############################################################################################
dist2 <- function(X,Y)
{
if(missing(Y)) Y <- X
if(!is.matrix(X)) X <- as.matrix(X)
if(!is.matrix(Y)) Y <- as.matrix(Y)
dimx<-dim(X)
dimy<-dim(Y)
sqrx<-diag(X%*%t(X))
sqry<-diag(Y%*%t(Y))
(sqrx %*% matrix(1,1,dimy[1])) + (matrix(1,dimx[1],1) %*% sqry) -
(2*X%*%t(Y))
}
###########################################################################################
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at
stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._