From 1abb80be687802fb3c77e65455bc9513ef5c9a66 Mon Sep 17 00:00:00 2001 From: HugoMVale <57530119+HugoMVale@users.noreply.github.com> Date: Thu, 8 Jun 2023 12:32:50 +0200 Subject: [PATCH 1/7] improve readability and equation rendering --- doc/specs/stdlib_stats_distribution_normal.md | 26 ++++++++++--------- 1 file changed, 14 insertions(+), 12 deletions(-) diff --git a/doc/specs/stdlib_stats_distribution_normal.md b/doc/specs/stdlib_stats_distribution_normal.md index 20d8b5956..8378050bc 100644 --- a/doc/specs/stdlib_stats_distribution_normal.md +++ b/doc/specs/stdlib_stats_distribution_normal.md @@ -14,13 +14,13 @@ 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. @@ -44,7 +44,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 as the same type of `scale` and `loc`. ### Example @@ -62,11 +62,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 (ref. "Probability and Random Processes with Applications to Signal Processing and Communications", 2nd ed., Scott L. Miller and Donald Childers, 2012, p.197): -$$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 @@ -88,7 +88,7 @@ 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 of the same type as the input arguments. ### Example @@ -106,11 +106,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 (ref. "Probability and Random Processes with Applications to Signal Processing and Communications", 2nd ed., Scott L. Miller and Donald Childers, 2012, p.197): -$$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 From 99ce6e004b421776d7f23c7c3be125774ea9efc4 Mon Sep 17 00:00:00 2001 From: HugoMVale <57530119+HugoMVale@users.noreply.github.com> Date: Thu, 8 Jun 2023 14:19:14 +0200 Subject: [PATCH 2/7] clean + add examples of call with sigma<=0 and do concurrent --- .../example_normal_pdf.f90 | 61 ++++++++++++------- 1 file changed, 38 insertions(+), 23 deletions(-) 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 From 9adb14e34968ee0ef0c9700b9e7883b45c3c4aa0 Mon Sep 17 00:00:00 2001 From: HugoMVale <57530119+HugoMVale@users.noreply.github.com> Date: Thu, 8 Jun 2023 14:51:20 +0200 Subject: [PATCH 3/7] rename file for consistency, clean/hamornize + add examples of call with sigma<=0 and do concurrent --- .../stats_distribution_normal/CMakeLists.txt | 2 +- .../example_norm_cdf.f90 | 45 -------------- .../example_normal_cdf.f90 | 60 +++++++++++++++++++ 3 files changed, 61 insertions(+), 46 deletions(-) delete mode 100644 example/stats_distribution_normal/example_norm_cdf.f90 create mode 100644 example/stats_distribution_normal/example_normal_cdf.f90 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 + From fa4758f174ff4b8998c9fcf25c18e589e12f6a97 Mon Sep 17 00:00:00 2001 From: HugoMVale <57530119+HugoMVale@users.noreply.github.com> Date: Thu, 8 Jun 2023 15:25:19 +0200 Subject: [PATCH 4/7] specify that sigma must be positive, otherwise result is NaN --- doc/specs/stdlib_stats_distribution_normal.md | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/doc/specs/stdlib_stats_distribution_normal.md b/doc/specs/stdlib_stats_distribution_normal.md index 8378050bc..bd9c77bf1 100644 --- a/doc/specs/stdlib_stats_distribution_normal.md +++ b/doc/specs/stdlib_stats_distribution_normal.md @@ -14,7 +14,7 @@ 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 ($\mu$). The `scale` specifies the standard deviation ($\sigma$). +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)$. @@ -82,13 +82,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 the arguments, and of the same type as the input arguments. +The result is a scalar or an array, with a shape conformable to the arguments, and of the same type as the input arguments. If `scale` is non-positive, the result is `NaN`. ### Example @@ -128,13 +128,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, as the same type of input arguments. +The result is a scalar or an array, with a shape conformable to arguments, as the same type of input arguments. If `scale` is non-positive, the result is `NaN`. ### Example From e1a85fb861bbd0ce69f908d06492eea00f37a8cb Mon Sep 17 00:00:00 2001 From: HugoMVale <57530119+HugoMVale@users.noreply.github.com> Date: Thu, 8 Jun 2023 16:04:29 +0200 Subject: [PATCH 5/7] add NaN spec to rvs_normal --- doc/specs/stdlib_stats_distribution_normal.md | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/doc/specs/stdlib_stats_distribution_normal.md b/doc/specs/stdlib_stats_distribution_normal.md index bd9c77bf1..bba8c2a95 100644 --- a/doc/specs/stdlib_stats_distribution_normal.md +++ b/doc/specs/stdlib_stats_distribution_normal.md @@ -36,7 +36,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 +44,7 @@ Elemental function (passing both `loc` and `scale`). ### Return value -The result is a scalar or rank-1 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 @@ -88,7 +88,7 @@ All three arguments must have the same type. ### Return value -The result is a scalar or an array, with a shape conformable to the arguments, and of the same type as the input arguments. If `scale` is non-positive, the result is `NaN`. +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 @@ -134,7 +134,7 @@ 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. If `scale` is non-positive, the result is `NaN`. +The result is a scalar or an array, with a shape conformable to arguments, as the same type as the input arguments. If `scale` is non-positive, the result is `NaN`. ### Example From 268ca2bdd6b44d4166a84e95caf6a715d99fdb68 Mon Sep 17 00:00:00 2001 From: HugoMVale <57530119+HugoMVale@users.noreply.github.com> Date: Thu, 8 Jun 2023 16:25:33 +0200 Subject: [PATCH 6/7] fix typo --- doc/specs/stdlib_stats_distribution_normal.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_stats_distribution_normal.md b/doc/specs/stdlib_stats_distribution_normal.md index bba8c2a95..3dfbb41cb 100644 --- a/doc/specs/stdlib_stats_distribution_normal.md +++ b/doc/specs/stdlib_stats_distribution_normal.md @@ -134,7 +134,7 @@ 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 as the input arguments. If `scale` is non-positive, the result is `NaN`. +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 From 32dfb474b16e05a50992c48108a02d68570fbe53 Mon Sep 17 00:00:00 2001 From: HugoMVale <57530119+HugoMVale@users.noreply.github.com> Date: Thu, 8 Jun 2023 22:43:14 +0200 Subject: [PATCH 7/7] move refs to footnote --- doc/specs/stdlib_stats_distribution_normal.md | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/doc/specs/stdlib_stats_distribution_normal.md b/doc/specs/stdlib_stats_distribution_normal.md index bd9c77bf1..6607a020f 100644 --- a/doc/specs/stdlib_stats_distribution_normal.md +++ b/doc/specs/stdlib_stats_distribution_normal.md @@ -22,7 +22,8 @@ With two arguments, the function returns a normal distributed random variate $N( 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 @@ -64,7 +65,7 @@ The probability density function (pdf) of the single real variable normal distri $$f(x) = \frac{1}{\sigma \sqrt{2}} \exp{\left[-\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^{2}\right]}$$ -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 (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}} \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]}$$ @@ -108,7 +109,7 @@ Cumulative distribution function of the single real variable normal distribution $$F(x) = \frac{1}{2}\left [ 1+\text{erf}\left(\frac{x-\mu}{\sigma \sqrt{2}}\right) \right ]$$ -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 (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} \ \left[ 1+\text{erf}\left(\frac{x-\mu_x}{\sigma_x \sqrt{2}}\right) \right] \ @@ -141,3 +142,7 @@ The result is a scalar or an array, with a shape conformable to arguments, as th ```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