Add log-normal

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: monty
 Title: Monte Carlo Models
-Version: 0.3.14
+Version: 0.3.15
 Authors@R: c(person("Rich", "FitzJohn", role = c("aut", "cre"),
                     email = "rich.fitzjohn@gmail.com"),
              person("Wes", "Hinsley", role = "aut"),

NAMESPACE

@@ -43,6 +43,7 @@ export(monty_random_exponential_rate)
 export(monty_random_gamma_rate)
 export(monty_random_gamma_scale)
 export(monty_random_hypergeometric)
+export(monty_random_log_normal)
 export(monty_random_n_beta)
 export(monty_random_n_beta_binomial_ab)
 export(monty_random_n_beta_binomial_prob)
@@ -53,6 +54,7 @@ export(monty_random_n_exponential_rate)
 export(monty_random_n_gamma_rate)
 export(monty_random_n_gamma_scale)
 export(monty_random_n_hypergeometric)
+export(monty_random_n_log_normal)
 export(monty_random_n_negative_binomial_mu)
 export(monty_random_n_negative_binomial_prob)
 export(monty_random_n_normal)

R/dsl-distributions.R

@@ -168,6 +168,17 @@ distr_hypergeometric <- distribution(
     mean = quote(k * n1 / (n1 + n2))),
   cpp = list(density = "hypergeometric", sample = "hypergeometric"))
 
+distr_log_normal <- distribution(
+  name = "LogNormal",
+  density = function(x, meanlog, sdlog) dlnorm(x, meanlog, sdlog, log = TRUE),
+  domain = c(0, Inf),
+  expr = list(
+    density = quote(-(log(x) - meanlog)^2 / (2 * sdlog^2) - log(2 * pi) / 2 - 
+                      log(sdlog) - log(x)),
+    mean = quote(exp(meanlog + sdlog^2 / 2))),
+  sample = "monty_random_log_normal",
+  cpp = list(density = "log_normal", sample = "log_normal"))
+
 distr_negative_binomial_prob <- distribution(
   name = "NegativeBinomial",
   variant = "prob",
@@ -269,6 +280,7 @@ dsl_distributions <- local({
     distr_gamma_rate,
     distr_gamma_scale,
     distr_hypergeometric,
+    distr_log_normal,
     distr_negative_binomial_prob,
     distr_negative_binomial_mu,
     distr_normal,

R/random.R

@@ -637,3 +637,27 @@ monty_random_weibull <- function(shape, scale, state) {
 monty_random_n_weibull <- function(n_samples, shape, scale, state) {
   cpp_monty_random_n_weibull(n_samples, shape, scale, state)
 }
+
+
+##' Sample from a log-normal distribution
+##'
+##' @title Sample from log-normal
+##'
+##' @param meanlog The mean of the distribution on the log scale
+##'
+##' @param sdlog The standard deviation of the distribution on the log scale
+##'
+##' @inheritParams monty_random_real
+##' @inherit monty_random_real return
+##'
+##' @export
+monty_random_log_normal <- function(meanlog, sdlog, state) {
+  cpp_monty_random_log_normal(meanlog, sdlog, state)
+}
+
+
+##' @export
+##' @rdname monty_random_log_normal
+monty_random_n_log_normal <- function(n_samples, meanlog, sdlog, state) {
+  cpp_monty_random_n_log_normal(n_samples, meanlog, sdlog, state)
+}

_pkgdown.yml

@@ -115,6 +115,8 @@ reference:
     - monty_random_n_hypergeometric
     - monty_random_truncated_normal
     - monty_random_n_truncated_normal
+    - monty_random_log_normal
+    - monty_random_n_log_normal
     - monty_random_weibull
     - monty_random_n_weibull
   - subsection: Legacy interface

inst/include/monty/random/density.hpp

@@ -293,5 +293,12 @@ __host__ __device__ T weibull(T x, T shape, T scale, bool log) {
   return maybe_log(ret, log);
 }
 
+template <typename T>
+__host__ __device__ T log_normal(T x, T mulog, T sdlog, bool log) {
+  const auto d = normal(monty::math::log(x), mulog, sdlog);
+
+  return log ? (d - monty::math::log(x)) : (monty::math::exp(d) / x);
+}
+
 }
 }

inst/include/monty/random/log_normal.hpp

@@ -0,0 +1,40 @@
+#pragma once
+
+#include <cmath>
+
+#include "monty/random/generator.hpp"
+#include "monty/random/normal.hpp"
+
+namespace monty {
+namespace random {
+
+namespace {
+
+template <typename real_type>
+void log_normal_validate(real_type meanlog, real_type sdlog) {
+  if (!std::isfinite(meanlog) || !std::isfinite(sdlog) || sdlog <= 0) {
+    char buffer[256];
+    snprintf(buffer, 256,
+             "Invalid call to log_normal with meanlog = %g, sdlog = %g",
+             meanlog, sdlog);
+    monty::utils::fatal_error(buffer);
+  }
+}
+
+}
+
+template <typename real_type, typename rng_state_type>
+real_type log_normal(rng_state_type& rng_state, real_type meanlog, real_type sdlog) {
+#ifdef __CUDA_ARCH__
+  static_assert("log_normal() not implemented for GPU targets");
+#endif
+    log_normal_validate(meanlog, sdlog);
+
+    if (rng_state.deterministic) {
+      return monty::math::exp(meanlog + monty::math::pow(sdlog, 2) / 2);
+    }
+    return monty::math::exp(normal(rng_state, meanlog, sdlog));
+}
+
+}
+}

inst/include/monty/random/random.hpp

@@ -10,6 +10,7 @@
 #include "monty/random/exponential.hpp"
 #include "monty/random/gamma.hpp"
 #include "monty/random/hypergeometric.hpp"
+#include "monty/random/log_normal.hpp"
 #include "monty/random/multinomial.hpp"
 #include "monty/random/negative_binomial.hpp"
 #include "monty/random/truncated_normal.hpp"

src/random.cpp

@@ -639,6 +639,26 @@ cpp11::doubles cpp_monty_random_n_weibull(size_t n_samples,
                                  shape, scale, "shape", "scale");
 }
 
+// log_normal
+[[cpp11::register]]
+cpp11::doubles cpp_monty_random_log_normal(cpp11::doubles meanlog,
+                                           cpp11::doubles sdlog,
+                                           cpp11::sexp ptr) {
+  const auto fn = [](auto& state, auto meanlog, auto sdlog) { return monty::random::log_normal<double>(state, meanlog, sdlog); };
+  return monty_random_sample_1_2(fn, ptr, "log_normal",
+                                 meanlog, sdlog, "meanlog", "sdlog");
+}
+
+[[cpp11::register]]
+cpp11::doubles cpp_monty_random_n_log_normal(size_t n_samples,
+                                             cpp11::doubles meanlog,
+                                             cpp11::doubles sdlog,
+                                             cpp11::sexp ptr) {
+  const auto fn = [](auto& state, auto meanlog, auto sdlog) { return monty::random::log_normal<double>(state, meanlog, sdlog); };
+  return monty_random_sample_n_2(fn, n_samples, ptr, "log_normal",
+                                 meanlog, sdlog, "meanlog", "sdlog");
+}
+
 //// 3-arg functions
 
 // beta_binomial_prob

tests/testthat/test-rng.R

@@ -1404,3 +1404,44 @@ test_that("weibull random numbers prevent bad inputs", {
     monty_random_weibull(5.1, -1.1, r),
     "Invalid call to Weibull with shape = 5.1, scale = -1.1")
 })
