# @author: joao lopes # @workplace: Instituto Gulbenkian de Ciencia # @date: 12th May 2011 # Function to perform the regression step on a file with rejection step results. # It also estimate the mode and the 0.95 credible interval of the obtained posterior distributions. # Creates .txt file with results in the directory from where it is run # @arg data_file - file with summary of 'real' data (.trg) # @arg rej_file - file with rejection step results (.rej) get_modes <- function(data_file , rej_file){ #import locfit library library(locfit) #import Mark Beaumont's scripts source("loc2plot_d.r") source("make_pd2.r") #demographic parameters' priors (minimum and maximum values) mintev1 <- 0 maxtev1 <- 50000 mintev2 <- 0 maxtev2 <- 1200000 minNe1 <- 0 maxNe1 <- 5000 minNe2 <- 0 maxNe2 <- 5000 minNe3 <- 0 maxNe3 <- 5000 minNeA1 <- 0 maxNeA1 <- 10000 minNeA2 <- 0 maxNeA2 <- 10000 minmig1 <- 0 maxmig1 <- 0.05 minmig2 <- 0 maxmig2 <- 0.05 minmig3 <- 0 maxmig3 <- 0.05 minmigA1 <- 0 maxmigA1 <- 0.05 #credible interval to use credible <- 0.05 prob <- 0.05 #import the .rej file abc.rej <- data.matrix(read.table(rej_file)) #import the .trg files target <- data.matrix(read.table(data_file)) #created a .txt file to save results write("mode 0.025 0.975 \n", "estimates.txt", append=F) #MUTATION RATE #save mutation rate to estimate.txt write("mutation rate estimate", "estimates.txt", append=T) write(loc1statsx(x=abc.rej[,3], prob=prob), "estimates.txt", append=T) write(" ", "estimates.txt", append=T) #SPLITTING TIMES #save Tev1 to estimate.txt write("Tev1 estimate", "estimates.txt", append=T) write(loc1statsx(x=abc.rej[,7], prob=prob), "estimates.txt", append=T) #save Tev2 to estimate.txt write("Tev2 estimate", "estimates.txt", append=T) write(loc1statsx(x=abc.rej[,8], prob=prob), "estimates.txt", append=T) #save Tev3 to estimate.txt #Ne ESTIMATES #save Ne1 to estimate.txt write("Ne1 estimate", "estimates.txt", append=T) write(loc1statsx(x=abc.rej[,9], prob=prob), "estimates.txt", append=T) #save Ne2 to estimate.txt write("Ne2 estimate", "estimates.txt", append=T) write(loc1statsx(x=abc.rej[,10], prob=prob), "estimates.txt", append=T) #save Ne3 to estimate.txt write("Ne3 estimate", "estimates.txt", append=T) write(loc1statsx(x=abc.rej[,11], prob=prob), "estimates.txt", append=T) #ANCIENT Ne ESTIMATES #save NeA1 to estimate.txt write("NeA1 estimate", "estimates.txt", append=T) write(loc1statsx(x=abc.rej[,12], prob=prob), "estimates.txt", append=T) #save NeA2 to estimate.txt write("NeA2 estimate", "estimates.txt", append=T) write(loc1statsx(x=abc.rej[,13], prob=prob), "estimates.txt", append=T) #MIGRATION ESTIMATES #save Mig1 to estimate.txt write("Mig1 estimate", "estimates.txt", append=T) write(loc1statsx(x=abc.rej[,14], prob=prob), "estimates.txt", append=T) #save Mig2 to estimate.txt write("Mig2 estimate", "estimates.txt", append=T) write(loc1statsx(x=abc.rej[,15], prob=prob), "estimates.txt", append=T) #save Mig3 to estimate.txt write("Mig3 estimate", "estimates.txt", append=T) write(loc1statsx(x=abc.rej[,16], prob=prob), "estimates.txt", append=T) #ANCIENT MIGRATION ESTIMATES #save MigA1 to estimate.txt write("MigA1 estimate", "estimates.txt", append=T) write(loc1statsx(x=abc.rej[,17], prob=prob), "estimates.txt", append=T) write(" ", "estimates.txt", append=T) write(" ", "estimates.txt", append=T) write(" ", "estimates.txt", append=T) write(" ", "estimates.txt", append=T) write(" ", "estimates.txt", append=T) write("if you like my R scripting, I have a used car in my driveway to sell you", "estimates.txt", append=T) }