#################################################################### # # STAT 718A Practical 6 # ##################################################################### Q1. We shall investigate the use of Iterative Conditional Modes (ICM) and MCMC for segmenting a noisy image into two classes. Download the program mcmc-practical6.r from the course web page and store in your local space. Run the program a few times to get a feel for what it does. You will see that there is a true image, which then has Gaussian noise added with standard deviation 1. The ICM estimate is calculate using an Ising model prior with parameter beta = 1.5. Then a Metropolis Hastings algorithm is run which simulates from the posterior distribution. The final simulated value can also be thought of as an estimate of the classes in the image. Measures of goodness of fit are printed out for each estimate, where the error is the sum of square differences from the estimate and the truth. ################# Now investigate the effect of the Ising prior. The prior here is the eight neighbor prior with parameter beta. Edit the file and change beta = 1.5 to other values, e.g. 0, 0.5, 1, 2. Which beta parameter is best here? What is the effect of increasing or decreasing beta? ################## Now change the true image to the disk pattern, by commenting and uncommenting the `true' image part of the code. Investigate which prior beta value is best here. ################## Now change the true image to the stripes pattern, by commenting and uncommenting the `true' image part of the code. Investigate which prior beta value is best here. ################## Q2. Shape analysis using maximum likelihood estimation. #type library(shapes) #then copy and paste the function in isomle.txt into R. #This will update the isomle function in the shapes package. type ans1 <- isomle( qcet2.dat ) ans2 <- isomle( qset2.dat ) plotshapes(ans1\$mshape,ans2\$mshape) #Likelihood ratio test k<-dim(qcet2.dat)[1] pool.dat <- array( 0, c(k,2,53) ) pool.dat[,,1:30] <- qcet2.dat pool.dat[,,31:53] <- qset2.dat ans3 <- isomle( pool.dat ) plotshapes( ans1\$mshape) points( ans2\$mshape, col=2) points( ans3\$mshape, col=3) # look at variability print( c(ans1\$tau , ans2\$tau , ans3\$tau ) ) #H0: means are equal and tau parameters are equal L0 <- ans3\$loglike L1 <- ans1\$loglike + ans2\$loglike #-2 (loglikelihood ratio) statistics test <- -2*(L0 - L1) print(test) #Compute p-value print( 1 - pchisq( test, df=(2*k-3) ) ) Carry out a similar analysis of the Schizophrenia data where the groups are: A <- schizophrenia.dat[,,1:14] #controls B <- schizophrenia.dat[,,15:28] #patients Is there a significant different between the patients and controls ############################################# Now try some alternative tests....we shall describe these in detail in the class next week (on Wednesday). ############################################## # Goodall F test (assumes isotropy) testmeanshapes(qset2.dat,qcet2.dat,Hotelling=FALSE) # Hotelling T^2 test testmeanshapes(qset2.dat,qcet2.dat,Hotelling=TRUE) # Permutation test resampletest(qset2.dat,qcet2.dat,resamples=100,permutation=TRUE) ################################################ # Goodall F test (assumes isotropy) testmeanshapes(A,B,Hotelling=FALSE) # Hotelling T^2 test testmeanshapes(A,B,Hotelling=TRUE) # Permutation test resampletest(A,B,resamples=100,permutation=TRUE)