diff --git a/doc/specs/index.md b/doc/specs/index.md index 389fcde0d..7a7a6f143 100644 --- a/doc/specs/index.md +++ b/doc/specs/index.md @@ -25,6 +25,7 @@ This is and index/directory of the specifications (specs) for each new module/fe - [sorting](./stdlib_sorting.html) - Sorting of rank one arrays - [stats](./stdlib_stats.html) - Descriptive Statistics - [stats_distributions_uniform](./stdlib_stats_distribution_uniform.html) - Uniform Probability Distribution + - [stats_distributions_normal](./stdlib_stats_distribution_normal.html) - Normal Probability Distribution - [string\_type](./stdlib_string_type.html) - Basic string support - [strings](./stdlib_strings.html) - String handling and manipulation routines - [stringlist_type](./stdlib_stringlist_type.html) - 1-Dimensional list of strings diff --git a/doc/specs/stdlib_stats_distribution_normal.md b/doc/specs/stdlib_stats_distribution_normal.md new file mode 100644 index 000000000..860aec2dd --- /dev/null +++ b/doc/specs/stdlib_stats_distribution_normal.md @@ -0,0 +1,272 @@ +--- +title: stats_distribution_normal +--- + +# Statistical Distributions -- Normal Distribution Module + +[TOC] + +## `rvs_normal` - normal distribution random variates + +### Status + +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. + +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 three arguments, the function returns a rank one array of normal distributed random variates. + +Note: the algorithm used for generating normal random variates is fundamentally limited to double precision. + +### Syntax + +`result = [[stdlib_stats_distribution_normal(module):rvs_normal(interface)]]([loc, scale] [[, array_size]])` + +### Class + +Function + +### Arguments + +`array_size`: optional argument has `intent(in)` and is a scalar of type `integer`. + +`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`. + +`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`. + +### Example + +```fortran +program demo_normal_rvs + use stdlib_random, only: random_seed + use stdlib_stats_distribution_normal, only: norm => rvs_normal + + implicit none + real :: a(2,3,4), b(2,3,4) + complx :: loc, scale + integer :: seed_put, seed_get + + seed_put = 1234567 + call random_seed(seed_put, seed_get) + + print *, norm( ) !single standard normal random variate + +! 0.563655198 + + print *, norm(1.0, 2.0) + !normal random variate mu=1.0, sigma=2.0 + +! -0.633261681 + + print *, norm(0.0, 1.0, 10) !an array of 10 standard norml random variates + +! -3.38123664E-02 -0.190365672 0.220678389 -0.424612164 -0.249541596 +! 0.865260184 1.11086845 -0.328349441 1.10873628 1.27049923 + + a(:,:,:) = 1.0 + b(:,:,:) = 1.0 + print *, norm(a,b) ! a rank 3 random variates array + +!0.152776539 -7.51764774E-02 1.47208166 0.180561781 1.32407105 +! 1.20383692 0.123445868 -0.455737948 -0.469808221 1.60750175 +! 1.05748117 0.720934749 0.407810807 1.48165631 2.31749439 +! 0.414566994 3.06084275 1.86505437 1.36338580 7.26878643E-02 +! 0.178585172 1.39557445 0.828021586 0.872084975 + + loc = (-1.0, 2.0) + scale = (2.0, 1.0) + print *, norm(loc, scale) + !single complex normal random variate with real part of mu=-1, sigma=2; + !imagainary part of mu=2.0 and sigma=1.0 + +! (1.22566295,2.12518454) + +end program demo_normal_rvs +``` + +## `pdf_normal` - normal distribution probability density function + +### Status + +Experimental + +### Description + +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}}$$ + +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): + +$$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}]}$$ + +### Syntax + +`result = [[stdlib_stats_distribution_normal(module):pdf_normal(interface)]](x, loc, scale)` + +### Class + +Elemental function + +### Arguments + +`x`: has `intent(in)` and is a scalar of type `real` or `complex`. + +`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`. + +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. + +### Example + +```fortran +program demo_normal_pdf + use stdlib_random, only : random_seed + use stdlib_stats_distribution_normal, only : norm_pdf=>pdf_normal, & + norm => rvs_normal + + implicit none + real :: x(3,4,5),a(3,4,5),b(3,4,5) + complx :: loc, scale + integer :: seed_put, seed_get + + 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 + + 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 + + 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 + + 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 + +end program demo_normal_pdf +``` + +## `cdf_normal` - normal distribution cumulative distribution function + +### Status + +Experimental + +### Description + +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 ]$$ + +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): + +$$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}})]$$ + +### Syntax + +`result = [[stdlib_stats_distribution_normal(module):cdf_normal(interface)]](x, loc, scale)` + +### Class + +Elemental function + +### Arguments + +`x`: has `intent(in)` and is a scalar of type `real` or `complex`. + +`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`. + +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. + +### Example + +```fortran +program demo_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) + complx :: 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 demo_norm_cdf + +``` diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 44c304c02..6a2f242cc 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -27,6 +27,7 @@ set(fppFiles stdlib_stats_moment_mask.fypp stdlib_stats_moment_scalar.fypp stdlib_stats_distribution_uniform.fypp + stdlib_stats_distribution_normal.fypp stdlib_stats_var.fypp stdlib_quadrature.fypp stdlib_quadrature_trapz.fypp diff --git a/src/Makefile.manual b/src/Makefile.manual index 2bfcb7091..7b5f75938 100644 --- a/src/Makefile.manual +++ b/src/Makefile.manual @@ -29,6 +29,7 @@ SRCFYPP = \ stdlib_stats_moment_mask.fypp \ stdlib_stats_moment_scalar.fypp \ stdlib_stats_distribution_uniform.fypp \ + stdlib_stats_distribution_normal.fypp \ stdlib_stats_var.fypp \ stdlib_math.fypp \ stdlib_math_linspace.fypp \ @@ -154,6 +155,11 @@ stdlib_stats_distribution_uniform.o: \ stdlib_kinds.o \ stdlib_error.o \ stdlib_random.o +stdlib_stats_distribution_normal.o: \ + stdlib_kinds.o \ + stdlib_error.o \ + stdlib_random.o \ + stdlib_stats_distribution_uniform.o stdlib_random.o: \ stdlib_kinds.o \ stdlib_error.o diff --git a/src/stdlib_stats_distribution_normal.fypp b/src/stdlib_stats_distribution_normal.fypp new file mode 100644 index 000000000..e21b7af14 --- /dev/null +++ b/src/stdlib_stats_distribution_normal.fypp @@ -0,0 +1,332 @@ +#:include "common.fypp" +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES +module stdlib_stats_distribution_normal + use stdlib_kinds, only : sp, dp, xdp, qp, int32 + use stdlib_error, only : error_stop + use stdlib_random, only : dist_rand + use stdlib_stats_distribution_uniform, only : uni=>rvs_uniform + + implicit none + private + + real(dp), parameter :: HALF = 0.5_dp, ONE = 1.0_dp, TWO = 2.0_dp + integer :: kn(0:127) + real(dp) :: wn(0:127), fn(0:127) + logical :: zig_norm_initialized = .false. + + public :: rvs_normal + public :: pdf_normal + public :: cdf_normal + + + + interface rvs_normal + !! version: experimental + !! + !! Normal Distribution Random Variates + !! ([Specification](../page/specs/stdlib_stats_distribution_normal.html# + !! rvs_normal-normal-distribution-random-variates)) + !! + module procedure rvs_norm_0_rsp !0 dummy variable + + #:for k1, t1 in RC_KINDS_TYPES + module procedure rvs_norm_${t1[0]}$${k1}$ !2 dummy variables + #:endfor + + #:for k1, t1 in RC_KINDS_TYPES + module procedure rvs_norm_array_${t1[0]}$${k1}$ !3 dummy variables + #:endfor + end interface rvs_normal + + + + interface pdf_normal + !! version: experimental + !! + !! Normal Distribution Probability Density Function + !! ([Specification](../page/specs/stdlib_stats_distribution_normal.html# + !! pdf_normal-normal-distribution-probability-density-function)) + !! + #:for k1, t1 in RC_KINDS_TYPES + module procedure pdf_norm_${t1[0]}$${k1}$ + #:endfor + end interface pdf_normal + + + + interface cdf_normal + !! version: experimental + !! + !! Normal Distribution Cumulative Distribution Function + !! ([Specification](../page/specs/stdlib_stats_distribution_normal.html# + !! cdf_normal-normal-distribution-cumulative-distribution-function)) + !! + #:for k1, t1 in RC_KINDS_TYPES + module procedure cdf_norm_${t1[0]}$${k1}$ + #:endfor + end interface cdf_normal + + + + + +contains + + subroutine zigset + ! Marsaglia & Tsang generator for random normals & random exponentials. + ! Translated from C by Alan Miller (amiller@bigpond.net.au), released as public + ! domain (https://jblevins.org/mirror/amiller/) + ! + ! Marsaglia, G. & Tsang, W.W. (2000) `The ziggurat method for generating + ! random variables', J. Statist. Software, v5(8). + ! + ! This is an electronic journal which can be downloaded from: + ! http://www.jstatsoft.org/v05/i08 + ! + ! Latest version - 1 January 2001 + ! + real(dp), parameter :: M1 = 2147483648.0_dp, vn = 0.00991256303526217_dp + real(dp) :: dn, tn, q + integer :: i + + dn = 3.442619855899_dp + tn = dn + !tables for random normals + q = vn * exp(HALF * dn * dn) + kn(0) = int((dn / q) * M1, kind = int32) + kn(1) = 0 + wn(0) = q / M1 + wn(127) = dn / M1 + fn(0) = ONE + fn(127) = exp( -HALF * dn * dn ) + do i = 126, 1, -1 + dn = sqrt( -TWO * log( vn / dn + exp( -HALF * dn * dn ) ) ) + kn(i+1) = int((dn / tn) * M1, kind = int32) + tn = dn + fn(i) = exp(-HALF * dn * dn) + wn(i) = dn / M1 + end do + zig_norm_initialized = .true. + end subroutine zigset + + + + #:for k1, t1 in REAL_KINDS_TYPES + function rvs_norm_0_${t1[0]}$${k1}$( ) result(res) + ! + ! Standard normal random vairate (0,1) + ! + ${t1}$ :: res + ${t1}$, parameter :: r = 3.442619855899_${k1}$, rr = 1.0_${k1}$ / r + ${t1}$ :: x, y + integer :: hz, iz + + if(.not. zig_norm_initialized) call zigset + iz = 0 + hz = dist_rand(1_int32) !32bit random integer + iz = iand( hz, 127 ) !random integer in [0, 127] + if( abs( hz ) < kn(iz) ) then + res = hz * wn(iz) + else + L1: do + L2: if( iz == 0 ) then + do + x = -log( uni(1.0_${k1}$) ) * rr + y = -log( uni(1.0_${k1}$) ) + if( y + y >= x * x ) exit + end do + res = r + x + if( hz <= 0 ) res = -res + exit L1 + end if L2 + x = hz * wn(iz) + if( fn(iz) + uni(1.0_${k1}$) * (fn(iz-1) - fn(iz)) < & + exp(-HALF * x * x) ) then + res = x + exit L1 + end if + hz = dist_rand(1_int32) + iz = iand( hz, 127 ) + if( abs( hz ) < kn(iz) ) then + res = hz * wn(iz) + exit L1 + end if + end do L1 + end if + end function rvs_norm_0_${t1[0]}$${k1}$ + + #:endfor + + + + #:for k1, t1 in REAL_KINDS_TYPES + function rvs_norm_${t1[0]}$${k1}$(loc, scale) result(res) + ! + ! Normal random variate (loc, scale) + ! + ${t1}$, intent(in) :: loc, scale + ${t1}$ :: res + + if(scale == 0._${k1}$) call error_stop("Error(rvs_norm): Normal" & + //" distribution scale parameter must be non-zero") + res = rvs_norm_0_${t1[0]}$${k1}$( ) + res = res * scale + loc + end function rvs_norm_${t1[0]}$${k1}$ + + #:endfor + + + + #:for k1, t1 in CMPLX_KINDS_TYPES + function rvs_norm_${t1[0]}$${k1}$(loc, scale) result(res) + ! + ! Normally distributed complex. The real part and imaginary part are & + ! independent of each other. + ! + ${t1}$, intent(in) :: loc, scale + ${t1}$ :: res + real(${k1}$) :: tr, ti + + tr = rvs_norm_r${k1}$(loc % re, scale % re) + ti = rvs_norm_r${k1}$(loc % im, scale % im) + res = cmplx(tr, ti, kind=${k1}$) + end function rvs_norm_${t1[0]}$${k1}$ + + #:endfor + + + + #:for k1, t1 in REAL_KINDS_TYPES + function rvs_norm_array_${t1[0]}$${k1}$(loc, scale, array_size) result(res) + ${t1}$, intent(in) :: loc, scale + integer, intent(in) :: array_size + ${t1}$ :: res(array_size) + ${t1}$, parameter :: r = 3.442619855899_${k1}$, rr = 1.0_${k1}$ / r + ${t1}$ :: x, y, re + integer :: hz, iz, i + + if(scale == 0._${k1}$) call error_stop("Error(rvs_norm_array): Normal" & + //"distribution scale parameter must be non-zero") + if(.not. zig_norm_initialized) call zigset + do i = 1, array_size + iz = 0 + hz = dist_rand(1_int32) + iz = iand( hz, 127 ) + if( abs( hz ) < kn(iz) ) then + re = hz * wn(iz) + else + L1: do + L2: if( iz == 0 ) then + do + x = -log( uni(1.0_${k1}$) ) * rr + y = -log( uni(1.0_${k1}$) ) + if( y + y >= x * x ) exit + end do + re = r + x + if( hz <= 0 ) re = -re + exit L1 + end if L2 + x = hz * wn(iz) + if( fn(iz) + uni(1.0_${k1}$) * (fn(iz-1) - fn(iz)) < & + exp(-HALF * x * x) ) then + re = x + exit L1 + end if + + hz = dist_rand(1_int32) + iz = iand( hz, 127 ) + if( abs( hz ) < kn(iz) ) then + re = hz * wn(iz) + exit L1 + end if + end do L1 + end if + res(i) = re * scale + loc + end do + end function rvs_norm_array_${t1[0]}$${k1}$ + + #:endfor + + + + #:for k1, t1 in CMPLX_KINDS_TYPES + function rvs_norm_array_${t1[0]}$${k1}$(loc, scale, array_size) result(res) + ${t1}$, intent(in) :: loc, scale + integer, intent(in) :: array_size + integer :: i + ${t1}$ :: res(array_size) + real(${k1}$) :: tr, ti + + do i = 1, array_size + tr = rvs_norm_r${k1}$(loc % re, scale % re) + ti = rvs_norm_r${k1}$(loc % im, scale % im) + res(i) = cmplx(tr, ti, kind=${k1}$) + end do + end function rvs_norm_array_${t1[0]}$${k1}$ + + #:endfor + + + + #:for k1, t1 in REAL_KINDS_TYPES + impure elemental function pdf_norm_${t1[0]}$${k1}$(x, loc, scale) result(res) + ! + ! Normal distribution probability density function + ! + ${t1}$, intent(in) :: x, loc, scale + ${t1}$ :: res + ${t1}$, parameter :: sqrt_2_pi = sqrt(2.0_${k1}$ * acos(-1.0_${k1}$)) + + if(scale == 0._${k1}$) call error_stop("Error(pdf_norm): Normal" & + //"distribution scale parameter must be non-zero") + res = exp(- 0.5_${k1}$ * ((x - loc) / scale) * (x - loc) / scale) / & + (sqrt_2_Pi * scale) + end function pdf_norm_${t1[0]}$${k1}$ + + #:endfor + + + + #:for k1, t1 in CMPLX_KINDS_TYPES + impure elemental function pdf_norm_${t1[0]}$${k1}$(x, loc, scale) result(res) + ${t1}$, intent(in) :: x, loc, scale + real(${k1}$) :: res + + res = pdf_norm_r${k1}$(x % re, loc % re, scale % re) + res = res * pdf_norm_r${k1}$(x % im, loc % im, scale % im) + end function pdf_norm_${t1[0]}$${k1}$ + + #:endfor + + + + #:for k1, t1 in REAL_KINDS_TYPES + impure elemental function cdf_norm_${t1[0]}$${k1}$(x, loc, scale) result(res) + ! + ! Normal distribution cumulative distribution function + ! + ${t1}$, intent(in) :: x, loc, scale + ${t1}$ :: res + ${t1}$, parameter :: sqrt_2 = sqrt(2.0_${k1}$) + + if(scale == 0._${k1}$) call error_stop("Error(cdf_norm): Normal" & + //"distribution scale parameter must be non-zero") + res = erfc(- (x - loc) / (scale * sqrt_2)) / 2.0_${k1}$ + end function cdf_norm_${t1[0]}$${k1}$ + + #:endfor + + + + #:for k1, t1 in CMPLX_KINDS_TYPES + impure elemental function cdf_norm_${t1[0]}$${k1}$(x, loc, scale) result(res) + ${t1}$, intent(in) :: x, loc, scale + real(${k1}$) :: res + + res = cdf_norm_r${k1}$(x % re, loc % re, scale % re) + res = res * cdf_norm_r${k1}$(x % im, loc % im, scale % im) + end function cdf_norm_${t1[0]}$${k1}$ + + #:endfor + +end module stdlib_stats_distribution_normal diff --git a/src/tests/stats/CMakeLists.txt b/src/tests/stats/CMakeLists.txt index 3df736de5..b3b70bd4d 100644 --- a/src/tests/stats/CMakeLists.txt +++ b/src/tests/stats/CMakeLists.txt @@ -6,6 +6,7 @@ set(fppFiles test_mean_f03.fypp test_median.fypp test_distribution_uniform.fypp + test_distribution_normal.fypp ) fypp_f90("${fyppFlags}" "${fppFiles}" outFiles) @@ -20,6 +21,7 @@ ADDTEST(var) ADDTEST(varn) ADDTEST(random) ADDTEST(distribution_uniform) +ADDTEST(distribution_normal) if(DEFINED CMAKE_MAXIMUM_RANK) if(${CMAKE_MAXIMUM_RANK} GREATER 7) diff --git a/src/tests/stats/Makefile.manual b/src/tests/stats/Makefile.manual index 181ba61c7..15797d509 100644 --- a/src/tests/stats/Makefile.manual +++ b/src/tests/stats/Makefile.manual @@ -1,7 +1,8 @@ SRCFYPP = \ test_mean.fypp \ test_median.fypp \ - test_distribution_uniform.fypp + test_distribution_uniform.fypp \ + test_distribution_normal.fypp SRCGEN = $(SRCFYPP:.fypp=.f90) @@ -10,4 +11,5 @@ $(SRCGEN): %.f90: %.fypp ../../common.fypp PROGS_SRC = $(SRCGEN) test_mean.f90 test_moment.f90 test_var.f90 test_random.f90 + include ../Makefile.manual.test.mk diff --git a/src/tests/stats/test_distribution_normal.fypp b/src/tests/stats/test_distribution_normal.fypp new file mode 100644 index 000000000..72e2746a4 --- /dev/null +++ b/src/tests/stats/test_distribution_normal.fypp @@ -0,0 +1,278 @@ + +#:include "common.fypp" +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES +program test_distribution_normal + use stdlib_kinds, only : sp, dp, xdp, qp + use stdlib_error, only : check + use stdlib_random, only : random_seed + use stdlib_stats_distribution_uniform, only : uni => rvs_uniform + use stdlib_stats_distribution_normal, only : nor_rvs => rvs_normal, & + nor_pdf => pdf_normal, nor_cdf => cdf_normal + + implicit none + #:for k1, t1 in REAL_KINDS_TYPES + ${t1}$, parameter :: ${k1}$tol = 1000 * epsilon(1.0_${k1}$) + #:endfor + logical :: warn = .true. + integer :: put, get + + put = 12345678 + call random_seed(put, get) + + call test_normal_random_generator + + + #:for k1, t1 in RC_KINDS_TYPES + call test_nor_rvs_${t1[0]}$${k1}$ + #:endfor + + + + #:for k1, t1 in RC_KINDS_TYPES + call test_nor_pdf_${t1[0]}$${k1}$ + #:endfor + + + + #:for k1, t1 in RC_KINDS_TYPES + call test_nor_cdf_${t1[0]}$${k1}$ + #:endfor + + + + + +contains + + subroutine test_normal_random_generator + + integer, parameter :: num = 10000000, array_size = 1000 + integer :: i, j, freq(0:array_size) + real(dp) :: chisq, expct + + print *, "" + print *, "Test normal random generator with chi-squared" + freq = 0 + do i = 1, num + j = array_size * (1 + erf(nor_rvs(0.0, 1.0) / sqrt(2.0))) / 2.0 + freq(j) = freq(j) + 1 + end do + chisq = 0.0_dp + expct = num / array_size + do i = 0, array_size - 1 + chisq = chisq + (freq(i) - expct) ** 2 / expct + end do + write(*,*) "The critical values for chi-squared with 1000 d. of f. is" & + //" 1143.92" + write(*,*) "Chi-squared for normal random generator is : ", chisq + call check((chisq < 1143.9), & + msg = "normal randomness failed chi-squared test", warn = warn) + end subroutine test_normal_random_generator + + + + + #:for k1, t1 in RC_KINDS_TYPES + subroutine test_nor_rvs_${t1[0]}$${k1}$ + ${t1}$ :: res(10), loc, scale + integer, parameter :: k = 5 + integer :: i, n + integer :: seed, get + #:if t1[0] == "r" + #! for real type + ${t1}$, parameter :: ans(10) = & + [2.66708039318040679432897377409972250_${k1}$, & + 2.36030794936128329730706809641560540_${k1}$, & + 1.27712218793084242296487218482070602_${k1}$, & + -2.39132544130814794769435138732660562_${k1}$, & + 1.72566595106028652928387145948363468_${k1}$, & + -1.50621775537767632613395107910037041_${k1}$, & + 2.13518835158352082714827702147886157_${k1}$, & + -0.636788253742142318358787633769679815_${k1}$, & + 2.48600787778845799813609573902795091_${k1}$, & + -3.03711473837981227319460231228731573_${k1}$] + #:else + #! for complex type + ${t1}$, parameter :: ans(10) = & + [(2.12531029488530509574673033057479188_${k1}$, & + 1.46507698734032082432676702410390135_${k1}$), & + (1.08284164094813181722365413861552952_${k1}$, & + 0.277168639672963013076412153168348595_${k1}$), & + (1.41924946329521489696290359461272601_${k1}$, & + 0.498445561155580918466512230224907398_${k1}$), & + (1.72639126368764062036120776610914618_${k1}$, & + 0.715802936564464420410303091557580046_${k1}$), & + (1.98950590834134349860207180427096318_${k1}$, & + 0.115721315405046931701349421928171068_${k1}$), & + (-1.16929014824793620075382705181255005_${k1}$, & + 0.250744737486995217246033007540972903_${k1}$), & + (1.57160542831869509683428987045772374_${k1}$, & + 0.638282596371312238581197107123443857_${k1}$), & + (-1.36106107654239116833139178197598085_${k1}$, & + 0.166259201494369124318950525776017457_${k1}$), & + (1.13403096805387920698038328737311531_${k1}$, & + 1.04232618148691447146347854868508875_${k1}$), & + (-1.68220535920475811053620418533682823_${k1}$, & + 1.63361446685040256898702182297711261_${k1}$)] + #:endif + + print *, "Test normal_distribution_rvs_${t1[0]}$${k1}$" + seed = 25836914 + call random_seed(seed, get) + #:if t1[0] == "r" + #! for real type + loc = 0.5_${k1}$; scale = 2.0_${k1}$ + #:else + #! for complex type + loc = (0.5_${k1}$, 1.0_${k1}$); scale = (1.5_${k1}$, 0.5_${k1}$) + #:endif + + do i = 1, k + res(i) = nor_rvs(loc, scale) ! 2 dummies + end do + res(6:10) = nor_rvs(loc, scale, k) ! 3 dummies + call check(all(abs(res - ans) < ${k1}$tol), & + msg="normal_distribution_rvs_${t1[0]}$${k1}$ failed", warn=warn) + end subroutine test_nor_rvs_${t1[0]}$${k1}$ + + #:endfor + + + + + #:for k1, t1 in RC_KINDS_TYPES + subroutine test_nor_pdf_${t1[0]}$${k1}$ + + ${t1}$ :: x1, x2(3,4), loc, scale + integer, parameter :: k = 5 + integer :: i, n, seed, get + real(${k1}$) :: res(3,5) + #:if t1[0] == "r" + #! for real type + real(${k1}$), parameter :: ans(15) = & + [0.215050766989949083210785218076278553_${k1}$, & + 0.215050766989949083210785218076278553_${k1}$, & + 0.215050766989949083210785218076278553_${k1}$, & + 0.200537626692880596839299818454439236_${k1}$, & + 5.66161527403268434368575024104022496E-0002_${k1}$, & + 0.238986957612021514867582359518579138_${k1}$, & + 0.265935969411942911029638292783132425_${k1}$, & + 0.262147558654079961109031890374902943_${k1}$, & + 0.249866408914952245533320687656894701_${k1}$, & + 3.98721117498705317877792757313696510E-0002_${k1}$, & + 0.265902369803533466897906694845094995_${k1}$, & + 0.161311603170650092038944290133124635_${k1}$, & + 0.249177740354276111998717092437037695_${k1}$, & + 0.237427217242213206474603807278971527_${k1}$, & + 0.155696086384122017518186260628090478_${k1}$] + #:else + #! for complex type + real(${k1}$), parameter :: ans(15) = & + [0.129377311291944176372137325120411497_${k1}$, & + 0.129377311291944176372137325120411497_${k1}$, & + 0.129377311291944176372137325120411497_${k1}$, & + 4.05915662853246811934977653001971736E-0002_${k1}$, & + 0.209143395418940756076861773161637924_${k1}$, & + 2.98881041363874672676853084975547667E-0002_${k1}$, & + 0.128679412679649127469385460133445161_${k1}$, & + 0.177484732473055532384223611956177231_${k1}$, & + 3.82205306942578982084957100753849738E-0002_${k1}$, & + 7.09915714309796034515515428785324918E-0002_${k1}$, & + 4.56126582912124629544443072483362352E-0002_${k1}$, & + 6.57454133967021123696499056531595921E-0002_${k1}$, & + 0.165161039915667041643464172210282279_${k1}$, & + 3.86104822953520989775015755966798359E-0002_${k1}$, & + 0.196922947431391188040943672442575686_${k1}$] + #:endif + + print *, "Test normal_distribution_pdf_${t1[0]}$${k1}$" + seed = 741852963 + call random_seed(seed, get) + #:if t1[0] == "r" + #! for real type + loc = -0.5_${k1}$; scale = 1.5_${k1}$ + #:else + #! for complex type + loc = (-0.5_${k1}$, 0.5_${k1}$); scale = (0.5_${k1}$, 1.5_${k1}$) + #:endif + + x1 = nor_rvs(loc, scale) + x2 = reshape(nor_rvs(loc, scale, 12), [3,4]) + res(:,1) = nor_pdf(x1, loc, scale) + res(:, 2:5) = nor_pdf(x2, loc, scale) + call check(all(abs(res - reshape(ans, [3,5])) < ${k1}$tol), & + msg="normal_distribution_pdf_${t1[0]}$${k1}$ failed", warn=warn) + end subroutine test_nor_pdf_${t1[0]}$${k1}$ + + #:endfor + + + + + #:for k1, t1 in RC_KINDS_TYPES + subroutine test_nor_cdf_${t1[0]}$${k1}$ + + ${t1}$ :: x1, x2(3,4), loc, scale + integer :: i, n + integer :: seed, get + real(${k1}$) :: res(3,5) + #:if t1[0] == "r" + #!for real type + real(${k1}$), parameter :: ans(15) = & + [7.50826305038441048487991102776953948E-0002_${k1}$, & + 7.50826305038441048487991102776953948E-0002_${k1}$, & + 7.50826305038441048487991102776953948E-0002_${k1}$, & + 0.143119834108717983250834016885129863_${k1}$, & + 0.241425421525703182028420560765471735_${k1}$, & + 0.284345878626039240974266199229875972_${k1}$, & + 0.233239836366015928845367994433532757_${k1}$, & + 0.341059506137219171082517155967522896_${k1}$, & + 0.353156850199835111081038166086606192_${k1}$, & + 6.81066766396638231790017005897813244E-0002_${k1}$, & + 4.38792331441682923984716366123285346E-0002_${k1}$, & + 0.763679637882860826030745070304416929_${k1}$, & + 0.363722187587355040667876190724308059_${k1}$, & + 0.868187114884980488672309198087692444_${k1}$, & + 0.626506799809652872401992867475200722_${k1}$] + #:else + #! for complex type + real(${k1}$), parameter :: ans(15) = & + [1.07458136221563368133842063954746170E-0002_${k1}$, & + 1.07458136221563368133842063954746170E-0002_${k1}$, & + 1.07458136221563368133842063954746170E-0002_${k1}$, & + 6.86483236063879585051085536740820057E-0002_${k1}$, & + 7.95486634025192048896990048539218724E-0002_${k1}$, & + 2.40523393996423661445007940057223384E-0002_${k1}$, & + 3.35096768781160662250307446207445131E-0002_${k1}$, & + 0.315778916661119434962814841317323376_${k1}$, & + 0.446311293878359175362094845206410428_${k1}$, & + 0.102010220821382542292905161748120877_${k1}$, & + 7.66919007012121545175655727052974512E-0002_${k1}$, & + 0.564690968410069125818268877247699603_${k1}$, & + 0.708769523556518785240723539383512333_${k1}$, & + 6.40553790808161720088070925562830659E-0002_${k1}$, & + 5.39999153072107729358158443133850711E-0002_${k1}$] + #:endif + + print *, "Test normal_distribution_cdf_${t1[0]}$${k1}$" + seed = 369147582 + call random_seed(seed, get) + #:if t1[0] == "r" + #! for real type + loc = -1.0_${k1}$; scale = 2.0_${k1}$ + #:else + #! for complex type + loc = (-1.0_${k1}$, 1.0_${k1}$); scale = (1.7_${k1}$, 2.4_${k1}$) + #:endif + + x1 = nor_rvs(loc, scale) + x2 = reshape(nor_rvs(loc, scale, 12), [3,4]) + res(:,1) = nor_cdf(x1, loc, scale) + res(:, 2:5) = nor_cdf(x2, loc, scale) + call check(all(abs(res - reshape(ans,[3,5])) < ${k1}$tol), & + msg="normal_distribution_cdf_${t1[0]}$${k1}$ failed", warn=warn) + end subroutine test_nor_cdf_${t1[0]}$${k1}$ + + #:endfor + +end program test_distribution_normal