Mercurial > repos > greg > ks_distribution

comparison ks_distribution.R @ 4:a91bd45aa8b1 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author greg
date Wed, 31 May 2017 07:55:32 -0400
parents 5ace8af6edb6
children 22cae2172406
comparison
equal deleted inserted replaced
3:e293a5736ae9 4:a91bd45aa8b1
15 15
16 get_num_components = function(components_data) 16 get_num_components = function(components_data)
17 { 17 {
18 # Get the max of the number_comp column. 18 # Get the max of the number_comp column.
19 number_comp = components_data[, 3] 19 number_comp = components_data[, 3]
20 num_components <- max(number_comp, na.rm=TRUE); 20 num_components <- max(number_comp, na.rm=TRUE)
21 num_components 21 return(num_components)
22 } 22 }
23 23
24 get_pi_mu_var = function(components_data, num_components) 24 get_pi_mu_var = function(components_data, num_components)
25 { 25 {
26 # FixMe: enhance this to generically handle any integer value for num_components. 26 # FixMe: enhance this to generically handle any integer value for num_components.
27 if (num_components == 1) 27 if (num_components == 1)
28 { 28 {
29 pi <- c(components_data[1, 9]); 29 pi <- c(components_data[1, 9])
30 mu <- c(components_data[1, 7]); 30 mu <- c(components_data[1, 7])
31 var <- c(components_data[1, 8]); 31 var <- c(components_data[1, 8])
32 } 32 }
33 else if (num_components == 2) 33 else if (num_components == 2)
34 { 34 {
35 pi <- c(components_data[2, 9], components_data[3, 9]); 35 pi <- c(components_data[2, 9], components_data[3, 9])
36 mu <- c(components_data[2, 7], components_data[3, 7]); 36 mu <- c(components_data[2, 7], components_data[3, 7])
37 var <- c(components_data[2, 8], components_data[3, 8]); 37 var <- c(components_data[2, 8], components_data[3, 8])
38 } 38 }
39 else if (num_components == 3) 39 else if (num_components == 3)
40 { 40 {
41 pi <- c(components_data[4, 9], components_data[5, 9], components_data[6, 9]); 41 pi <- c(components_data[4, 9], components_data[5, 9], components_data[6, 9])
42 mu <- c(components_data[4, 7], components_data[5, 7], components_data[6, 7]); 42 mu <- c(components_data[4, 7], components_data[5, 7], components_data[6, 7])
43 var <- c(components_data[4, 8], components_data[5, 8], components_data[6, 8]); 43 var <- c(components_data[4, 8], components_data[5, 8], components_data[6, 8])
44 } 44 }
45 else if (num_components == 4) 45 else if (num_components == 4)
46 { 46 {
47 pi <- c(components_data[7, 9], components_data[8, 9], components_data[9, 9], components_data[10, 9]); 47 pi <- c(components_data[7, 9], components_data[8, 9], components_data[9, 9], components_data[10, 9])
48 mu <- c(components_data[7, 7], components_data[8, 7], components_data[9, 7], components_data[10, 7]); 48 mu <- c(components_data[7, 7], components_data[8, 7], components_data[9, 7], components_data[10, 7])
49 var <- c(components_data[7, 8], components_data[8, 8], components_data[9, 8], components_data[10, 8]); 49 var <- c(components_data[7, 8], components_data[8, 8], components_data[9, 8], components_data[10, 8])
50 } 50 }
51 return = c(pi, mu, var) 51 else if (num_components == 5)
52 return 52 {
53 pi <- c(components_data[11, 9], components_data[12, 9], components_data[13, 9], components_data[14, 9], components_data[15, 9])
54 mu <- c(components_data[11, 7], components_data[12, 7], components_data[13, 7], components_data[14, 7], components_data[15, 7])
55 var <- c(components_data[11, 8], components_data[12, 8], components_data[13, 8], components_data[14, 8], components_data[15, 8])
56 }
57 else if (num_components == 6)
58 {
59 pi <- c(components_data[16, 9], components_data[17, 9], components_data[18, 9], components_data[19, 9], components_data[20, 9], components_data[21, 9])
60 mu <- c(components_data[16, 7], components_data[17, 7], components_data[18, 7], components_data[19, 7], components_data[20, 7], components_data[21, 7])
61 var <- c(components_data[16, 8], components_data[17, 8], components_data[18, 8], components_data[19, 8], components_data[20, 8], components_data[21, 8])
62 }
63 results = c(pi, mu, var)
64 return(results)
53 } 65 }
54 66
55 plot_ks<-function(kaks_input, output, pi, mu, var) 67 plot_ks<-function(kaks_input, output, pi, mu, var, max_ks)
56 { 68 {
57 # Start PDF device driver to save charts to output. 69 # Start PDF device driver to save charts to output.
58 pdf(file=output, bg="white") 70 pdf(file=output, bg="white")
71 kaks <- read.table(file=kaks_input, header=T)
72 max_ks <- max(kaks$Ks, na.rm=TRUE)
59 # Change bin width 73 # Change bin width
60 bin <- 0.05 * seq(0, 40); 74 max_bin_range <- as.integer(max_ks / 0.05)
61 kaks <- read.table(file=kaks_input, header=T); 75 bin <- 0.05 * seq(0, max_bin_range)
62 kaks <- kaks[kaks$Ks<2,]; 76 kaks <- kaks[kaks$Ks<max_ks,];
63 h.kst <- hist(kaks$Ks, breaks=bin, plot=F); 77 h.kst <- hist(kaks$Ks, breaks=bin, plot=F)
64 nc <- h.kst$counts; 78 nc <- h.kst$counts
65 vx <- h.kst$mids; 79 vx <- h.kst$mids
66 ntot <- sum(nc); 80 ntot <- sum(nc)
67 # Set margin for plot bottom, left top, right. 81 # Set margin for plot bottom, left top, right.
68 par(mai=c(0.5, 0.5, 0, 0)); 82 par(mai=c(0.5, 0.5, 0, 0))
69 # Plot dimension in inches. 83 # Plot dimension in inches.
70 par(pin=c(2.5, 2.5)); 84 par(pin=c(2.5, 2.5))
71 g <- calculate_fitted_density(pi, mu, var); 85 g <- calculate_fitted_density(pi, mu, var)
72 h <- ntot * 2.5 / sum(g); 86 h <- ntot * 2.5 / sum(g)
73 vx <- seq(1, 100) * 0.02; 87 vx <- seq(1, 100) * 0.02
74 ymax <- max(nc) + 5; 88 ymax <- max(nc) + 5
75 barplot(nc, space=0.25, offset=0, width=0.04, xlim=c(0,2), ylim=c(0, ymax)); 89 barplot(nc, space=0.25, offset=0, width=0.04, xlim=c(0, max_ks), ylim=c(0, ymax), col="lightpink1", border="lightpink3")
76 # Add x-axis. 90 # Add x-axis.
77 axis(1); 91 axis(1)
78 color <- c('green', 'blue', 'black', 'red'); 92 color <- c('red', 'yellow','green','black','blue', 'darkorange' )
79 for (i in 1:length(mu)) 93 for (i in 1:length(mu))
80 { 94 {
81 lines(vx, g[,i] * h, lwd=2, col=color[i]); 95 lines(vx, g[,i] * h, lwd=2, col=color[i])
82 } 96 }
83 }; 97 }
84 98
85 calculate_fitted_density <- function(pi, mu, var) 99 calculate_fitted_density <- function(pi, mu, var)
86 { 100 {
87 comp <- length(pi); 101 comp <- length(pi)
88 var <- var/mu^2; 102 var <- var/mu^2
89 mu <- log(mu); 103 mu <- log(mu)
90 #calculate lognormal density 104 # Calculate lognormal density.
91 vx <- seq(1, 100) * 0.02; 105 vx <- seq(1, 100) * 0.02
92 fx <- matrix(0, 100, comp); 106 fx <- matrix(0, 100, comp)
93 for (i in 1:100) 107 for (i in 1:100)
94 { 108 {
95 for (j in 1:comp) 109 for (j in 1:comp)
96 { 110 {
97 fx[i, j] <- pi[j] * dlnorm(vx[i], meanlog=mu[j], sdlog=(sqrt(var[j]))); 111 fx[i, j] <- pi[j] * dlnorm(vx[i], meanlog=mu[j], sdlog=(sqrt(var[j])))
98 }; 112 if (is.nan(fx[i,j])) fx[i,j]<-0
99 }; 113 }
100 fx; 114 }
115 return(fx)
101 } 116 }
102 117
103 # Read in the components data. 118 # Read in the components data.
104 components_data <- read.delim(opt$components_input, header=TRUE); 119 components_data <- read.delim(opt$components_input, header=TRUE)
105 # Get the number of components. 120 # Get the number of components.
106 num_components <- get_num_components(components_data) 121 num_components <- get_num_components(components_data)
107 122
108 # Set pi, mu, var. 123 # Set pi, mu, var.
109 items <- get_pi_mu_var(components_data, num_components); 124 items <- get_pi_mu_var(components_data, num_components)
110 pi <- items[1]; 125 pi <- items[1:3]
111 mu <- items[2]; 126 mu <- items[4:6]
112 var <- items[3]; 127 var <- items[7:9]
113 128
114 # Plot the output. 129 # Plot the output.
115 plot_ks(opt$kaks_input, opt$output, pi, mu, var); 130 plot_ks(opt$kaks_input, opt$output, pi, mu, var, max_ks)