Adding packages

R/aedseo.R

@@ -23,7 +23,7 @@
 #' fraction of observables in the window of size k that are allowed to be NA.
 #' @param season_weeks A numeric vector of length 2, `c(start,end)`, with the
 #' start and end weeks of the seasons to stratify the observations by. Must
-#' spand the new year; ex: `season = c(21,20)`. Default, `NULL`, is no
+#' spand the new year; ex: `season_weeks = c(21,20)`. Default, `NULL`, is no
 #' stratification by season.
 #'
 #' @return A `aedseo` object containing:

R/fit_peak.R

@@ -95,7 +95,7 @@ fit_peak <- function(
     init_params <- switch(family,
       weibull = log(c(1.5, mean(observations))),
       lnorm = c(mean(log(observations)),
-                sd(log(observations))),
+                stats::sd(log(observations))),
       exp = log(1.5)
     )
     return(init_params)
@@ -104,27 +104,30 @@ fit_peak <- function(
   # The weighted negative loglikelihood function
   nll <- function(par, weighted_observations, family = "weibull") {
     switch(family,
-      weibull = -sum(dweibull(weighted_observations$observation,
-                              shape = exp(par[1]), scale = exp(par[2]),
-                              log = TRUE) * weighted_observations$weight),
-      lnorm = -sum(dlnorm(weighted_observations$observation,
-                          meanlog =  par[1], sdlog = par[2],
-                          log = TRUE) * weighted_observations$weight),
-      exp =  -sum(dexp(weighted_observations$observation, rate = exp(par[1]),
-                       log = TRUE) * weighted_observations$weight)
+      weibull = -sum(stats::dweibull(weighted_observations$observation,
+                                     shape = exp(par[1]), scale = exp(par[2]),
+                                     log = TRUE) *
+                       weighted_observations$weight),
+      lnorm = -sum(stats::dlnorm(weighted_observations$observation,
+                                 meanlog =  par[1], sdlog = par[2],
+                                 log = TRUE) * weighted_observations$weight),
+      exp =  -sum(stats::dexp(weighted_observations$observation,
+                              rate = exp(par[1]),
+                              log = TRUE) * weighted_observations$weight)
     )
   }
 
   # Run optimisation for weighted observations
-  optim_obj <- optim(par = init_par_fun(family = family,
-                                        observations =
-                                          weighted_observations$observation),
-                     fn = nll,
-                     weighted_observations = weighted_observations,
-                     family = family,
-                     method = optim_method,
-                     lower = lower_optim,
-                     upper = upper_optim)
+  optim_obj <-
+    stats::optim(par = init_par_fun(family = family,
+                                    observations =
+                                      weighted_observations$observation),
+                 fn = nll,
+                 weighted_observations = weighted_observations,
+                 family = family,
+                 method = optim_method,
+                 lower = lower_optim,
+                 upper = upper_optim)
 
   # Back-transform optimized parameters to their original scale if needed.
   # This is done by exponentiating the parameters, as they were

tests/testthat/test-fit_peak.R

@@ -0,0 +1,23 @@
+test_that("Can run the fit_peak fun without error", {
+
+  # Weibull fit
+  obs <- 10
+  season <- c("2018/2019", "2019/2020", "2020/2021")
+  season_num_rev <- rev(seq(from = 1, to = length(season)))
+  observations <- rep(rnorm(10, obs), length(season))
+
+  peak_input <- tibble(
+    observation = observations,
+    weight = 0.8^rep(season_num_rev, each = obs),
+    season = season[rep(seq(from = 1, to = length(season)), each = obs)]
+  )
+
+  expect_no_error(fit_peak(weighted_observations = peak_input))
+
+  # Exp fit
+  expect_no_error(fit_peak(weighted_observations = peak_input,
+                           family = "exp",
+                           optim_method = "Brent",
+                           lower_optim = 0,
+                           upper_optim = 1000))
+})