[R] Re: Gaussian Mixtures Models

[Noel Yvonnick]

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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._