Merge pull request #718 from HugoMVale/stats-distribution-normal

Improve doc + examples of `stdlib_stats_distribution_normal`

Author: jvdp1

doc/specs/stdlib_stats_distribution_normal.md

@@ -14,15 +14,16 @@ Experimental
 
 ### Description
 
-A normal continuous random variate distribution, also known as Gaussian, or Gauss or Laplace-Gauss distribution. The location `loc` specifies the mean or expectation. The `scale` specifies the standard deviation. 
+A normal continuous random variate distribution, also known as Gaussian, or Gauss or Laplace-Gauss distribution. The location `loc` specifies the mean or expectation ($\mu$). The `scale` specifies the standard deviation ($\sigma$).
 
-Without argument the function returns a standard normal distributed random variate N(0,1).
+Without argument, the function returns a standard normal distributed random variate $N(0,1)$.
 
-With two arguments, the function returns a normal distributed random variate N(loc, scale^2). For complex arguments, the real and imaginary parts are independent of each other.
+With two arguments, the function returns a normal distributed random variate $N(\mu=\text{loc}, \sigma^2=\text{scale}^2)$. For complex arguments, the real and imaginary parts are independent of each other.
 
-With three arguments, the function returns a rank one array of normal distributed random variates.
+With three arguments, the function returns a rank-1 array of normal distributed random variates.
 
-Note: the algorithm used for generating normal random variates is fundamentally limited to double precision.
+@note
+The algorithm used for generating exponential random variates is fundamentally limited to double precision.[^1]
 
 ### Syntax
 
@@ -36,15 +37,15 @@ Elemental function (passing both `loc` and `scale`).
 
 `loc`: optional argument has `intent(in)` and is a scalar of type `real` or `complex`.
 
-`scale`: optional argument has `intent(in)` and is a scalar of type `real` or `complex`.
+`scale`: optional argument has `intent(in)` and is a positive scalar of type `real` or `complex`.
 
 `array_size`: optional argument has `intent(in)` and is a scalar of type `integer`.
 
 `loc` and `scale` arguments must be of the same type.
 
 ### Return value
 
-The result is a scalar or rank one array, with a size of `array_size`, and as the same type of `scale` and `loc`.
+The result is a scalar or rank-1 array, with a size of `array_size`, and the same type as `scale` and `loc`. If `scale` is non-positive, the result is `NaN`.
 
 ### Example
 
@@ -62,11 +63,11 @@ Experimental
 
 The probability density function (pdf) of the single real variable normal distribution:
 
-$$f(x) = \frac{1}{\sigma \sqrt{2}} e^{-\frac{1}{2}(\frac{x-\mu}{\sigma})^{2}}$$
+$$f(x) = \frac{1}{\sigma \sqrt{2}} \exp{\left[-\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^{2}\right]}$$
 
