Raw Values

library(ggplot2)
library(reshape2)

# Load ECan data (stored locally)
data <- read.csv('RangioraWinter2016.csv', stringsAsFactors=FALSE)
names(data) <- c('date', 'time', 'wind.speed', 'wind.dir', 'wind.dir.std',
                 'wind.speed.std', 'wind.max', 'co', 'temp.ground',
                 'temp.2m', 'temp.6m', 'pm10', 'pm2.5', 'pm.course')
data$date <- as.POSIXct(data$date, format = '%m/%d/%Y %H:%M', tz='Etc/GMT-12')
data$datel <- as.POSIXlt(data$date)
data$time <- NULL # this is included in date

# Stipulate start and end times of initial and final colocations
data.1.start <- as.POSIXct('2016-07-11 00:00 NZST')
data.1.end <- as.POSIXct('2016-07-26 00:00 NZST')
data.2.start <- as.POSIXct('2016-10-03 00:00 NZST')
data.2.end <- as.POSIXct('2016-10-19 00:00 NZST')

# Split into initial and final colocation datasets
data.1 <- data[which(data$date<data.1.end & data$date>data.1.start), ]
data.2 <- data[which(data$date<data.2.end & data$date>data.2.start), ]

# Remove NAs
data.1 <- data.1[which(!is.na(data.1$pm2.5)), ]
data.2 <- data.2[which(!is.na(data.2$pm2.5)), ]

# Add 'dataset' field
data.1[,'dataset'] <- rep('data.1', nrow(data.1))
data.2[,'dataset'] <- rep('data.2', nrow(data.2))

PlotDensity <- function() {
  # Plot the densities of data.1 and data.2
  
  # Get densities
  dens.1 <- density(x=data.1$pm2.5)
  dens.2 <- density(x=data.2$pm2.5)
  dens.1 <- data.frame(pm2.5=dens.1$x, density=dens.1$y,
                       dataset=rep("data.1", length(dens.1$x)))
  dens.2 <- data.frame(pm2.5=dens.2$x, density=dens.2$y,
                       dataset=rep("data.2", length(dens.2$x)))
  densities <- rbind(dens.1, dens.2)
  
  # Plot
  plt <- ggplot() + 
    geom_line(data=densities, aes(pm2.5, density, colour=dataset))
  
  return(plt)
}

PlotHist <- function() {
  # Plot a histogram of data.1 and data.2
  
  data <- rbind(data.1, data.2)
  
  plt <- ggplot() +
    geom_histogram(aes(x=pm2.5, y=..density.., fill=dataset), binwidth=1, 
                   data=data.1, alpha=0.5) +
    geom_histogram(aes(x=pm2.5, y=..density.., fill=dataset), binwidth=1, 
                   data=data.2, alpha=0.5)
  
  return (plt)
}

PlotHist()

Positivise all values

Weibull distributions only have non-zero density for positive values. We therefore need to truncate all negative values at some small positive value.

# Truncate PM concentrations at zero
data.1[which(data.1$pm2.5<=0), 'pm2.5'] <- 1e-10
data.2[which(data.2$pm2.5<=0), 'pm2.5'] <- 1e-10

# Plot
PlotHist()

Initial eyeball of Weibull parameters

AddWeibull <- function(plt, params.list) {
  # Add Weibull curves to an existing plot
  #
  # Args:
  #   plt: An existing ggplot2 plot
  #   params.list: A list of lists of Weibull distribution parameters of the
  #                form, each including a $shape, $scale, and $dataset field
  #
  # Returns:
  #   Another ggplot2 plot
  
  # The range of x values used
  x <- x <- seq(0, 100, length.out=100)
  
  # Create a dataframe for the density plot
  wei.dens <- data.frame(pm2.5=NA, density=NA, dataset=NA)
  for (params in params.list) {
    wei.dens.2 <- data.frame(pm2.5=x, 
                            density=dweibull(x, shape=params$shape, 
                                             scale=params$scale),
                            dataset=rep(params$dataset, length(x)))
    wei.dens <- rbind(wei.dens, wei.dens.2)
  }
  
  # Strip initial NAs
  wei.dens <- wei.dens[2:nrow(wei.dens), ]
  
  plt <- plt + 
    geom_line(data=wei.dens, aes(pm2.5, density, colour=dataset), size=1)
  
  return(plt)
}

# Eyeballed parameters
wei.iball.1 <- list(shape=1.1, scale=15, dataset='data.1')
wei.iball.2 <- list(shape=2, scale=6, dataset='data.2')
wei.iball <- list(wei.iball.1, wei.iball.2)

# Plot
plt <- PlotHist()
AddWeibull(plt, wei.iball)

Fitted Weibull distribution

library(MASS)

# Fit Weibull distributions to data
wei.params.1 <- as.list(fitdistr(data.1$pm2.5, densfun=dweibull, 
                                 start=list(shape=1.1, scale=15)))$estimate
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
wei.params.1$dataset <- 'data.1'
## Warning in wei.params.1$dataset <- "data.1": Coercing LHS to a list
wei.params.2 <- as.list(fitdistr(data.2$pm2.5, densfun=dweibull, 
                                 start=list(shape=2, scale=6)))$estimate
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
wei.params.2$dataset <- 'data.2'
## Warning in wei.params.2$dataset <- "data.2": Coercing LHS to a list
wei.params <- list(wei.params.1, wei.params.2)

# Print fitted Weibull parameters
wei.params.1
## $shape
## [1] 0.2899302
## 
## $scale
## [1] 6.466588
## 
## $dataset
## [1] "data.1"
wei.params.2
## $shape
## [1] 0.5902768
## 
## $scale
## [1] 4.379898
## 
## $dataset
## [1] "data.2"
# The size of the two datasets
nrow(data.1)
## [1] 2080
nrow(data.2)
## [1] 2272
# Plot
plt <- PlotHist()
AddWeibull(plt, wei.params)

```