### Code for calculating posterior probabilities of K for Structure runs
# 5/18/09
# Jennie Lavine and Mike Antolin


setwd('C:/Documents and Settings/slcv/Desktop')

#read in data - this is ln(P(D)) for K=10:13, 4 runs each
dat10<--c(2569.9, 2576.9,2576.3, 2572.0)
dat11<--c(2558.5,2575.2,2593.9,2554.3)
dat12<--c(2539.0,2533.7,2548.8 ,2548.2)
dat13<--c(2580.2,2583.9,2562.8,2569.1)

data<-as.data.frame(cbind('10'=dat10, '11'=dat11, '12'=dat12, '13'=dat13) )
 ## define some matrices
 stand <- array(0, c(4,4))
 expstand<-array(0, c(4,4))
 post_prob <- array(0, c(4,4))

# Use each set of runs as a replicate
for (i in 1:4)  {
stand[,i]<-data[,i]-min(data[,i])     # Standardize the by subtracting off the lowest value
expstand[,i]<-exp(stand[,i])      #exponentiate 
post_prob[,i]<-expstand[,i]/sum(expstand[,i])   # divide by sum of all model likelihoods to calculate posterior probabilities
}
k<-c(10,10,10,10,11,11,11,11,12,12,12,12,13,13,13,13)
P<-as.vector(post_prob)
#plot with a smooth line
plot(c(10,10,10,10,11,11,11,11,12,12,12,12,13,13,13,13),as.vector(post_prob))

plot(post_prob,ylab="Posterior P(K)",xlab='Number of Clusters (K)', pch=16,col=4)   #plot
lines(lowess(post_prob),col=4)