+
+
+test_that("can draw log-normal random numbers", {
+  meanlog <- 1.5
+  sdlog <- 0.5
+  n <- 10000000
+
+  ans1 <- 
+    monty_random_n_log_normal(n, meanlog, sdlog, monty_rng_create(seed = 1))
+  ans2 <- 
+    monty_random_n_log_normal(n, meanlog, sdlog, monty_rng_create(seed = 1))
+  expect_identical(ans1, ans2)
+
+  expect_equal(mean(ans1), exp(meanlog + sdlog^2 / 2), tolerance = 1e-3)
+  true_var <- (exp(sdlog^2) - 1) * exp(2 * meanlog + sdlog^2)
+  expect_equal(var(ans1), true_var, tolerance = 1e-2)
+})
+
+
+test_that("deterministic log-normal returns mean", {
+  n_reps <- 10
+  meanlog <- as.numeric(sample(seq(-10, 10), n_reps, replace = TRUE))
+  sdlog <- 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_log_normal, meanlog, sdlog, MoreArgs = list(rng)),
+    exp(meanlog + sdlog^2 / 2))
+  expect_equal(monty_rng_state(rng), state)
+})
+
+
+test_that("log-normal random numbers prevent bad inputs", {
+  r <- monty_rng_create(seed = 1)
+
+  expect_error(
+    monty_random_log_normal(1.1, -5.1, r),
+    "Invalid call to log_normal with meanlog = 1.1, sdlog = -5.1")
+})