############################################################# # STAT 718A Computer Practical 3 ############################################################# # Part 1 library(shapes) # First of all we'll read in some 3D data from a file # and carry out some analysis #download and save the file "monkey.csv" from the course webpage # read it in tem <- read.csv("C:/directoryname/monkey.csv") #type the name `tem' to look at it tem # The dataset consists of n=19 individuals, with k=7 landmarks # in m=3 dimensions from one side of the skull. We need to get # the data into an array of size 7 x 3 x 19 # convert from dataframe to a matrix tem<-as.matrix(tem) dim(tem) #Note the first column has a landmark label (number) in it, which we don't need. tem1 <- tem[,2:58] # convert to an array monk <- array(tem1, c(7, 3, 19) ) # just check it looks OK... e.g. right dimension, looks like the file coordinates # are in the right order? monk # Let's plot the landmarks shapes3d( monk ) # you can rotated the plot with the mouse to visualize the 3D nature of the plot # here's a 2D projection of the data plotshapes( monk ) # and another projection plotshapes( monk , orthproj=c(1,3) ) # carry out Procrustes analysis proc <- procGPA( monk ) # plot Procrustes rotated data shapes3d( proc\$rotated ) # Note the 3d commands use the R package RGL which in general adds # to the current plot shapes3d( proc\$rotated , col=4, joinline=c(1:7) ) # To clear the current plot you can just quit the plotting window. # Investigate the principal components of shape # using the commands # shapepca( ) # proc\$percent gives the percentage of variability # Also, investigate centroid size and riemannian distance # proc\$size , proc\$rho # What are the main features of the analysis? Are there any unusual # observations? If so can you see how they differ? # Carry out the analysis without any unusual observations. ########################################################################### # Part 2 # Let's investigate how the Procrustes mean performs for some # simulated data n<-10 sigma<-0.1 mu<- matrix( c(-0.5,0,0.5,0,0,sqrt(3)/2) , 2, 3) mu <- t(mu) x <- array( 0, c(3,2,n) ) plotshapes (mu) for (i in 1:n){ x[,,i] <- mu + matrix( sigma*rnorm(6) , 3 ,2) } plotshapes(x) # Bookstein mean bb <- bookstein2d( x ) plotshapes(x) ans <- procGPA( x ) #convert Procrustes mean to Bookstein SV for plotting pp <- bookstein.shpv( ans\$mshape ) lines( pp[c(1,2,3,1),] , col=2) # True mean lines( mu[c(1,2,3,1),] , col=1) # Add the Bookstein mean lines( bb\$mshape[c(1,2,3,1),] , col=3) ############################################# # Let's investigate mean square error from 100 Monte Carlo # simulations n<-30 sigma<-0.1 x <- array( 0, c(3,2,n) ) biasbook<-mu*0 biasprocf<-mu*0 biasprocp<-mu*0 msebook<-0 mseprocf<-0 mseprocp<-0 for (isim in 1:100){ for (i in 1:n){ x[,,i] <- mu + matrix( sigma*rnorm(6) , 3 ,2) } bb <- bookstein2d( x ) ansf <- procGPA( x ) # full Procrustes ansp <- procGPA( x , scale=FALSE) # partial Procrustes biasbook <- biasbook + bb\$mshape - mu biasprocf <- biasprocf + bookstein.shpv(ansf\$mshape) - mu biasprocp <- biasprocp + bookstein.shpv(ansp\$mshape) - mu msebook <- msebook + norm(bb\$mshape - mu)**2 mseprocf <- mseprocf + norm(bookstein.shpv(ansf\$mshape) - mu)**2 mseprocp <- mseprocp + norm(bookstein.shpv(ansp\$mshape) - mu)**2 } biasbook <- biasbook/100 biasprocf <- biasprocf/100 biasprocp <- biasprocp/100 msebook <- msebook/100 mseprocf <- mseprocf/100 mseprocp <- mseprocp/100 # Norm of the Bias print( c(norm(biasbook), norm(biasprocf), norm(biasprocp) ) ) # MSE print( c( msebook, mseprocf, mseprocp) ) # Investigate the performance of the estimators # for different n and different sigma. Which is best ?