diff --git a/doc/specs/stdlib_stats_distribution_normal.md b/doc/specs/stdlib_stats_distribution_normal.md index 20d8b5956..7dfddaefe 100644 --- a/doc/specs/stdlib_stats_distribution_normal.md +++ b/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,7 +37,7 @@ 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`. @@ -44,7 +45,7 @@ Elemental function (passing both `loc` and `scale`). ### 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). \ No newline at end of file diff --git a/example/stats_distribution_normal/CMakeLists.txt b/example/stats_distribution_normal/CMakeLists.txt index 955f00ee1..c4c57cad9 100644 --- a/example/stats_distribution_normal/CMakeLists.txt +++ b/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) diff --git a/example/stats_distribution_normal/example_norm_cdf.f90 b/example/stats_distribution_normal/example_norm_cdf.f90 deleted file mode 100644 index a67450277..000000000 --- a/example/stats_distribution_normal/example_norm_cdf.f90 +++ /dev/null @@ -1,45 +0,0 @@ -program example_norm_cdf - use stdlib_random, only: random_seed - use stdlib_stats_distribution_normal, only: norm_cdf => cdf_normal, & - norm => rvs_normal - - implicit none - real :: x(2, 3, 4), a(2, 3, 4), b(2, 3, 4) - complex :: loc, scale - integer :: seed_put, seed_get - - seed_put = 1234567 - call random_seed(seed_put, seed_get) - - print *, norm_cdf(1.0, 0.0, 1.0) ! a standard normal cumulative at 1.0 - -! 0.841344714 - - print *, norm_cdf(2.0, -1.0, 2.0) -! a cumulative at 2.0 with mu=-1 sigma=2 - -! 0.933192849 - - x = reshape(norm(0.0, 1.0, 24), [2, 3, 4]) -! standard normal random variates array - - a(:, :, :) = 0.0 - b(:, :, :) = 1.0 - print *, norm_cdf(x, a, b) ! standard normal cumulative array - -! 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 - - loc = (1.0, 0.0) - scale = (0.5, 1.0) - print *, norm_cdf((0.5, -0.5), loc, scale) -!complex normal cumulative distribution at (0.5,-0.5) with real part of - !mu=1.0, sigma=0.5 and imaginary part of mu=0.0, sigma=1.0 - -!4.89511043E-02 - -end program example_norm_cdf - diff --git a/example/stats_distribution_normal/example_normal_cdf.f90 b/example/stats_distribution_normal/example_normal_cdf.f90 new file mode 100644 index 000000000..6ac4f5092 --- /dev/null +++ b/example/stats_distribution_normal/example_normal_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 + diff --git a/example/stats_distribution_normal/example_normal_pdf.f90 b/example/stats_distribution_normal/example_normal_pdf.f90 index 2dc62a370..82bbe68e1 100644 --- a/example/stats_distribution_normal/example_normal_pdf.f90 +++ b/example/stats_distribution_normal/example_normal_pdf.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