In jrminter/particlesizeR: Convenience functions for particle size analysis

library(knitr)
options(width = 72)

First, let's load the required packages.

library(particlesizeR)
library(mixtools)
library(ggplot2)

Now, let's generate a two Gaussian mixture. Because we consider this a particle size distribution, we need to require the values to be positive.

samples <- gen_two_gaussian_mixture(10000, 0.5, 3.0, 1.0, 12.0, 2.0)
samples <- samples[samples > 0]

Let's now construct a histogram object. This will give us access to arrays of midpoints and density values.

library(ggplot2)

df <- data.frame(samples=samples)
plt <- ggplot(df, aes(x=samples)) + geom_histogram(bins=250)
print(plt)

Let's look at what is in a ggplot2 histogram, starting at the beginning:

pg <- ggplot_build(plt)
head(pg$data[[1]])

Note that we have access to a vector of midpoints (x) as well as a vector of lower boundaries (xmin) and a vector of upper boundaries (xmax) as well as the density and the ncount.

Next, look at the end of the data...

tail(pg$data[[1]])

Eventually, we want to plot the results. We will build a plot using the newer ggplot2 package. We will also generate a series of x (diameter) values to model.

First, we will generate the data

truth <- gen_two_gaussian_dist(0.5, 3.0, 1.0, 12.0, 2.0, 0, 20, 0.1)

verbose = FALSE

if (verbose == TRUE){
  print(c(min(truth), max(truth)))
  print(max(truth$dist))
  print(max(pg$data[[1]]$density))
  print(head(truth))
}

pt_line_size <- 1
df_work <- data.frame(diam=pg$data[[1]]$x, density=pg$data[[1]]$density)
plt <- ggplot(df_work, aes(diam,density)) + 
       geom_point(size=pt_line_size, colour='darkblue' ) +
       geom_line(data=truth, aes(x,dist), colour='darkred', size=pt_line_size) +
       ylab("density") +
       xlab("diameter [nm]") +
       ggtitle("Diameter distribution of two Gaussian distributions") +
       theme(axis.text=element_text(size=12), axis.title=element_text(size=14),
              plot.title = element_text(hjust = 0.5)) + # center the title
      NULL

print(plt)

We want to compute the mixture model. We will make a guess at the parameters from the histogram.

out <- normalmixEM(samples, lambda=0.5, mu=c( 2, 9.), sigma=c(.8, 2.))

We can now compute the fitted values:

x <- seq(from=0, to=20, by=0.1)
fit <- out$lambda[1]*dnorm(x, out$mu[1], out$sigma[1]) + 
       out$lambda[2]*dnorm(x, out$mu[2],out$sigma[2])

df_fit <- data.frame(ecd=x, fit=fit)

And make a final plot.

``````r pt_line_size <- 1 final_plt <- ggplot() + geom_point(data=df_work, aes(diam,density), colour="darkblue", size=pt_line_size) + geom_line(data=df_fit, aes(ecd, fit), colour="darkred", size=pt_line_size) + xlab("size [nm]") + ylab("density") + scale_x_continuous(breaks = seq(from = 0, to = 20, by = 1), limits = c(0, 20)) + scale_y_continuous(breaks = seq(from = 0, to = 0.25, by = 0.05), limits = c(0, 0.25)) + ggtitle("Diameter distribution from two Gaussian distributions") + labs(caption = '[e-mail] / @jrminter') + theme(axis.text=element_text(size=12), axis.title=element_text(size=12), plot.title=element_text(hjust = 0.5)) + # center the title NULL print(final_plt)

```

jrminter/particlesizeR documentation built on June 5, 2019, 11:05 p.m.