# x <- matrix( c(
# rnorm(50, mean=.3, sd=.5),
# rnorm(50, mean=.16, sd=.1),
# rnorm(50, mean=.4, sd=.3),
# rnorm(50, mean=.6, sd=.2),
# rnorm(50, mean=.2, sd=.2)
# ), ncol=2)

membership <- function(x, classes)
	{
	
	#setup display
	par(mfrow=c(3,classes))
	
	#compute classes
	#k-means, use 10 starts so results converge
	cl <- kmeans(x, classes, nstart=10)
	
	#c-means
	cc <- cmeans(x,classes)
	
	#pam with standardization
	c.pam <- pam(x, classes, stand=TRUE)
	
	#display k-means classes
	for(i in 1:classes)
		{
		#compute ellipsoid hulls for c-means clusters
		exy <- ellipsoidhull( matrix(c(x[,1][which(cl$cluster == i)], x[,2][which(cl$cluster == i)]), ncol=2 ) )
	
		plot(x, col = cl$cluster, main="K-means: Hard Classes", cex=0.5, pch=16, xlab="", ylab="", cex.main=.8)
		points(cl$centers, col = 1:classes, pch = 8, cex=2)
		lines(predict(exy))
		}
		
# 	#display hard c-means classes
# 	for(i in 1:classes)
# 		{
# 		#compute ellipsoid hulls for c-means clusters
# 		exy <- ellipsoidhull( matrix(c(x[,1][which(cc$cluster == i)], x[,2][which(cc$cluster == i)]), ncol=2 ) )
# 	
# 		plot(x, col = cc$cluster, main="C-means: Hard Classes", cex=0.5, pch=16, xlab="", ylab="", cex.main=.8)
# 		points(cc$centers, col = 1:classes, pch = 8, cex=2)
# 		lines(predict(exy))
# 		}
		
	#display pam classes
	for(i in 1:classes)
		{
		#compute ellipsoid hulls for c-means clusters
		exy <- ellipsoidhull( matrix(c(x[,1][which(c.pam$clustering == i)], x[,2][which(c.pam$clustering == i)]), ncol=2 ) )
	
		plot(x, col = c.pam$clustering, main="PAM with Standardization", cex=0.5, pch=16, xlab="", ylab="", cex.main=.8)
		points(c.pam$medoids, col = 1:classes, pch = 8, cex=2)
		lines(predict(exy))
		}
	
	#display c-means fuzzy membership
	for(i in 1:classes)
		{
		#compute ellipsoid hulls for c-means clusters
		exy <- ellipsoidhull( matrix(c(x[,1][which(cc$cluster == i)], x[,2][which(cc$cluster == i)]), ncol=2 ) )
		
		plot(x, main=paste("C-means: Fuzzy Class Membership ", i ), type="n", xlab="", ylab="", cex.main=.8)
		points(x, col=gray(1-cc$membership[,i]), cex=0.5, pch=16)
# 		points(x, col=rgb(cc$membership[,i], 0, 0), cex=0.5, pch=16)
		lines(predict(exy))
		}
	
	}