Add Weibull

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: monty
 Title: Monte Carlo Models
-Version: 0.3.13
+Version: 0.3.14
 Authors@R: c(person("Rich", "FitzJohn", role = c("aut", "cre"),
                     email = "[email protected]"),
              person("Wes", "Hinsley", role = "aut"),

NAMESPACE

@@ -60,13 +60,15 @@ export(monty_random_n_poisson)
 export(monty_random_n_real)
 export(monty_random_n_truncated_normal)
 export(monty_random_n_uniform)
+export(monty_random_n_weibull)
 export(monty_random_negative_binomial_mu)
 export(monty_random_negative_binomial_prob)
 export(monty_random_normal)
 export(monty_random_poisson)
 export(monty_random_real)
 export(monty_random_truncated_normal)
 export(monty_random_uniform)
+export(monty_random_weibull)
 export(monty_rng)
 export(monty_rng_create)
 export(monty_rng_distributed_pointer)

R/dsl-distributions.R

@@ -246,6 +246,17 @@ distr_uniform <- distribution(
   },
   cpp = list(density = "uniform", sample = "uniform"))
 
+distr_weibull <- distribution(
+  name = "Weibull",
+  density = function(x, shape, scale) dweibull(x, shape, scale, log = TRUE),
+  domain = c(0, Inf),
+  sample = "monty_random_weibull",
+  expr = list(
+    density = quote(log(shape) + (shape - 1) * log(x) - shape * log(scale) - 
+                      (x / scale)^shape),
+    mean = quote(scale * gamma(1 + 1 / k))),
+  cpp = list(density = "weibull", sample = "weibull"))
+
 dsl_distributions <- local({
   d <- list(
     distr_beta,
@@ -263,7 +274,8 @@ dsl_distributions <- local({
     distr_normal,
     distr_poisson,
     distr_truncated_normal,
-    distr_uniform)
+    distr_uniform,
+    distr_weibull)
   split(d, vapply(d, "[[", "", "name"))
 })
 

R/random.R

@@ -613,3 +613,27 @@ monty_random_n_truncated_normal <- function(n_samples, mean, sd, min, max,
                                             state) {
   cpp_monty_random_n_truncated_normal(n_samples, mean, sd, min, max, state)
 }
+
+
+##' Sample from a Weibull distribution
+##'
+##' @title Sample from Weibull
+##'
+##' @param shape Shape
+##'
+##' @param scale Scale
+##'
+##' @inheritParams monty_random_real
+##' @inherit monty_random_real return
+##'
+##' @export
+monty_random_weibull <- function(shape, scale, state) {
+  cpp_monty_random_weibull(shape, scale, state)
+}
+
+
+##' @export
+##' @rdname monty_random_weibull
+monty_random_n_weibull <- function(n_samples, shape, scale, state) {
+  cpp_monty_random_n_weibull(n_samples, shape, scale, state)
+}

_pkgdown.yml

@@ -115,6 +115,8 @@ reference:
     - monty_random_n_hypergeometric
     - monty_random_truncated_normal
     - monty_random_n_truncated_normal
+    - monty_random_weibull
+    - monty_random_n_weibull
   - subsection: Legacy interface
     desc: >-
       Functions imported from dust; these will be retired soon

inst/include/monty/random/density.hpp

@@ -282,5 +282,16 @@ __host__ __device__ T cauchy(T x, T location, T scale, bool log) {
   return maybe_log(ret, log);
 }
 
+template <typename T>
+__host__ __device__ T weibull(T x, T shape, T scale, bool log) {
+#ifndef __CUDA_ARCH__
+  static_assert(std::is_floating_point<T>::value,
+                "weibull should only be used with real types");
+#endif
+  const auto ret = monty::math::log(shape) + (shape - 1) * monty::math::log(x) -
+    shape * monty::math::log(scale) - (x / scale)^shape;
+  return maybe_log(ret, log);
+}
+
 }
 }

inst/include/monty/random/random.hpp

@@ -16,6 +16,7 @@
 #include "monty/random/normal.hpp"
 #include "monty/random/poisson.hpp"
 #include "monty/random/uniform.hpp"
+#include "monty/random/weibull.hpp"
 
 #include "monty/random/version.hpp"
 

inst/include/monty/random/weibull.hpp

