Mercurial > repos > greg > ks_distribution

comparison ks_distribution.R @ 7:22cae2172406 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author greg
date Tue, 06 Jun 2017 09:02:08 -0400
parents a91bd45aa8b1
children 214e2710c51e
comparison
equal deleted inserted replaced
6:9831319a19fb 7:22cae2172406
62 } 62 }
63 results = c(pi, mu, var) 63 results = c(pi, mu, var)
64 return(results) 64 return(results)
65 } 65 }
66 66
67 plot_ks<-function(kaks_input, output, pi, mu, var, max_ks) 67 plot_ks<-function(kaks_input, output, pi, mu, var)
68 { 68 {
69 # Start PDF device driver to save charts to output. 69 # Start PDF device driver to save charts to output.
70 pdf(file=output, bg="white") 70 pdf(file=output, bg="white")
71 kaks <- read.table(file=kaks_input, header=T) 71 kaks <- read.table(file=kaks_input, header=T)
72 max_ks <- max(kaks$Ks, na.rm=TRUE) 72 max_ks <- max(kaks$Ks, na.rm=TRUE)
73 # Change bin width 73 # Change bin width
74 max_bin_range <- as.integer(max_ks / 0.05) 74 max_bin_range <- as.integer(max_ks / 0.05)
75 bin <- 0.05 * seq(0, max_bin_range) 75 bin <- 0.05 * seq(0, (max_bin_range + 1 ))
76 kaks <- kaks[kaks$Ks<max_ks,]; 76 kaks <- kaks[kaks$Ks<max_ks,]
77 h.kst <- hist(kaks$Ks, breaks=bin, plot=F) 77 h.kst <- hist(kaks$Ks, breaks=bin, plot=F)
78 nc <- h.kst$counts 78 nc <- h.kst$counts
79 vx <- h.kst$mids 79 vx <- h.kst$mids
80 ntot <- sum(nc) 80 ntot <- sum(nc)
81 # Set margin for plot bottom, left top, right. 81 # Set margin for plot bottom, left top, right.
82 par(mai=c(0.5, 0.5, 0, 0)) 82 par(mai=c(0.5, 0.5, 0, 0))
83 # Plot dimension in inches. 83 # Plot dimension in inches.
84 par(pin=c(2.5, 2.5)) 84 par(pin=c(3.0, 3.0))
85 g <- calculate_fitted_density(pi, mu, var) 85 g <- calculate_fitted_density(pi, mu, var, max_ks)
86 h <- ntot * 2.5 / sum(g) 86 h <- ntot * 1.5 / sum(g)
87 vx <- seq(1, 100) * 0.02 87 vx <- seq(1, 100) * (max_ks / 100)
88 ymax <- max(nc) + 5 88 ymax <- max(nc)
89 barplot(nc, space=0.25, offset=0, width=0.04, xlim=c(0, max_ks), ylim=c(0, ymax), col="lightpink1", border="lightpink3") 89 barplot(nc, space=0.25, offset=0, width=0.04, xlim=c(0, max_ks), ylim=c(0, ymax), col="lightpink1", border="lightpink3")
90 # Add x-axis. 90 # Add x-axis.
91 axis(1) 91 axis(1)
92 color <- c('red', 'yellow','green','black','blue', 'darkorange' ) 92 color <- c('red', 'yellow','green','black','blue', 'darkorange' )
93 for (i in 1:length(mu)) 93 for (i in 1:length(mu))
94 { 94 {
95 lines(vx, g[,i] * h, lwd=2, col=color[i]) 95 lines(vx, g[,i] * h, lwd=2, col=color[i])
96 } 96 }
97 } 97 }
98 98
99 calculate_fitted_density <- function(pi, mu, var) 99 calculate_fitted_density <- function(pi, mu, var, max_ks)
100 { 100 {
101 comp <- length(pi) 101 comp <- length(pi)
102 var <- var/mu^2 102 var <- var/mu^2
103 mu <- log(mu) 103 mu <- log(mu)
104 # Calculate lognormal density. 104 # Calculate lognormal density.
105 vx <- seq(1, 100) * 0.02 105 vx <- seq(1, 100) * (max_ks / 100)
106 fx <- matrix(0, 100, comp) 106 fx <- matrix(0, 100, comp)
107 for (i in 1:100) 107 for (i in 1:100)
108 { 108 {
109 for (j in 1:comp) 109 for (j in 1:comp)
110 { 110 {
120 # Get the number of components. 120 # Get the number of components.
121 num_components <- get_num_components(components_data) 121 num_components <- get_num_components(components_data)
122 122
123 # Set pi, mu, var. 123 # Set pi, mu, var.
124 items <- get_pi_mu_var(components_data, num_components) 124 items <- get_pi_mu_var(components_data, num_components)
125 pi <- items[1:3] 125 if (num_components == 1)
126 mu <- items[4:6] 126 {
127 var <- items[7:9] 127 pi <- items[1]
128 mu <- items[2]
129 var <- items[3]
130 }
131 if (num_components == 2)
132 {
133 pi <- items[1:2]
134 mu <- items[3:4]
135 var <- items[5:6]
136 }
137 if (num_components == 3)
138 {
139 pi <- items[1:3]
140 mu <- items[4:6]
141 var <- items[7:9]
142 }
143 if (num_components == 4)
144 {
145 pi <- items[1:4]
146 mu <- items[5:8]
147 var <- items[9:12]
148 }
149 if (num_components == 5)
150 {
151 pi <- items[1:5]
152 mu <- items[6:10]
153 var <- items[11:15]
154 }
155 if (num_components == 6)
156 {
157 pi <- items[1:6]
158 mu <- items[7:12]
159 var <- items[13:18]
160 }
128 161
129 # Plot the output. 162 # Plot the output.
130 plot_ks(opt$kaks_input, opt$output, pi, mu, var, max_ks) 163 plot_ks(opt$kaks_input, opt$output, pi, mu, var)