diff ks_distribution.R @ 7:22cae2172406 draft

Uploaded
author greg
date Tue, 06 Jun 2017 09:02:08 -0400
parents a91bd45aa8b1
children 214e2710c51e
line wrap: on
line diff
--- a/ks_distribution.R	Wed May 31 12:14:32 2017 -0400
+++ b/ks_distribution.R	Tue Jun 06 09:02:08 2017 -0400
@@ -64,7 +64,7 @@
     return(results)
 }
 
-plot_ks<-function(kaks_input, output, pi, mu, var, max_ks)
+plot_ks<-function(kaks_input, output, pi, mu, var)
 {
     # Start PDF device driver to save charts to output.
     pdf(file=output, bg="white")
@@ -72,8 +72,8 @@
     max_ks <- max(kaks$Ks, na.rm=TRUE)
     # Change bin width
     max_bin_range <- as.integer(max_ks / 0.05)
-    bin <- 0.05 * seq(0, max_bin_range)
-    kaks <- kaks[kaks$Ks<max_ks,];
+    bin <- 0.05 * seq(0, (max_bin_range + 1 ))
+    kaks <- kaks[kaks$Ks<max_ks,]
     h.kst <- hist(kaks$Ks, breaks=bin, plot=F)
     nc <- h.kst$counts
     vx <- h.kst$mids
@@ -81,11 +81,11 @@
     # Set margin for plot bottom, left top, right.
     par(mai=c(0.5, 0.5, 0, 0))
     # Plot dimension in inches.
-    par(pin=c(2.5, 2.5))
-    g <- calculate_fitted_density(pi, mu, var)
-    h <- ntot * 2.5 / sum(g)
-    vx <- seq(1, 100) * 0.02
-    ymax <- max(nc) + 5
+    par(pin=c(3.0, 3.0))
+    g <- calculate_fitted_density(pi, mu, var, max_ks)
+    h <- ntot * 1.5 / sum(g)
+    vx <- seq(1, 100) * (max_ks / 100)
+    ymax <- max(nc)
     barplot(nc, space=0.25, offset=0, width=0.04, xlim=c(0, max_ks), ylim=c(0, ymax), col="lightpink1", border="lightpink3")
     # Add x-axis.
     axis(1)
@@ -96,13 +96,13 @@
     }
 }
 
-calculate_fitted_density <- function(pi, mu, var)
+calculate_fitted_density <- function(pi, mu, var, max_ks)
 {
     comp <- length(pi)
     var <- var/mu^2
     mu <- log(mu)
     # Calculate lognormal density.
-    vx <- seq(1, 100) * 0.02
+    vx <- seq(1, 100) * (max_ks / 100)
     fx <- matrix(0, 100, comp)
     for (i in 1:100)
     {
@@ -122,9 +122,42 @@
 
 # Set pi, mu, var.
 items <- get_pi_mu_var(components_data, num_components)
-pi <- items[1:3]
-mu <- items[4:6]
-var <- items[7:9]
+if (num_components == 1)
+{
+	pi <- items[1]
+	mu <- items[2]
+	var <- items[3]
+}
+if (num_components == 2)
+{
+	pi <- items[1:2]
+	mu <- items[3:4]
+	var <- items[5:6]
+}
+if (num_components == 3)
+{
+	pi <- items[1:3]
+	mu <- items[4:6]
+	var <- items[7:9]
+}
+if (num_components == 4)
+{
+	pi <- items[1:4]
+	mu <- items[5:8]
+	var <- items[9:12]
+}
+if (num_components == 5)
+{
+	pi <- items[1:5]
+	mu <- items[6:10]
+	var <- items[11:15]
+}
+if (num_components == 6)
+{
+	pi <- items[1:6]
+	mu <- items[7:12]
+	var <- items[13:18]
+}
 
 # Plot the output.
-plot_ks(opt$kaks_input, opt$output, pi, mu, var, max_ks)
+plot_ks(opt$kaks_input, opt$output, pi, mu, var)