-For complex varible (x + y i) with independent real x and imaginary y parts, the joint probability density function is the product of corresponding marginal pdf of real and imaginary pdf (ref. "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 varible $z=(x + y i)$ with independent real $x$ and imaginary $y$ parts, the joint probability density function is the product of the the corresponding real and imaginary marginal pdfs:[^2]
 
-$$f(x + y \mathit{i}) = f(x) f(y) = \frac{1}{2\sigma_{x}\sigma_{y}} e^{-\frac{1}{2}[(\frac{x-\mu}{\sigma_{x}})^{2}+(\frac{y-\nu}{\sigma_{y}})^{2}]}$$
+$$f(x + y \mathit{i}) = f(x) f(y) = \frac{1}{2\sigma_{x}\sigma_{y}} \exp{\left[-\frac{1}{2}\left(\left(\frac{x-\mu_x}{\sigma_{x}}\right)^{2}+\left(\frac{y-\mu_y}{\sigma_{y}}\right)^{2}\right)\right]}$$
 
 ### Syntax
 
@@ -82,13 +83,13 @@ Elemental function
 
 `loc`: has `intent(in)` and is a scalar of type `real` or `complex`.
 
-`scale`: has `intent(in)` and is a scalar of type `real` or `complex`.
+`scale`: has `intent(in)` and is a positive scalar of type `real` or `complex`.
 
 All three arguments must have the same type.
 
 ### Return value
 
-The result is a scalar or an array, with a shape conformable to arguments, and as 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 `scale` is non-positive, the result is `NaN`.
 
 ### Example
 
@@ -106,11 +107,13 @@ Experimental
 
 Cumulative distribution function of the single real variable normal distribution:
 
-$$F(x)=\frac{1}{2}\left [ 1+erf(\frac{x-\mu}{\sigma \sqrt{2}}) \right ]$$
+$$F(x) = \frac{1}{2}\left [ 1+\text{erf}\left(\frac{x-\mu}{\sigma \sqrt{2}}\right) \right ]$$
 
-For the 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 (ref. "Probability and Random Processes with Applications to Signal Processing and Communications", 2nd ed., Scott L. Miller and Donald Childers, 2012, p.197):
+For the complex variable $z=(x + y i)$ with independent real $x$ and imaginary $y$ parts, the joint cumulative distribution function is the product of the corresponding real and imaginary marginal cdfs:[^2]
 
-$$F(x+y\mathit{i})=F(x)F(y)=\frac{1}{4} [1+erf(\frac{x-\mu}{\sigma_{x} \sqrt{2}})] [1+erf(\frac{y-\nu}{\sigma_{y} \sqrt{2}})]$$
+$$ F(x+y\mathit{i})=F(x)F(y)=\frac{1}{4} \
+\left[ 1+\text{erf}\left(\frac{x-\mu_x}{\sigma_x \sqrt{2}}\right) \right] \
+\left[ 1+\text{erf}\left(\frac{y-\mu_y}{\sigma_y \sqrt{2}}\right) \right] $$
 
 ### Syntax
 
@@ -126,16 +129,20 @@ Elemental function
 
 `loc`: has `intent(in)` and is a scalar of type `real` or `complex`.
 
-`scale`: has `intent(in)` and is a scalar of type `real` or `complex`.
+`scale`: has `intent(in)` and is a positive scalar of type `real` or `complex`.
 
 All three arguments must have the same type.
 
 ### Return value
 
-The result is a scalar or an array, with a shape conformable to arguments, as 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 `scale` is non-positive, the result is `NaN`.
 
 ### Example
 
 ```fortran
 {!example/stats_distribution_normal/example_norm_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_normal/CMakeLists.txt

@@ -1,3 +1,3 @@
 ADD_EXAMPLE(normal_pdf)
 ADD_EXAMPLE(normal_rvs)
-ADD_EXAMPLE(norm_cdf)
+ADD_EXAMPLE(normal_cdf)

example/stats_distribution_normal/example_norm_cdf.f90

@@ -0,0 +1,60 @@
+program example_normal_cdf
+  use stdlib_random, only: random_seed
+  use stdlib_stats_distribution_normal, only: norm_cdf => cdf_normal, &
+                                              norm => rvs_normal
+
+  implicit none
+  real, dimension(2, 3, 4) :: x, mu, sigma
+  real :: xsum
+  complex :: loc, scale
+  integer :: seed_put, seed_get, i
+
+  seed_put = 1234567
+  call random_seed(seed_put, seed_get)
+
+  ! standard normal cumulative probability at x=0.0
+  print *, norm_cdf(0.0, 0.0, 1.0)
+  ! 0.500000000
+
+  ! cumulative probability at x=2.0 with mu=-1.0 sigma=2.0
+  print *, norm_cdf(2.0, -1.0, 2.0)
+  ! 0.933192849
+
+  ! cumulative probability at x=1.0 with mu=1.0 and sigma=-1.0 (out of range) 
+  print *, norm_cdf(1.0, 1.0, -1.0)
+  ! NaN
+
+  ! standard normal random variates array
+  x = reshape(norm(0.0, 1.0, 24), [2, 3, 4])
+
+  ! standard normal cumulative array
+  mu(:, :, :) = 0.0
+  sigma(:, :, :) = 1.0
+  print *, norm_cdf(x, mu, sigma)        
+  !  0.713505626     0.207069695  0.486513376  0.424511284  0.587328553
+  !  0.335559726     0.401470929  0.806552052  0.866687536  0.371323735
+  !  0.866228044     0.898046613  0.198435277  0.141147852  0.681565762
+  !  0.206268221     0.627057910  0.580759525  0.190364420  7.27325380E-02
+  !  7.08068311E-02  0.728241026  0.522919059  0.390097380
+
+   ! cumulative probability array where sigma<=0.0 for certain elements 
+  print *, norm_cdf([1.0, 1.0, 1.0], [1.0, 1.0, 1.0], [1.0, 0.0, -1.0])
+  ! 0.500000000  NaN NaN
+
+  ! `cdf_normal` is pure and, thus, can be called concurrently 
+  xsum = 0.0
+  do concurrent (i=1:size(x,3))
+    xsum = xsum + sum(norm_cdf(x(:,:,i), mu(:,:,i), sigma(:,:,i)))
+  end do
+  print *, xsum
+  ! 11.3751936
+
+  ! complex normal cumulative distribution at x=(0.5,-0.5) with real part of
+  ! mu=1.0, sigma=0.5 and imaginary part of mu=0.0, sigma=1.0
+  loc = (1.0, 0.0)
+  scale = (0.5, 1.0)
+  print *, norm_cdf((0.5, -0.5), loc, scale)
+  ! 4.89511043E-02
+
+end program example_normal_cdf
+

example/stats_distribution_normal/example_normal_cdf.f90

@@ -4,41 +4,56 @@ program example_normal_pdf
                                               norm => rvs_normal
 
   implicit none
-  real :: x(3, 4, 5), a(3, 4, 5), b(3, 4, 5)
+  real, dimension(3, 4, 5) :: x, mu, sigma
+  real :: xsum
   complex :: loc, scale
-  integer :: seed_put, seed_get
+  integer :: seed_put, seed_get, i
 
   seed_put = 1234567
   call random_seed(seed_put, seed_get)
 
-  print *, norm_pdf(1.0, 0., 1.) !a probability density at 1.0 in standard normal
-
-! 0.241970733
+  ! probability density at x=1.0 in standard normal
+  print *, norm_pdf(1.0, 0., 1.) 
+  ! 0.241970733
 
+  ! probability density at x=2.0 with mu=-1.0 and sigma=2.0 
   print *, norm_pdf(2.0, -1.0, 2.0)
-!a probability density at 2.0 with mu=-1.0 sigma=2.0
+  ! 6.47588000E-02
 
-!6.47588000E-02
+  ! probability density at x=1.0 with mu=1.0 and sigma=-1.0 (out of range) 
+  print *, norm_pdf(1.0, 1.0, -1.0)
+  ! NaN
 
+  ! standard normal random variates array
   x = reshape(norm(0.0, 1.0, 60), [3, 4, 5])
-! standard normal random variates array
-
-  a(:, :, :) = 0.0
-  b(:, :, :) = 1.0
-  print *, norm_pdf(x, a, b)  ! standard normal probability density array
-
-!  0.340346158  0.285823315  0.398714304  0.391778737  0.389345556
-!  0.364551932  0.386712372  0.274370432  0.215250477  0.378006011
-!  0.215760440  0.177990928  0.278640658  0.223813817  0.356875211
-!  0.285167664  0.378533930  0.390739858  0.271684974  0.138273031
-!  0.135456234  0.331718773  0.398283750  0.383706540
-
+  
+  ! standard normal probability density array
+  mu(:, :, :) = 0.0
+  sigma(:, :, :) = 1.0
+  print *, norm_pdf(x, mu, sigma)  
+  !  0.340346158  0.285823315  0.398714304  0.391778737  0.389345556
+  !  0.364551932  0.386712372  0.274370432  0.215250477  0.378006011
+  !  0.215760440  0.177990928  0.278640658  0.223813817  0.356875211
+  !  0.285167664  0.378533930  0.390739858  0.271684974  0.138273031
+  !  0.135456234  0.331718773  0.398283750  0.383706540
+
+  ! probability density array where sigma<=0.0 for certain elements 
+  print *, norm_pdf([1.0, 1.0, 1.0], [1.0, 1.0, 1.0], [1.0, 0.0, -1.0])
+  ! 0.398942292  NaN NaN
+
+  ! `pdf_normal` is pure and, thus, can be called concurrently 
+  xsum = 0.0
+  do concurrent (i=1:size(x,3))
+    xsum = xsum + sum(norm_pdf(x(:,:,i), mu(:,:,i), sigma(:,:,i)))
+  end do
+  print *, xsum
+  ! 18.0754433
+
+  ! complex normal probability density function at x=(1.5,1.0) with real part
+  ! of mu=1.0, sigma=1.0 and imaginary part of mu=-0.5, sigma=2.0
   loc = (1.0, -0.5)
   scale = (1.0, 2.)
   print *, norm_pdf((1.5, 1.0), loc, scale)
-! a complex normal probability density function at (1.5,1.0) with real part
-  ! of mu=1.0, sigma=1.0 and imaginary part of mu=-0.5, sigma=2.0
-
-! 5.30100204E-02
+  ! 5.30100204E-02
 
 end program example_normal_pdf