16 add the lognormal distribution to supported distribution classes

DESCRIPTION

@@ -33,6 +33,7 @@ Imports:
     rstantools (>= 2.4.0),
     dplyr,
     magrittr,
+    matrixStats,
     purrr,
     stringr,
     R6 (>= 2.4.1)

NAMESPACE

@@ -8,6 +8,7 @@ export(check_numeric)
 export(distribution.binomial)
 export(distribution.exponential)
 export(distribution.gamma)
+export(distribution.lognormal)
 export(distribution.mixture)
 export(distribution.negative_binomial)
 export(distribution.normal)
@@ -36,6 +37,7 @@ importFrom(data.table,set)
 importFrom(data.table,setnames)
 importFrom(dplyr,if_else)
 importFrom(magrittr,"%>%")
+importFrom(matrixStats,colMaxs)
 importFrom(purrr,map_chr)
 importFrom(purrr,map_dbl)
 importFrom(rstan,sampling)

R/distribution_continuous.R

@@ -701,4 +701,137 @@ distribution.continuous.normal.class <- R6.class(
 distribution.normal <- function( mean, sd ){
   distribution.continuous.normal.class$new( mean = mean,
                                             sd   = sd )
+}
+
+################################################################################/
+#  distribution.continuous.lognormal
+################################################################################/
+#' Class: `distribution.continuous.lognormal.class`
+#' @description Derived class for a lognormally-distributed random variable.
+#'
+#' @param meanlog    the mean of log(X)
+#' @param sdlog      the standard deviation of log(X)
+#' @param x          vector of quantiles.
+#' @param q          vector of quantiles.
+#' @param p          vector of probabilities.
+#' @param n          number of observations. If `length( n ) > 1`, the length is
+#'   taken to be the number required.
+#' @param log        logical; if TRUE, probabilities p are given as `log(p)`.
+#' @param log.p      logical; if TRUE, probabilities p are given as `log(p)`.
+#' @param lower.tail logical; if TRUE (default), probabilities are 
+#'   \eqn{P[ X \leq x ]}, otherwise, \eqn{P[X>x]}.
+#' 
+#' @field interfaces The list of available class interfaces
+#' @field mean      the mean of the distribution
+#' @field sd        the standard deviation of the distribution
+#' @field var       the variance of the distribution
+distribution.continuous.lognormal.class <- R6.class(
+  classname = "distribution.continuous.lognormal.class",
+  inherit   = distribution.continuous.class,
+  private   = list(
+    .name    = "lognormal",
+    .param_names = c( "meanlog",  "sdlog" ),
+    .check_params = function( params ){
+      # Check that params contains all elements of private$.param_names
+      super$.check_params( params )
+
+      if ( !is.numeric( params$meanlog ) )
+        stop( "`params$meanlog` must be a numeric value.")
+      if ( params$sdlog < 0 || !is.numeric( params$sdlog ) )
+        stop( "`params$sdlog` must be a non-negative numeric value.")
+
+      return( NULL )
+    }
+  ),
+  public = list(
+    ############################################################################/
+    # initialize
+    ############################################################################/
+    #' @description Create a new object of class
+    #'   `distribution.continuous.normal.class`
+    initialize = function( meanlog, sdlog ){
+      super$initialize( support = c( 0, Inf ) )
+      self$params <- list( meanlog = meanlog,
+                           sdlog   = sdlog )
+    },
+    ############################################################################/
+    # density
+    ############################################################################/
+    #' @description Density function for a lognormal random variable with mean log(X)
+    #'   `$params$meanlog` and standard deviation log(X) `$params$sdlog`.
+    d = function( x, log = FALSE ){
+      stats::dlnorm( x, meanlog = private$.params$meanlog, sdlog = private$.params$sdlog,
+                    log = log )
+    },
+    ############################################################################/
+    # distribution function
+    ############################################################################/
+    #' @description Cumulative density function for a lognormal random variable
+    #'   with mean log(X) `$params$meanlog` and standard deviation log(X) `$params$sdlog`.
+    p = function( q, lower.tail = TRUE, log.p = FALSE ){
+      stats::plnorm( q, meanlog = private$.params$meanlog, sdlog = private$.params$sdlog,
+                    lower.tail = lower.tail, log.p = log.p )
+    },
+    ############################################################################/
+    # quantile function
+    ############################################################################/
+    #' @description Quantile function for a lognormal random variable with mean
+    #'  log(X)  `$params$meanlog` and standard deviation log(X) `$params$sdlog`.
+    q = function( p, lower.tail = TRUE, log.p = FALSE ){
+      stats::qlnorm( p, meanlog = private$.params$meanlog, sdlog = private$.params$sdlog,
+                    lower.tail = lower.tail, log.p = log.p )
+    },
+    ############################################################################/
+    # random deviates
+    ############################################################################/
+    #' @description Generates random deviates for a lognormal random variable with
+    #'   mean log(X) `$params$meanlog` and standard deviation logx(X) `$params$sdlog` .
+    r = function( n ){
+      stats::rlnorm( n, meanlog = private$.params$meanlog, sdlog = private$.params$sdlog )
+    }
+  ),
+  active = list(
+    ############################################################################/
+    # mean
+    ############################################################################/
+    mean = function( val ){
+      if( !missing( val ) )
+        stop( "cannot set `$mean`" )
+      return( exp( private$.params$meanlog + 0.5 * private$.params$sdlog^2 ) )
+    },
+    ############################################################################/
+    # standard deviation
+    ############################################################################/
+    sd = function( val ){
+      if( !missing( val ) )
+        stop( "cannot set `$sd`" )
+      return( sqrt( self$var ) )
+    },
+    ############################################################################/
+    # variance
+    ############################################################################/
+    var = function( val ){
+      if( !missing( val ) )
+        stop( "cannot set `$var`" )
+      sigma2 <- private$.params$sdlog^2
+      mu     <- private$.params$meanlog
+      return( ( exp( sigma2) - 1 ) * exp( 2 * mu + sigma2 ) )
+    }
+  )
+)
+
+#' distribution.lognormal
+#' 
+#' Constructor function for an object of class `distribution.continuous.lognormal.class`
+#' 
+#' @param meanlog The mean of log(X)
+#' @param sdlog   The standard deviation log(X)
+#' 
+#' @returns An object of class [[distribution.continuous.lognormal.class]]
+#'
+#' @seealso [Mastiff-Distributions]
+#' @export
+distribution.lognormal <- function( meanlog, sdlog ){
+  distribution.continuous.lognormal.class$new( meanlog = meanlog,
+                                               sdlog   = sdlog )
 }

R/distribution_mixture.R

@@ -7,11 +7,22 @@
 #' Class: `distribution.mixture.class`
 #' @description Class to describe the mixture of distributions
 #'
+#' @param distributions list of distributions in the mixture
+#' @param weights       vector of weights of the distributions in the mixture
+#' @param x          vector of quantiles.
 #' @param q          vector of quantiles.
 #' @param p          vector of probabilities.
+#' @param n          number of observations. If `length( n ) > 1`, the length is
+#'   taken to be the number required.
+#' @param log        logical; if TRUE, probabilities p are given as `log(p)`.
 #' @param log.p      logical; if TRUE, probabilities p are given as `log(p)`.
 #' @param lower.tail logical; if TRUE (default), probabilities are \eqn{P[ X \leq x ]},
 #'   otherwise, \eqn{P[X>x]}.
+#'   
+#' @field n_distributions the number of distributions in the mixture
+#' @field distributions   the distributions in the mixture
+#' @field weights         the weights of the distributions in the mixture
+#' @field interfaces      the list of available class interfaces
 #' 
 #' @include R6_class.R
 #' @include distribution_R6_class.R
@@ -21,7 +32,16 @@ distribution.mixture.class <- R6.class(
   private   = list(
     .distributions = NULL,
     .n_distributions = NULL,
-    .weights = NULL
+    .weights = NULL,
+    .logExpAdd = function( ds, w ) {
+      eps = 1e-16
+      nw  = ncol( ds )
+      ns  = nrow( ds ) 
+      ds  <- t( ds ) + matrix( log( w + eps ), nrow = nw, ncol = ns ) 
+      max <- matrixStats::colMaxs( ds )
+      ds  <- ds - matrix( rep( max, each = nw ), nrow = nw, ncol = ns )  
+      return( log( colSums( exp( ds ) ) ) + max )
+    }
   ),
   public = list(
     ############################################################################/
@@ -46,9 +66,14 @@ distribution.mixture.class <- R6.class(
     #' @description Density function for a random variable of the mixture
     d = function( x, log = FALSE ){
       ds <- matrix( unlist( lapply( private$.distributions, 
-        function( d ) d$d( x ) ) ), 
+        function( d ) d$d( x, log = log ) ) ), 
         ncol = self$n_distributions )  
-      return( ( ds %*% self$weights )[, 1 ] )
+
+      if( log ) {
+        return( private$.logExpAdd( ds, self$weights ) )
+      } else {
+        return( ( ds %*% self$weights )[, 1 ] )
+      }
     },
     ##############################################################################/
     # cumulative distribution function
@@ -58,18 +83,19 @@ distribution.mixture.class <- R6.class(
       ps <- matrix( unlist( lapply( private$.distributions, 
         function( d ) d$p( q, lower.tail = lower.tail, log.p = log.p ) ) ), 
         ncol = self$n_distributions ) 
-      return( ( ps %*% self$weights )[, 1 ] )
+
+      if( log.p ) {
+        return( private$.logExpAdd( ps, self$weights ) )
+      } else {
+        return( ( ps %*% self$weights )[, 1 ] )
+      }
     },
     ##############################################################################/
     # quantile function
     ##############################################################################/
     #' @description Evaluates the quantile function of the mixture
     q = function( p, lower.tail = TRUE, log.p = FALSE ){
       super$q( p, lower.tail = lower.tail, log.p = log.p )      
-      # qs <- matrix( unlist( lapply( private$.distributions, 
-      #   function( d ) d$q( p, lower.tail = lower.tail, log.p = log.p ) ) ), 
-      #   ncol = self$n_distributions  ) 
-      # return( ( qs %*% self$weights )[, 1 ] )
     },
     ############################################################################/
     # random deviates
@@ -90,20 +116,15 @@ distribution.mixture.class <- R6.class(
     }
   ),
   active = list(
-    support = function( val ){
-      private$.staticReturn( val, "support" )
-    },
     ############################################################################/
     # n_distributions
     ############################################################################/
-    #' @description Number of distributions in the mixture
     n_distributions = function( val ){
       private$.staticReturn( val, "n_distributions" )
     },
     ############################################################################/
     # distributions
     ############################################################################/
-    #' @description The distributions of the mixture
     distributions = function( val ){
       if( missing( val ) )
         return( private$.staticReturn( val, "distributions" ) )
@@ -121,7 +142,6 @@ distribution.mixture.class <- R6.class(
     ############################################################################/
     # weights
     ############################################################################/
-    #' @description The weights of each distribution in the mixture
     weights = function( new = NA ){
       if( length( new) == 1 ) 
         if( is.na( new ) ) {