@@ -0,0 +1,41 @@
+#pragma once
+
+#include <cmath>
+
+#include "monty/random/generator.hpp"
+#include "monty/random/numeric.hpp"
+#include "monty/random/math.hpp"
+
+namespace monty {
+namespace random {
+
+namespace {
+
+template <typename real_type>
+void weibull_validate(real_type shape, real_type scale) {
+  if (shape < 0.0 || scale < 0.0) {
+    char buffer[256];
+    snprintf(buffer, 256,
+             "Invalid call to Weibull with shape = %g, scale = %g",
+             shape, scale);
+    monty::utils::fatal_error(buffer);
+  }
+}
+}
+
+template <typename real_type, typename rng_state_type>
+real_type weibull(rng_state_type& rng_state, real_type shape, real_type scale) {
+#ifdef __CUDA_ARCH__
+  static_assert("weibull() not implemented for GPU targets");
+#endif
+  weibull_validate(shape, scale);
+
+  if (rng_state.deterministic) {
+    return scale * std::tgamma(1 + 1 / shape);
+  }
+  const real_type u = random_real<real_type>(rng_state);
+  return scale * monty::math::pow(-monty::math::log(1 - u), 1 / shape);
+}
+
+}
+}

src/random.cpp

@@ -619,6 +619,26 @@ cpp11::doubles cpp_monty_random_n_uniform(size_t n_samples,
                                  min, max, "min", "max");
 }
 
+// weibull
+[[cpp11::register]]
+cpp11::doubles cpp_monty_random_weibull(cpp11::doubles shape,
+                                        cpp11::doubles scale,
+                                        cpp11::sexp ptr) {
+  const auto fn = [](auto& state, auto shape, auto scale) { return monty::random::weibull<double>(state, shape, scale); };
+  return monty_random_sample_1_2(fn, ptr, "weibull",
+                                 shape, scale, "shape", "scale");
+}
+
+[[cpp11::register]]
+cpp11::doubles cpp_monty_random_n_weibull(size_t n_samples,
+                                          cpp11::doubles shape,
+                                          cpp11::doubles scale,
+                                          cpp11::sexp ptr) {
+  const auto fn = [](auto& state, auto shape, auto scale) { return monty::random::weibull<double>(state, shape, scale); };
+  return monty_random_sample_n_2(fn, n_samples, ptr, "weibull",
+                                 shape, scale, "shape", "scale");
+}
+
 //// 3-arg functions
 
 // beta_binomial_prob

tests/testthat/test-rng.R

@@ -1360,3 +1360,47 @@ test_that("can generate from truncated normal from lower tail", {
   })
   expect_gt(sum(res > 0.05), 5)
 })
+
+
+test_that("can draw weibull random numbers", {
+  shape <- 5
+  scale <- 3
+  n <- 10000000
+
+  ans1 <- monty_random_n_weibull(n, shape, scale, monty_rng_create(seed = 1))
+  ans2 <- monty_random_n_weibull(n, shape, scale, monty_rng_create(seed = 1))
+  expect_identical(ans1, ans2)
+
+  expect_equal(mean(ans1), scale * gamma(1 + 1 / shape), tolerance = 1e-3)
+  true_var <- scale^2 * (gamma(1 + 2 / shape) - gamma(1 + 1 / shape)^2)
+  expect_equal(var(ans1), true_var, tolerance = 1e-3)
+})
+
+
+test_that("deterministic weibull returns mean", {
+  n_reps <- 10
+  shape <- as.numeric(sample(10, n_reps, replace = TRUE))
+  scale <- as.numeric(sample(10, n_reps, replace = TRUE))
+
+  rng <- monty_rng_create(seed = 1, deterministic = TRUE)
+  state <- monty_rng_state(rng)
+
+  expect_equal(
+    mapply(monty_random_weibull, shape, scale, MoreArgs = list(rng)),
+    scale * gamma(1 + 1 / shape))
+  expect_equal(monty_rng_state(rng), state)
+})
+
+
+test_that("weibull random numbers prevent bad inputs", {
+  r <- monty_rng_create(seed = 1)
+  expect_equal(monty_random_weibull(0, 0, r), 0)
+  expect_equal(monty_random_weibull(Inf, Inf, r), Inf)
+
+  expect_error(
+    monty_random_weibull(-1.1, 5.1, r),
+    "Invalid call to Weibull with shape = -1.1, scale = 5.1")
+  expect_error(
+    monty_random_weibull(5.1, -1.1, r),
+    "Invalid call to Weibull with shape = 5.1, scale = -1.1")
+})