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


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) )

means = apply(data, 2,mean) #mean of each K value
stand<-means-min(means)     # Standardize the by subtracting off the lowest value
expstand<-exp(stand)      #exponentiate 
post_prob<-expstand/sum(estand)   # divide by sum of all model likelihoods to calculate posterior probabilities

#plot with a smooth line
plot(10:13,post_prob,ylab="Posterior P(K)",xlab='Number of Clusters (K)', pch=16,col=4)   #plot
lines(lowess(10:13,post_prob),col=4)