Mercurial > repos > greg > ks_distribution

comparison ks_distribution.R @ 44:69ac2edffdde draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author greg
date Mon, 26 Jun 2017 13:58:37 -0400
parents 5074cb26597a
children f8b303b9f455
comparison
equal deleted inserted replaced
43:5074cb26597a 44:69ac2edffdde
5 option_list <- list( 5 option_list <- list(
6 make_option(c("-c", "--components_input"), action="store", dest="components_input", help="Ks significant components input dataset"), 6 make_option(c("-c", "--components_input"), action="store", dest="components_input", help="Ks significant components input dataset"),
7 make_option(c("-k", "--kaks_input"), action="store", dest="kaks_input", help="KaKs analysis input dataset"), 7 make_option(c("-k", "--kaks_input"), action="store", dest="kaks_input", help="KaKs analysis input dataset"),
8 make_option(c("-n", "--number_comp"), action="store", dest="number_comp", type="integer", help="Number of significant components in the Ks distribution"), 8 make_option(c("-n", "--number_comp"), action="store", dest="number_comp", type="integer", help="Number of significant components in the Ks distribution"),
9 make_option(c("-o", "--output"), action="store", dest="output", help="Output dataset"), 9 make_option(c("-o", "--output"), action="store", dest="output", help="Output dataset"),
10 make_option(c("-r", "--colors"), action="store", default="", help="List of component colors") 10 make_option(c("-r", "--colors"), action="store", default=NA, help="List of component colors")
11 ) 11 )
12 12
13 parser <- OptionParser(usage="%prog [options] file", option_list=option_list) 13 parser <- OptionParser(usage="%prog [options] file", option_list=option_list)
14 args <- parse_args(parser, positional_arguments=TRUE) 14 args <- parse_args(parser, positional_arguments=TRUE)
15 opt <- args$options 15 opt <- args$options
16 16
17 set_component_colors = function(colors, number_comp) 17 set_component_colors = function(colors, number_comp)
18 { 18 {
19 # Handle colors for components. 19 # Handle colors for components.
20 if (nchar(colors) == 0) 20 if (is.na(colors) == 0) {
21 {
22 # Randomly specify colors for components. 21 # Randomly specify colors for components.
23 specified_colors <- c('red', 'yellow', 'green', 'black', 'blue', 'darkorange') 22 specified_colors <- c('red', 'yellow', 'green', 'black', 'blue', 'darkorange')
24 } 23 } else {
25 else
26 {
27 # Handle selected colors for components. 24 # Handle selected colors for components.
28 specified_colors <- c() 25 specified_colors <- c()
29 items <- strsplit(colors, ",") 26 items <- strsplit(colors, ",")
30 for (item in items) 27 for (item in items) {
31 {
32 specified_colors <- c(specified_colors, item) 28 specified_colors <- c(specified_colors, item)
33 } 29 }
34 num_colors_specified <- length(specified_colors) 30 num_colors_specified <- length(specified_colors)
35 if (num_colors_specified < number_comp) 31 if (num_colors_specified < number_comp) {
36 { 32 for (i in num_colors_specified:number_comp) {
37 for (i in num_colors_specified:number_comp) 33 if (!any(specified_colors == 'red')) {
38 {
39 if (!any(specified_colors == 'red'))
40 {
41 specified_colors <- c(specified_colors, 'red') 34 specified_colors <- c(specified_colors, 'red')
42 } 35 } else if (!any(specified_colors == 'yellow')) {
43 else if (!any(specified_colors == 'yellow'))
44 {
45 specified_colors <- c(specified_colors, 'yellow') 36 specified_colors <- c(specified_colors, 'yellow')
46 } 37 } else if (!any(specified_colors == 'green')) {
47 else if (!any(specified_colors == 'green'))
48 {
49 specified_colors <- c(specified_colors, 'green') 38 specified_colors <- c(specified_colors, 'green')
50 } 39 } else if (!any(specified_colors == 'black')) {
51 else if (!any(specified_colors == 'black'))
52 {
53 specified_colors <- c(specified_colors, 'black') 40 specified_colors <- c(specified_colors, 'black')
54 } 41 } else if (!any(specified_colors == 'blue')) {
55 else if (!any(specified_colors == 'blue'))
56 {
57 specified_colors <- c(specified_colors, 'blue') 42 specified_colors <- c(specified_colors, 'blue')
58 } 43 } else {
59 else
60 {
61 specified_colors <- c(specified_colors, 'darkorange') 44 specified_colors <- c(specified_colors, 'darkorange')
62 } 45 }
63 } 46 }
64 } 47 }
65 } 48 }
67 } 50 }
68 51
69 get_pi_mu_var = function(components_data, number_comp) 52 get_pi_mu_var = function(components_data, number_comp)
70 { 53 {
71 # FixMe: enhance this to generically handle any integer value for number_comp. 54 # FixMe: enhance this to generically handle any integer value for number_comp.
72 if (number_comp == 1) 55 if (number_comp == 1) {
73 {
74 pi <- c(components_data[1, 9]) 56 pi <- c(components_data[1, 9])
75 mu <- c(components_data[1, 7]) 57 mu <- c(components_data[1, 7])
76 var <- c(components_data[1, 8]) 58 var <- c(components_data[1, 8])
77 } 59 } else if (number_comp == 2) {
78 else if (number_comp == 2)
79 {
80 pi <- c(components_data[2, 9], components_data[3, 9]) 60 pi <- c(components_data[2, 9], components_data[3, 9])
81 mu <- c(components_data[2, 7], components_data[3, 7]) 61 mu <- c(components_data[2, 7], components_data[3, 7])
82 var <- c(components_data[2, 8], components_data[3, 8]) 62 var <- c(components_data[2, 8], components_data[3, 8])
83 } 63 } else if (number_comp == 3) {
84 else if (number_comp == 3)
85 {
86 pi <- c(components_data[4, 9], components_data[5, 9], components_data[6, 9]) 64 pi <- c(components_data[4, 9], components_data[5, 9], components_data[6, 9])
87 mu <- c(components_data[4, 7], components_data[5, 7], components_data[6, 7]) 65 mu <- c(components_data[4, 7], components_data[5, 7], components_data[6, 7])
88 var <- c(components_data[4, 8], components_data[5, 8], components_data[6, 8]) 66 var <- c(components_data[4, 8], components_data[5, 8], components_data[6, 8])
89 } 67 } else if (number_comp == 4) {
90 else if (number_comp == 4)
91 {
92 pi <- c(components_data[7, 9], components_data[8, 9], components_data[9, 9], components_data[10, 9]) 68 pi <- c(components_data[7, 9], components_data[8, 9], components_data[9, 9], components_data[10, 9])
93 mu <- c(components_data[7, 7], components_data[8, 7], components_data[9, 7], components_data[10, 7]) 69 mu <- c(components_data[7, 7], components_data[8, 7], components_data[9, 7], components_data[10, 7])
94 var <- c(components_data[7, 8], components_data[8, 8], components_data[9, 8], components_data[10, 8]) 70 var <- c(components_data[7, 8], components_data[8, 8], components_data[9, 8], components_data[10, 8])
95 } 71 } else if (number_comp == 5) {
96 else if (number_comp == 5)
97 {
98 pi <- c(components_data[11, 9], components_data[12, 9], components_data[13, 9], components_data[14, 9], components_data[15, 9]) 72 pi <- c(components_data[11, 9], components_data[12, 9], components_data[13, 9], components_data[14, 9], components_data[15, 9])
99 mu <- c(components_data[11, 7], components_data[12, 7], components_data[13, 7], components_data[14, 7], components_data[15, 7]) 73 mu <- c(components_data[11, 7], components_data[12, 7], components_data[13, 7], components_data[14, 7], components_data[15, 7])
100 var <- c(components_data[11, 8], components_data[12, 8], components_data[13, 8], components_data[14, 8], components_data[15, 8]) 74 var <- c(components_data[11, 8], components_data[12, 8], components_data[13, 8], components_data[14, 8], components_data[15, 8])
101 } 75 } else if (number_comp == 6) {
102 else if (number_comp == 6)
103 {
104 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]) 76 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])
105 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]) 77 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])
106 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]) 78 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])
107 } 79 }
108 results = c(pi, mu, var) 80 results = c(pi, mu, var)
132 vx <- seq(1, 100) * (max_ks / 100) 104 vx <- seq(1, 100) * (max_ks / 100)
133 ymax <- max(nc) 105 ymax <- max(nc)
134 barplot(nc, space=0.25, offset=0, width=0.04, xlim=c(0, max_ks), ylim=c(0, ymax), col="lightpink1", border="lightpink3") 106 barplot(nc, space=0.25, offset=0, width=0.04, xlim=c(0, max_ks), ylim=c(0, ymax), col="lightpink1", border="lightpink3")
135 # Add x-axis. 107 # Add x-axis.
136 axis(1) 108 axis(1)
137 for (i in 1:length(mu)) 109 for (i in 1:length(mu)) {
138 {
139 lines(vx, g[,i] * h, lwd=2, col=specified_colors[i]) 110 lines(vx, g[,i] * h, lwd=2, col=specified_colors[i])
140 } 111 }
141 } 112 }
142 113
143 calculate_fitted_density <- function(pi, mu, var, max_ks) 114 calculate_fitted_density <- function(pi, mu, var, max_ks)
146 var <- var/mu^2 117 var <- var/mu^2
147 mu <- log(mu) 118 mu <- log(mu)
148 # Calculate lognormal density. 119 # Calculate lognormal density.
149 vx <- seq(1, 100) * (max_ks / 100) 120 vx <- seq(1, 100) * (max_ks / 100)
150 fx <- matrix(0, 100, comp) 121 fx <- matrix(0, 100, comp)
151 for (i in 1:100) 122 for (i in 1:100) {
152 { 123 for (j in 1:comp) {
153 for (j in 1:comp)
154 {
155 fx[i, j] <- pi[j] * dlnorm(vx[i], meanlog=mu[j], sdlog=(sqrt(var[j]))) 124 fx[i, j] <- pi[j] * dlnorm(vx[i], meanlog=mu[j], sdlog=(sqrt(var[j])))
156 if (is.nan(fx[i,j])) 125 if (is.nan(fx[i,j])) {
157 {
158 fx[i,j]<-0 126 fx[i,j]<-0
159 } 127 }
160 } 128 }
161 } 129 }
162 return(fx) 130 return(fx)
169 # Set component colors. 137 # Set component colors.
170 specified_colors <- set_component_colors(opt$colors, number_comp) 138 specified_colors <- set_component_colors(opt$colors, number_comp)
171 139
172 # Set pi, mu, var. 140 # Set pi, mu, var.
173 items <- get_pi_mu_var(components_data, number_comp) 141 items <- get_pi_mu_var(components_data, number_comp)
174 if (number_comp == 1) 142 if (number_comp == 1) {
175 {
176 pi <- items[1] 143 pi <- items[1]
177 mu <- items[2] 144 mu <- items[2]
178 var <- items[3] 145 var <- items[3]
179 } 146 } else if (number_comp == 2) {
180 else if (number_comp == 2)
181 {
182 pi <- items[1:2] 147 pi <- items[1:2]
183 mu <- items[3:4] 148 mu <- items[3:4]
184 var <- items[5:6] 149 var <- items[5:6]
185 } 150 } else if (number_comp == 3) {
186 else if (number_comp == 3)
187 {
188 pi <- items[1:3] 151 pi <- items[1:3]
189 mu <- items[4:6] 152 mu <- items[4:6]
190 var <- items[7:9] 153 var <- items[7:9]
191 } 154 } else if (number_comp == 4) {
192 else if (number_comp == 4)
193 {
194 pi <- items[1:4] 155 pi <- items[1:4]
195 mu <- items[5:8] 156 mu <- items[5:8]
196 var <- items[9:12] 157 var <- items[9:12]
197 } 158 } else if (number_comp == 5) {
198 else if (number_comp == 5)
199 {
200 pi <- items[1:5] 159 pi <- items[1:5]
201 mu <- items[6:10] 160 mu <- items[6:10]
202 var <- items[11:15] 161 var <- items[11:15]
203 } 162 } else if (number_comp == 6) {
204 else if (number_comp == 6)
205 {
206 pi <- items[1:6] 163 pi <- items[1:6]
207 mu <- items[7:12] 164 mu <- items[7:12]
208 var <- items[13:18] 165 var <- items[13:18]
209 } 166 }
210 167