Plot birth pulse data: last capture in December 2015

There are 32 births overall; 5 in 2013, 15 in 2014, 12 in 2015. First birth in July 2013; last in April 2015 Observations started July 2013 (t=0.54), ended June 2015 (t=2.46). Set t0= 1 Jan 2013; birth times approximated as midway through the month.

Plot hazard of birth pulse: mu = 0.443

dark line: cv = 1; light gray lines: cv = 1.1-1.5 (dashed) or 0.6-0.8 (dotted) Birth pulse narrows as the cv becomes larger

#hazard function
get_hazard_birth
## function (cv, time) 
## {
##     tfrac <- time%%1
##     mu <- 0.443
##     if (cv < 1/sqrt(3)) {
##         alpha <- 1 + cv * sqrt(3)
##         beta <- (2 * cv * sqrt(3))/(1 + cv * sqrt(3))
##     }
##     if (cv >= 1/sqrt(3)) {
##         alpha <- 1.5 * (1 + (cv^2))
##         beta <- 0.75 * (1 + (cv^2))
##     }
##     if (tfrac < 0.5 & tfrac >= 0) {
##         val <- mu * alpha * max((1 - 2 * beta * tfrac), 0)
##     }
##     else if (tfrac >= 0.5 & tfrac <= 1) {
##         val <- mu * alpha * max((1 - 2 * beta * (1 - tfrac)), 
##             0)
##     }
##     else {
##         val <- "error in evaluating maximum, time value is wierd"
##     }
##     return(val)
## }

Plot probability of birth over time, functions to get likelihood

dark lines represent cv=1; Calculated as exp(get_log_birth_like(time, t0 = 0, sigma=1))

# calculate log likelihood, set initial time = August 2014 (t=1.5)
get_log_birth_like
## function (t, t0, cv) 
## {
##     mu <- 0.443
##     if (cv < 1/sqrt(3)) {
##         talpha <- 1 + cv * sqrt(3)
##         tbeta <- (2 * cv * sqrt(3))/(1 + cv * sqrt(3))
##     }
##     if (cv >= 1/sqrt(3)) {
##         talpha <- 1.5 * (1 + (cv^2))
##         tbeta <- 0.75 * (1 + (cv^2))
##     }
##     tfrac <- t%%1
##     d <- t0 + t
##     df <- d%%1
##     c <- t0
##     cf <- c%%1
##     if (tfrac < 0.5 & tfrac >= 0) {
##         hazard <- mu * talpha * max((1 - 2 * tbeta * tfrac), 
##             0)
##     }
##     else if (tfrac >= 0.5 & tfrac <= 1) {
##         hazard <- mu * talpha * max((1 - 2 * tbeta * (1 - tfrac)), 
##             0)
##     }
##     H1 <- (1 - (tbeta/2)) * (floor(d) - floor(c) - 1)
##     ifelse(cf < 0.5, H2 <- 0.5 * (1 - (tbeta/2)) + (0.5 - cf) * 
##         (1 - tbeta + (tbeta * (0.5 - cf))), H2 <- (1 - cf) * 
##         (1 - tbeta + (tbeta * (1 - cf))))
##     ifelse(df < 0.5, H3 <- df * (1 - (tbeta * df)), H3 <- 0.5 * 
##         (1 - tbeta/2) + (df - 0.5) * (1 - tbeta + tbeta * (df - 
##         0.5)))
##     H4 <- (1/(2 * tbeta)) * (floor(d) - floor(c) - 1)
##     if (cf < (1/(2 * tbeta))) {
##         H5 <- (1/(4 * tbeta)) + tbeta * (((1/(2 * tbeta)) - cf)^2)
##     }
##     else if (cf >= (1/(2 * tbeta)) & cf < 1 - (1/(2 * tbeta))) {
##         H5 <- 1/(4 * tbeta)
##     }
##     else {
##         H5 <- (1 - cf) * (1 - (tbeta * (1 - cf)))
##     }
##     if (df < (1/(2 * tbeta))) {
##         H6 <- df * (1 - tbeta * df)
##     }
##     else if (df >= (1/(2 * tbeta)) & df < 1 - (1/(2 * tbeta))) {
##         H6 <- 1/(4 * tbeta)
##     }
##     else {
##         H6 <- (1/(4 * tbeta)) + tbeta * (df - (1 - (1/(2 * tbeta))))^2
##     }
##     ifelse(cv < 1/sqrt(3), H <- mu * talpha * (H1 + H2 + H3), 
##         H <- mu * talpha * (H4 + H5 + H6))
##     Surv <- exp(-H)
##     val <- hazard * Surv
##     LL <- log(val)
##     return(LL)
## }
# minimize -loglikelihood; parameters cv
birthNLL = function(x, cv){ 
  val <- 0
  for (i in 1:length(x)){
    temp <- get_log_birth_like(t = x[i], t0 = 0.04166, cv = cv)
    val <- temp + val
  }
  return (-val) 
}  

# minimize -loglikelihood; parameteres cv and t0
birthNLL_t0cv = function(x, params){ 
  # params[1]= t0
  # params[2]= cv
  val <- 0
  for (i in 1:length(x)){
    temp <- get_log_birth_like(t = x[i], t0 = params[1], cv = params[2])
    val <- temp + val
  }
  return (-val)
}

Test: plot Likelihood function and estimate CV assuming t0=0

Set t0 = 0.0417 = mid-January previous year

Estimate CV and t0