Merge pull request #717 from HugoMVale/stats-distribution-exponential

Extension of PR#679 to `stdlib_stats_distribution_exponential`

Author: jvdp1

doc/specs/stdlib_stats_distribution_exponential.md

@@ -14,34 +14,37 @@ Experimental
 
 ### Description
 
-An exponential distribution is the distribution of time between events in a Poisson point process. The inverse scale parameter `lambda` specifies the average time between events, also called the rate of events.
+An exponential distribution is the distribution of time between events in a Poisson point process. The inverse scale parameter `lambda` specifies the average time between events ($\lambda$), also called the rate of events.
 
-Without argument the function returns a random sample from the standard exponential distribution `E(1)` with `lambda = 1`.
+Without argument, the function returns a random sample from the standard exponential distribution $E(\lambda=1)$.
 
-With a single argument, the function returns a random sample from the exponential distribution `E(lambda)`.
+With a single argument, the function returns a random sample from the exponential distribution $E(\lambda=\text{lambda})$.
 For complex arguments, the real and imaginary parts are sampled independently of each other.
 
-With two arguments the function returns a rank one array of exponentially distributed random variates.
+With two arguments, the function returns a rank-1 array of exponentially distributed random variates.
 
-Note: the algorithm used for generating exponetial random variates is fundamentally limited to double precision. Ref.: Marsaglia, G. & Tsang, W.W. (2000) `The ziggurat method for generating random variables', J. Statist. Software, v5(8).
+@note
+The algorithm used for generating exponential random variates is fundamentally limited to double precision.[^1]
 
 ### Syntax
 
 `result = [[stdlib_stats_distribution_exponential(module):rvs_exp(interface)]]([lambda] [[, array_size]])`
 
 ### Class
 
-Function
+Elemental function
 
 ### Arguments
 
-`lambda`: optional argument has `intent(in)` and is a scalar of type `real` or `complex`. The value of `lambda` has to be non-negative.
+`lambda`: optional argument has `intent(in)` and is a scalar of type `real` or `complex`.
+If `lambda` is `real`, its value must be positive. If `lambda` is `complex`, both the real and imaginary components must be positive.
 
 `array_size`: optional argument has `intent(in)` and is a scalar of type `integer` with default kind.
 
 ### Return value
 
-The result is a scalar or rank one array with a size of `array_size`, and has the same type of `lambda`.
+The result is a scalar or rank-1 array with a size of `array_size`, and the same type as `lambda`.
+If `lambda` is non-positive, the result is `NaN`.
 
 ### Example
 
@@ -57,15 +60,13 @@ Experimental
 
 ### Description
 
-The probability density function (pdf) of the single real variable exponential distribution:
+The probability density function (pdf) of the single real variable exponential distribution is:
 
-$$f(x)=\begin{cases} \lambda e^{-\lambda x} &x\geqslant 0 \\\\ 0 &x< 0\end{}$$
+$$f(x)=\begin{cases} \lambda e^{-\lambda x} &x\geqslant 0 \\\\ 0 &x< 0\end{cases}$$
 
-For a complex variable (x + y i) with independent real x and imaginary y parts, the joint probability density function
-is the product of the corresponding marginal pdf of real and imaginary pdf (for more details, see
-"Probability and Random Processes with Applications to Signal Processing and Communications", 2nd ed., Scott L. Miller and Donald Childers, 2012, p.197):
+For a complex variable $z=(x + y i)$ with independent real $x$ and imaginary $y$ parts, the joint probability density function is the product of the corresponding real and imaginary marginal pdfs:[^2]
 
-$$f(x+\mathit{i}y)=f(x)f(y)=\begin{cases} \lambda_{x} \lambda_{y} e^{-(\lambda_{x} x + \lambda_{y} y)} &x\geqslant 0, y\geqslant 0 \\\\ 0 &otherwise\end{}$$
+$$f(x+\mathit{i}y)=f(x)f(y)=\begin{cases} \lambda_{x} \lambda_{y} e^{-(\lambda_{x} x + \lambda_{y} y)} &x\geqslant 0, y\geqslant 0 \\\\ 0 &\text{otherwise}\end{cases}$$
 
 ### Syntax
 
@@ -80,12 +81,13 @@ Elemental function
 `x`: has `intent(in)` and is a scalar of type `real` or `complex`.
 
 `lambda`: has `intent(in)` and is a scalar of type `real` or `complex`.
+If `lambda` is `real`, its value must be positive. If `lambda` is `complex`, both the real and imaginary components must be positive.
 
 All arguments must have the same type.
 
 ### Return value
 
-The result is a scalar or an array, with a shape conformable to arguments, and has the same type of input arguments.
+The result is a scalar or an array, with a shape conformable to the arguments, and the same type as the input arguments. If `lambda` is non-positive, the result is `NaN`.
 
 ### Example
 
@@ -103,13 +105,11 @@ Experimental
 
 Cumulative distribution function (cdf) of the single real variable exponential distribution:
 
-$$F(x)=\begin{cases}1 - e^{-\lambda x} &x\geqslant 0 \\\\ 0 &x< 0\end{}$$
+$$F(x)=\begin{cases}1 - e^{-\lambda x} &x\geqslant 0 \\\\ 0 &x< 0\end{cases}$$
 
-For a complex variable (x + y i) with independent real x and imaginary y parts, the joint cumulative distribution
-function is the product of corresponding marginal cdf of real and imaginary cdf (for more details, see
-"Probability and Random Processes with Applications to Signal Processing and Communications", 2nd ed., Scott L. Miller and Donald Childers, 2012, p.197):
+For a complex variable $z=(x + y i)$ with independent real $x$ and imaginary $y$ parts, the joint cumulative distribution function is the product of corresponding real and imaginary marginal cdfs:[^2]
 
-$$F(x+\mathit{i}y)=F(x)F(y)=\begin{cases} (1 - e^{-\lambda_{x} x})(1 - e^{-\lambda_{y} y}) &x\geqslant 0, \;\; y\geqslant 0 \\\\ 0 &otherwise \end{}$$
+$$F(x+\mathit{i}y)=F(x)F(y)=\begin{cases} (1 - e^{-\lambda_{x} x})(1 - e^{-\lambda_{y} y}) &x\geqslant 0, \;\; y\geqslant 0 \\\\ 0 & \text{otherwise} \end{cases}$$
 
 ### Syntax
 
@@ -124,15 +124,20 @@ Elemental function
 `x`: has `intent(in)` and is a scalar of type `real` or `complex`.
 
 `lambda`: has `intent(in)` and is a scalar of type `real` or `complex`.
+If `lambda` is `real`, its value must be positive. If `lambda` is `complex`, both the real and imaginary components must be positive.
 
 All arguments must have the same type.
 
 ### Return value
 
-The result is a scalar or an array, with a shape conformable to arguments, and has the same type of input arguments.
+The result is a scalar or an array, with a shape conformable to the arguments, and the same type as the input arguments. If `lambda` is non-positive, the result is `NaN`.
 
 ### Example
 
 ```fortran
 {!example/stats_distribution_exponential/example_exponential_cdf.f90!}
 ```
+
+[^1] Marsaglia, George, and Wai Wan Tsang. "The ziggurat method for generating random variables." _Journal of statistical software_ 5 (2000): 1-7.
+
+[^2] Miller, Scott, and Donald Childers. _Probability and random processes: With applications to signal processing and communications_. Academic Press, 2012 (p. 197).

example/stats_distribution_exponential/example_exponential_cdf.f90

@@ -4,38 +4,61 @@ program example_exponential_cdf
                                                    rexp => rvs_exp
 
   implicit none
-  real :: x(2, 3, 4), a(2, 3, 4)
+  real, dimension(2, 3, 4) :: x, lambda
+  real :: xsum
   complex :: scale
-  integer :: seed_put, seed_get
+  integer :: seed_put, seed_get, i
 
   seed_put = 1234567
   call random_seed(seed_put, seed_get)
 
-  print *, exp_cdf(1.0, 1.0)  ! a standard exponential cumulative at 1.0
+  ! standard exponential cumulative distribution at x=1.0
+  print *, exp_cdf(1.0, 1.0)
+  ! 0.632120550
 
-! 0.632120550
+  ! cumulative distribution at x=2.0 with lambda=2
+  print *, exp_cdf(2.0, 2.0)
+  ! 0.981684387
 
-  print *, exp_cdf(2.0, 2.0) ! a cumulative at 2.0 with lambda=2
-
-! 0.981684387
+  ! cumulative distribution at x=2.0 with lambda=-1.0 (out of range)
+  print *, exp_cdf(2.0, -1.0) 
+  ! NaN
 
+   ! standard exponential random variates array
   x = reshape(rexp(0.5, 24), [2, 3, 4])
-! standard exponential random variates array
-  a(:, :, :) = 0.5
-  print *, exp_cdf(x, a)  ! a rank 3 array of standard exponential cumulative
 
-!  8.57694745E-02  9.70223546E-02  1.52170658E-02  2.95336246E-02
-!  0.107568979  0.196659625  2.97447443E-02  0.366151094  0.163051903
-!  3.36527228E-02  0.385267735  0.428375721  0.103964329  0.147119939
-!  0.192206264  0.330693483  0.179247737  2.92580128E-02  0.332765043
-!  0.472417951  0.500440359  8.56802464E-02  8.72612000E-03  3.55126858E-02
+  ! a rank-3 exponential cumulative distribution
+  lambda(:, :, :) = 0.5
+  print *, exp_cdf(x, lambda)
+  ! 0.301409245  0.335173965  5.94930053E-02  0.113003314
+  ! 0.365694344  0.583515942  0.113774836     0.838585377
+  ! 0.509324908  0.127967060  0.857194781     0.893231630
+  ! 0.355383813  0.470882893  0.574203610     0.799321830
+  ! 0.546216846  0.111995399  0.801794767     0.922525287
+  ! 0.937719882  0.301136374  3.44503522E-02  0.134661376 
+
+  ! cumulative distribution array where lambda<=0.0 for certain elements 
+  print *, exp_cdf([1.0, 1.0, 1.0], [1.0, 0.0, -1.0])
+  ! 0.632120550  NaN NaN
+
+  ! `cdf_exp` is pure and, thus, can be called concurrently 
+  xsum = 0.0
+  do concurrent (i=1:size(x,3))
+    xsum = xsum + sum(exp_cdf(x(:,:,i), lambda(:,:,i)))
+  end do
+  print *, xsum
+  ! 11.0886612
 
+  ! complex exponential cumulative distribution at (0.5,0.5) with real part of
+  ! lambda=0.5 and imaginary part of lambda=1.0
   scale = (0.5, 1.0)
   print *, exp_cdf((0.5, 0.5), scale)
-!complex exponential cumulative distribution at (0.5,0.5) with real part of
-!lambda=0.5 and imaginary part of lambda=1.0
+  ! 8.70351046E-02
 
-! 8.70351046E-02
+  ! As above, but with lambda%im < 0 
+  scale = (1.0, -2.0)
+  print *, exp_cdf((1.5, 1.0), scale)
+  ! NaN
 
 end program example_exponential_cdf
 

example/stats_distribution_exponential/example_exponential_pdf.f90

@@ -1,39 +1,62 @@
 program example_exponential_pdf
   use stdlib_random, only: random_seed
   use stdlib_stats_distribution_exponential, only: exp_pdf => pdf_exp, &
-                                                   rexp => rvs_exp
+                                                    rexp => rvs_exp
 
   implicit none
-  real :: x(2, 3, 4), a(2, 3, 4)
+  real, dimension(2, 3, 4) :: x, lambda
+  real :: xsum
   complex :: scale
-  integer :: seed_put, seed_get
+  integer :: seed_put, seed_get, i
 
   seed_put = 1234567
   call random_seed(seed_put, seed_get)
 
-  print *, exp_pdf(1.0, 1.0) !a probability density at 1.0 in standard expon
-
-! 0.367879450
-
-  print *, exp_pdf(2.0, 2.0) !a probability density at 2.0 with lambda=2.0
-
-! 3.66312787E-02
-
-  x = reshape(rexp(0.5, 24), [2, 3, 4]) ! standard expon random variates array
-  a(:, :, :) = 0.5
-  print *, exp_pdf(x, a)     ! a rank 3 standard expon probability density
-
-!  0.457115263  0.451488823  0.492391467  0.485233188  0.446215510
-!  0.401670188  0.485127628  0.316924453  0.418474048  0.483173639
-!  0.307366133  0.285812140  0.448017836  0.426440030  0.403896868
-!  0.334653258  0.410376132  0.485370994  0.333617479  0.263791025
-!  0.249779820  0.457159877  0.495636940  0.482243657
-
+  ! probability density at x=1.0 in standard exponential
+  print *, exp_pdf(1.0, 1.0)
+  ! 0.367879450
+
+  ! probability density at x=2.0 with lambda=2.0
+  print *, exp_pdf(2.0, 2.0) 
+  ! 3.66312787E-02
+
+  ! probability density at x=2.0 with lambda=-1.0 (out of range)
+  print *, exp_pdf(2.0, -1.0) 
+  ! NaN
+
+  ! standard exponential random variates array  
+  x = reshape(rexp(0.5, 24), [2, 3, 4])
+
+  ! a rank-3 exponential probability density
+  lambda(:, :, :) = 0.5
+  print *, exp_pdf(x, lambda)
+  ! 0.349295378      0.332413018     0.470253497     0.443498343      0.317152828
+  ! 0.208242029      0.443112582     8.07073265E-02  0.245337561      0.436016470
+  ! 7.14025944E-02   5.33841923E-02  0.322308093     0.264558554      0.212898195
+  ! 0.100339092      0.226891592     0.444002301     9.91026312E-02   3.87373678E-02
+  ! 3.11400592E-02   0.349431813     0.482774824     0.432669312     
+
+  ! probability density array where lambda<=0.0 for certain elements 
+  print *, exp_pdf([1.0, 1.0, 1.0], [1.0, 0.0, -1.0])
+  ! 0.367879450  NaN NaN
+
+  ! `pdf_exp` is pure and, thus, can be called concurrently 
+  xsum = 0.0
+  do concurrent (i=1:size(x,3))
+    xsum = xsum + sum(exp_pdf(x(:,:,i), lambda(:,:,i)))
+  end do
+  print *, xsum
+  ! 6.45566940
+
+  ! complex exponential probability density function at (1.5,1.0) with real part
+  ! of lambda=1.0 and imaginary part of lambda=2.0
   scale = (1.0, 2.)
   print *, exp_pdf((1.5, 1.0), scale)
-! a complex expon probability density function at (1.5,1.0) with real part
-!of lambda=1.0 and imaginary part of lambda=2.0
+  ! 6.03947677E-02
 
-! 6.03947677E-02
+  ! As above, but with lambda%re < 0 
+  scale = (-1.0, 2.)
+  print *, exp_pdf((1.5, 1.0), scale)
+  ! NaN
 
 end program example_exponential_pdf