#############################################################
#              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 ?