diff --git a/doc/specs/index.md b/doc/specs/index.md index 91284c2df..c52902473 100644 --- a/doc/specs/index.md +++ b/doc/specs/index.md @@ -17,6 +17,7 @@ This is and index/directory of the specifications (specs) for each new module/fe - [optval](./stdlib_optval.html) - Fallback value for optional arguments - [quadrature](./stdlib_quadrature.html) - Numerical integration - [stats](./stdlib_stats.html) - Descriptive Statistics + - [distributions](./stdlib_stats_distribution.html) - Random Distributions and Statistical Functions ## Missing specs diff --git a/doc/specs/stdlib_stats_distribution.md b/doc/specs/stdlib_stats_distribution.md new file mode 100644 index 000000000..703f3dc57 --- /dev/null +++ b/doc/specs/stdlib_stats_distribution.md @@ -0,0 +1,456 @@ +--- +title: stats_distribution +--- + +# Statistical Distributions + +[TOC] + +## `random_seed` - set or get a value of seed to the random distribution pseudorandom number generator + +### Status + +Experimental + +### Description + +Sets the seed value before calling random distribution for variates. + +### Syntax + +`call random_seed(put, get)` + +### Arguments + +`put`: argument has intent `in` and may be a scalar of type `integer` with kind of `int32`. + +`get`: argument has intent `out` and is a scalar of type `integer` with kind of `int32`. + +### Return value + +Return a scalar of type `integer` with kind of `int32` + +### Example + +```fortran +program demo_norm_seed + use stdlib_stats_distribution, only : random_seed + use iso_fortran_env, only : int32 + implicit none + integer(int32) :: seed, seed_value + + seed = 1234567 + call random_seed(seed, seed_value) ! set and get current value of seed +end program demo_norm_seed +``` + +## `uniform_distribution_rvs` - uniform distribution random variates + +### Status + +Experimental + +### Description + +Using the parameters `loc` and `scale`, one obtains the uniformly distributed random variates on [loc, loc + scale]. Return a scalar of type `real` pseudorandom number or a rank one array of such numbers with a size of `array_size`. + +### Syntax + +`result = [[stdlib_stats_distribution(module):uniform_distribution_rvs(interface)]](loc, scale, [array_size])` + +### Arguments + +`array_size`: optional argument has intent `in` and is a scalar of type `integer`. + +`loc`: has intent `in` and is a scalar of type `real`. + +`scale`: has intent `in` and is a scalar of type `real`. + +### Return value + +The result is a scalar of type `real`, or is a rank one array of type `real` with size of `array_size`. + +### Example + +```fortran +program demo_uniform_rvs + use stdlib_stats_distribution, only: uniform => uniform_distribution_rvs + implicit none + real :: loc, scale, rv(10) + + loc = 0.0 + scale = 1.0 + print *, uniform(loc, scale) ! single standard uniform random variate in [0, 1] + rv = uniform(loc, scale, 10) ! an array of 10 uniform random variates in [0, 1] + loc = -1.0 + scale = 2.0 + print *, uniform(loc, scale) ! single uniform random variate in [-1,1] + loc = -0.5 + scale = 1.0 + rv = uniform(loc, scale, 10) ! an array of 10 uniform variates in [-0.5, 0.5] +end program demo_uniform_rvs +``` + +## `uniform_distribution_pdf` - uniform probability density function + +### Status + +Experimental + +### Description + +The probability density function of the continuous uniform distribution. + +![equation](https://latex.codecogs.com/gif.latex?f(x)=\begin{cases}\frac{1}{scale}&loc\leqslant&space;xloc+scale\end{cases}) + +### Syntax + +`result = [[stdlib_stats_distribution(module):uniform_distribution_pdf(interface)]](x, loc, scale)` + +### Arguments + +`x`: has intent `in` and is a scalar or an array of type `real`. + +`loc`: has intent `in` and is a scalar of type `real`. + +`scale`: has intent `in` and is a scalar of type `real`. + +### Return value + +The result is a scalar of type `real` or an array of type `real` with a shape conformable to input `x`. + +### Example + +```fortran +program demo_uniform_pdf + use stdlib_stats_distribution, uniform_pdf => uniform_distribution_pdf, & + uniform => uniform_distribution_rvs + implicit none + real :: loc, scale, x(3,4) + + loc = 0.0 + scale = 1.0 + print *, uniform_pdf(0.5, loc, scale) ! a probability density at 0.5 in [0,1] + print *, uniform_pdf(0.7, -1.0, 2.0) ! a probability density at 0.7 in [-1,1] + x = reshape(uniform(loc, scale, 12),[3,4]) ! uniform random variates array in [0,1] + print *, uniform_pdf(x, loc, scale) ! probability density array in [0,1] +end program demo_uniform_pdf + +``` + +## `uniform_distribution_cdf` - uniform cumulative distribution function + +### Status + +Experimental + +### Description + +Cumulative distribution function of the uniform continuous distribution + +![equation](https://latex.codecogs.com/gif.latex?F(x)=\begin{cases}\0&x uniform_distribution_cdf, & + uniform => uniform_distribution_rvs + implicit none + real :: loc, scale, x(3,4) + +loc = 0.0 + scale = 1.0 + print *, uniform_cdf(0.5, loc, scale) ! a cumulative at 0.5 in [0,1] + print *, uniform_cdf(0.7, -1.0, 2.0) ! a cumulative at 0.7 in [-1,1] + x = reshape(uniform(loc, scale, 12),[3,4]) ! uniform random variates array in [0,1] + print *, uniform_cdf(x, loc, scale) ! cumulative array in [0,1] +end program demo_uniform_cdf + +``` + +## `normal_distribution_rvs` - 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. + +The function return a scalar or a rank one array of type `real` pseudorandom number. + +### Syntax + +`result = [[stdlib_stats_distribution(module):normal_distribution_rvs(interface)]](loc, scale [, array_size])` + +### Arguments + +`array_size`: optional argument has intent `in` and is a scalar of type `integer`. + +`loc`: has intent `in` and is a scalar of type `real`. + +`scale`: has intent `in` and is a scalar of type `real`. + +### Return value + +The result is a scalar of type `real`, or is a rank one array of type `real` with size of `array_size`. + +### Example + +```fortran +program demo_normal_rvs + use stdlib_stats_distribution, only: norm => normal_distribution_rvs + implicit none + real :: loc, scale, rv(10) + + loc = 0.0 + scale = 1.0 + print *, norm(loc, scale) ! single standard normal random variate + rv = norm(loc, scale, 10) ! an array of 10 standard norml random variates + print *, norm(-1.0, 2.0) ! single norma random variate with \mu=-1, \sigma=2 + rv = norm(1.0, 1.0, 10) ! an array of 10 normal variates with \mu=1, \sigma=1 +end program demo_normal_rvs +``` + +## `normal_distribution_pdf` - normal probability density function + +### Status + +Experimental + +### Description + +The probability density function of the continuous normal distribution. + +![equation](https://latex.codecogs.com/gif.latex?f(x)=\frac{1}{\sigma&space;\sqrt{2&space;\pi}}&space;e^{-\frac{1}{2}(\frac{x-\mu}{\sigma})^{2}}) + +### Syntax + +`result = [[stdlib_stats_distribution(module):normal_distribution_pdf(interface)]](x, loc, scale)` + +### Arguments + +`x`: has intent `in` and is a scalar or an array of type `real`. + +`loc`: has intent `in` and is a scalar of type `real`. + +`scale`: has intent `in` and is a scalar of type `real`. + +### Return value + +The result is a scalar of type `real` or an array of type `real` with a shape conformable to input `x`. + +### Example + +```fortran +program demo_normal_pdf + use stdlib_stats_distribution, norm_pdf => normal_distribution_pdf, & + norm => normal_distribution_rvs + implicit none + real :: loc, scale, x(3,4) + + loc = 0.0 + scale = 1.0 + print *, norm_pdf(1.0,loc, scale) ! a probability density at 1.0 in standard normal + print *, norm_pdf(2.0,-1.0, 2.0) ! a probability density at 2.0 with \mu=-1.0 \sigma=2.0 + x = reshape(norm(0.0, 1.0, 12),[3,4]) ! standard normal random variates array + print *, norm_pdf(x, 0.0, 1.0) ! standard normal probability density array +end program demo_normal_pdf +``` + +## `normal_distribution_cdf` - normal cumulative distribution function + +### Status + +Experimental + +### Description + +Cumulative distribution function of the normal continuous distribution + +![equation](https://latex.codecogs.com/gif.latex?F(X)=\frac{1}{2}\left&space;[&space;1&space;+&space;erf(\frac{x-\mu}{\sqrt{2}&space;\sigma})&space;\right&space;]) + +### Syntax + +`result = [[stdlib_stats_distribution(module):normal_distribution_cdf(interface)]](x, loc, scale)` + +### Arguments + +`x`: has intent `in` and is a scalar or an array of type `real`. + +`loc`: has intent `in` and is a scalar of type `real`. + +`scale`: has intent `in` and is a scalar of type `real`. + +### Return value + +The result is a scalar of type `real` or an array of type `real` with a shape conformable to input `x`. + +### Example + +```fortran +program demo_norm_cdf + use stdlib_stats_distribution, norm_cdf => normal_distribution_cdf, & + norm => normal_distribution_rvs + implicit none + real :: loc, scale, x(3,4) + + print *, norm_cdf(1.0, 0.0, 1.0) ! a standard normal cumulative at 1.0 + print *, norm_cdf(2.0, -1.0, 2.0) ! a cumulative at 2.0 with \mu=-1 \sigma=2 + x = reshape(norm(0.0, 1.0, 12),[3,4]) ! standard normal random variates array + print *, norm_cdf(x, 0.0, 1.0) ! standard normal cumulative array +end program demo_norm_cdf + +``` +## `binomial_distribution_rvs` - binomial distribution random variates + +### Status + +Experimental + +### Description + +A binomial discrete random variate distribution. It is used to characterize the number of successes in a sequence of n Bernoulli trials, each with a same probablity of p. + +For a single trial, binomial distribution is Bernoulli distribution. + +### Syntax + +`result = [[stdlib_stats_distribution(module):binomial_distribution_rvs(interface)]](n, p, [array_size])` + +### Arguments + +`n`: has intent `in` and is a scalar of type `integer`. + +`p`: has intent `in` and is a scalar of type `real` with single precision. + +`array_size`: optional argument has intent `in` and is a scalar of type `integer`. + +### Return value + +The result is a scalar of type `integer`, or is a rank one array of type `integer` with a size of `array_size`. + +### Example + +```fortran +program demo_binomial_rvs + use stdlib_stats_distribution, only: binom => binomial_distribution_rvs + implicit none + integer :: n + real :: p, rv(10) + + n = 20 + p = 0.3 + print *, binom(n, p) ! single binomial random variate + rv = binom(n, p, 10) ! an array of 10 standard binomial random variates +end program demo_binomial_rvs +``` + +## `binomial_distribution_pmf` - Binomial probability mass function + +### Status + +Experimental + +### Description + +The probability mass function of the discrete binomial distribution. + +![equation](https://latex.codecogs.com/gif.latex?\begin{align*}p(k)=&\binom{n}{k}p^{k}q^{n-k},k=0,1,2,\cdots,n\\&0\leqslant&space;p\leqslant&space;1,q=1-p\end{align}) + +### Syntax + +`result = [[stdlib_stats_distribution(module):binomial_distribution_pmf(interface)]](k, n, p)` + +### Arguments + +`k`: has intent `in` and is a scalar or an array of type `integer`. + +`n`: has intent `in` and is a scalar of type `integer`. + +`p`: has intent `in` and is a scalar of type `real`. + +### Return value + +The result is a scalar of type `real` or an array of type `real` with a shape conformable to input `k`. + +### Example + +```fortran +program demo_binom_pmf + use stdlib_stats_distribution, only: binom_pmf => binomial_distribution_pmf + implicit none + real :: p + integer :: n, m(2,3) + + n = 20 + p = 0.4 + m = reshape(source=[1,2,3,4,5,6], shape=[2,3]) + print *, binom_pmf(5, n, p) ! a probability density for 5 in binomial + print *, binom_pmf(m, n, p) ! binomial probability density array +end program demo_binom_pmf +``` + +## `binomial_distribution_cdf` - Binomial cumulative distribution function + +### Status + +Experimental + +### Description + +The cumuative distribution function of the discrete binomial distribution. + +![equation](https://latex.codecogs.com/gif.latex?\begin{align*}F(k)=&\sum_{i=0}^{k}\binom{n}{i}p^{i}q^{n-i},k=0,1,2,\cdots,n\\&0\leqslant&space;p\leqslant&space;1,q=1-p\end{align}) + +### Syntax + +`result = [[stdlib_stats_distribution(module):binomial_distribution_cdf(interface)]](k, n, p)` + +### Arguments + +`k`: has intent `in` and is a scalar or an array of type `integer`. + +`n`: has intent `in` and is a scalar of type `integer`. + +`p`: has intent `in` and is a scalar of type `real`. + +### Return value + +The result is a scalar of type `real` or an array of type `real` with a shape conformable to input `k`. + +### Example + +```fortran +program demo_binom_cdf + use stdlib_stats_distribution, only: binom_cdf => binomial_distribution_cdf + implicit none + real :: p + integer :: n, m(3,2) + + n = 20 + p = 0.4 + m = reshape(source=[1,2,3,4,5,6], shape=[3,2]) + print *, binom_cdf(5, n, p) ! total probability for k not greater than 5 + print *, binom_cdf(m, n, p) ! binomial cumulative probability array +end program demo_binom_cdf +``` diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index ea7403663..94c4a5f1d 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -15,6 +15,9 @@ set(fppFiles stdlib_quadrature.fypp stdlib_quadrature_trapz.fypp stdlib_quadrature_simps.fypp + stdlib_stats_distribution.fypp + stdlib_stats_distribution_rvs.fypp + stdlib_stats_distribution_implementation.fypp ) diff --git a/src/Makefile.manual b/src/Makefile.manual index 1c731b9bb..cbfab5671 100644 --- a/src/Makefile.manual +++ b/src/Makefile.manual @@ -12,7 +12,10 @@ SRC = f18estop.f90 \ stdlib_stats.f90 \ stdlib_stats_mean.f90 \ stdlib_stats_moment.f90 \ - stdlib_stats_var.f90 + stdlib_stats_var.f90 \ + stdlib_stats_distribution.f90 \ + stdlib_stats_distribution_rvs.f90 \ + stdlib_stats_distribution_implementation.f90 LIB = libstdlib.a @@ -61,7 +64,13 @@ stdlib_stats_var.o: \ stdlib_optval.o \ stdlib_kinds.o \ stdlib_stats.o - +stdlib_stats_distribution_rvs.o: stdlib_kinds.o +stdlib_stats_distribution.o: \ + stdlib_error.o \ + stdlib_kinds.o \ + stdlib_stats_distribution_rvs.o +stdlib_stats_distribution_implementation.o: \ + stdlib_stats_distribution.o # Fortran sources that are built from fypp templates stdlib_io.f90: stdlib_io.fypp stdlib_linalg.f90: stdlib_linalg.fypp @@ -71,3 +80,7 @@ stdlib_stats.f90: stdlib_stats.fypp stdlib_stats_mean.f90: stdlib_stats_mean.fypp stdlib_stats_moment.f90: stdlib_stats_moment.fypp stdlib_stats_var.f90: stdlib_stats_var.fypp +stdlib_stats_distribution_rvs.f90: stdlib_stats_distribution_rvs.fypp +stdlib_stats_distribution.f90: stdlib_stats_distribution.fypp +stdlib_stats_distribution_implementation.f90: \ + stdlib_stats_distribution_implementation.fypp \ No newline at end of file diff --git a/src/stdlib_stats_distribution.fypp b/src/stdlib_stats_distribution.fypp new file mode 100644 index 000000000..61dcb721e --- /dev/null +++ b/src/stdlib_stats_distribution.fypp @@ -0,0 +1,218 @@ +#:include "common.fypp" +module stdlib_stats_distribution + !! Stdlib_stats_distribution + !! This module contains a large number of probability distribution and statistical functions + !! + !! Continuouse Distributions + !! Uniform distribution + !! Normal distribution + !! + !! Discrete Distribution + !! Binormial distribution + !! + use stdlib_error, only : error_stop + use stdlib_kinds + use stdlib_stats_distribution_rvs + implicit none + private + #:for k1, t1 in REAL_KINDS_TYPES + ${t1}$, parameter :: sqrt_2_pi_${t1[0]}$${k1}$ = sqrt(2.0_${k1}$ * acos(-1.0_${k1}$)) + ${t1}$, parameter :: sqrt_2_${t1[0]}$${k1}$ = sqrt(2.0_${k1}$) + #:endfor + + public :: random_seed + public :: uniform_distribution_rvs, uniform_distribution_pdf, uniform_distribution_cdf + public :: normal_distribution_rvs, normal_distribution_pdf, normal_distribution_cdf + public :: binomial_distribution_rvs, binomial_distribution_pmf, binomial_distribution_cdf + + interface random_seed + !! version: experimental + !! + !! Set or get seed value of probability distributions for pseudorandom bumber generator + !! ([Specification](../page/specs/stdlib_stats_distribution.html#description)) + !! + procedure random_distribution_seed + end interface random_seed + + interface uniform_distribution_rvs + !! version: experimental + !! + !! Uniform Distribution Random Variates + !! + !! ([Specification](../page/specs/stdlib_stats_distribution.html#description)) + !! + + #:for k1, t1 in REAL_KINDS_TYPES + module function unif_dist_rvs_${t1[0]}$${k1}$(loc, scale) result(res) + ${t1}$ :: res + ${t1}$, intent(in) :: loc, scale + end function unif_dist_rvs_${t1[0]}$${k1}$ + + #:endfor + #:for k1, t1 in REAL_KINDS_TYPES + #:for k2, t2 in INT_KINDS_TYPES + module function unif_dist_rvs_array_${t1[0]}$${k1}$_${k2}$(loc, scale, & + array_size) result(res) + ${t2}$, intent(in) :: array_size + ${t1}$, allocatable :: res(:) + ${t1}$, intent(in) :: loc, scale + end function unif_dist_rvs_array_${t1[0]}$${k1}$_${k2}$ + + #:endfor + #:endfor + end interface uniform_distribution_rvs + + interface uniform_distribution_pdf + !! version: experimental + !! + !! Uniform Distribution Probability Density Function + !! ([Specification](../page/specs/stdlib_stats_distribution.html#description)) + !! + #:for k1, t1 in REAL_KINDS_TYPES + elemental module function unif_dist_pdf_${t1[0]}$${k1}$(x, loc, scale) result(res) + ${t1}$, intent(in) :: x + ${t1}$ :: res + ${t1}$, intent(in) :: loc, scale + end function unif_dist_pdf_${t1[0]}$${k1}$ + + #:endfor + end interface uniform_distribution_pdf + + interface uniform_distribution_cdf + !! version: experimentl + !! + !! Uniform Distribution Cumulative Distribution Function + !! ([Specification](../page/specs/stdlib_stats_distribution.html#description)) + !! + #:for k1, t1 in REAL_KINDS_TYPES + elemental module function unif_dist_cdf_${t1[0]}$${k1}$(x, loc, scale) result(res) + ${t1}$, intent(in) :: x + ${t1}$ :: res + ${t1}$, intent(in) :: loc, scale + end function unif_dist_cdf_${t1[0]}$${k1}$ + + #:endfor + end interface uniform_distribution_cdf + + + interface normal_distribution_rvs + !! version: experimental + !! + !! Normal Distribution Random Variates + !! + !! ([Specification](../page/specs/stdlib_stats_distribution.html#description)) + !! + #:for k1, t1 in REAL_KINDS_TYPES + module function norm_dist_rvs_${t1[0]}$${k1}$(loc, scale) result(res) + ${t1}$ :: res + ${t1}$, intent(in) :: loc, scale + end function norm_dist_rvs_${t1[0]}$${k1}$ + + #:endfor + #:for k1, t1 in REAL_KINDS_TYPES + #:for k2, t2 in INT_KINDS_TYPES + module function norm_dist_rvs_array_${t1[0]}$${k1}$_${k2}$(loc, scale, & + array_size) result(res) + ${t2}$, intent(in) :: array_size + ${t1}$, allocatable :: res(:) + ${t1}$, intent(in) :: loc, scale + end function norm_dist_rvs_array_${t1[0]}$${k1}$_${k2}$ + + #:endfor + #:endfor + end interface normal_distribution_rvs + + interface normal_distribution_pdf + !! version: experimental + !! + !! Normal Distribution Probability Density Function + !! ([Specification](../page/specs/stdlib_stats_distribution.html#description)) + !! + #:for k1, t1 in REAL_KINDS_TYPES + elemental module function norm_dist_pdf_${t1[0]}$${k1}$(x, & + loc, scale) result(res) + ${t1}$, intent(in) :: x + ${t1}$ :: res + ${t1}$, intent(in) :: loc, scale + end function norm_dist_pdf_${t1[0]}$${k1}$ + + #:endfor + end interface normal_distribution_pdf + + interface normal_distribution_cdf + !! version: experimentl + !! + !! Normal Distribution Cumulative Distribution Function + !! ([Specification](../page/specs/stdlib_stats_distribution.html#description)) + !! + #:for k1, t1 in REAL_KINDS_TYPES + elemental module function norm_dist_cdf_${t1[0]}$${k1}$(x, & + loc, scale) result(res) + ${t1}$, intent(in) :: x + ${t1}$ :: res + ${t1}$, intent(in) :: loc, scale + end function norm_dist_cdf_${t1[0]}$${k1}$ + + #:endfor + end interface normal_distribution_cdf + + + interface binomial_distribution_rvs + !! version: experimental + !! + !! Binomial Distribution Random Variates + !! + !! ([Specification](../page/specs/stdlib_stats_distribution.html#description)) + !! + #:for k1, t1 in INT_KINDS_TYPES + module function binom_dist_rvs_${t1[0]}$${k1}$(n, p) result(res) + ${t1}$, intent(in) :: n + ${t1}$ :: res + real, intent(in) :: p + end function binom_dist_rvs_${t1[0]}$${k1}$ + + #:endfor + #:for k1, t1 in INT_KINDS_TYPES + module function binom_dist_rvs_array_${t1[0]}$${k1}$(n, p, & + array_size) result(res) + ${t1}$, intent(in) :: array_size, n + ${t1}$, allocatable :: res(:) + real, intent(in) :: p + end function binom_dist_rvs_array_${t1[0]}$${k1}$ + + #:endfor + end interface binomial_distribution_rvs + + interface binomial_distribution_pmf + !! version: experimental + !! + !! Normal Distribution Probability Density Function + !! ([Specification](../page/specs/stdlib_stats_distribution.html#description)) + !! + #:for k1, t1 in INT_KINDS_TYPES + elemental module function binom_dist_pmf_${t1[0]}$${k1}$(k, n, p) result(res) + ${t1}$, intent(in) :: k, n + real, intent(in) :: p + real :: res + end function binom_dist_pmf_${t1[0]}$${k1}$ + + #:endfor + end interface binomial_distribution_pmf + + interface binomial_distribution_cdf + !! version: experimentl + !! + !! Normal Distribution Cumulative Distribution Function + !! ([Specification](../page/specs/stdlib_stats_distribution.html#description)) + !! + #:for k1, t1 in INT_KINDS_TYPES + elemental module function binom_dist_cdf_${t1[0]}$${k1}$(k, n, p) result(res) + ${t1}$, intent(in) :: k, n + real, intent(in) :: p + real :: res + end function binom_dist_cdf_${t1[0]}$${k1}$ + + #:endfor + end interface binomial_distribution_cdf + +end module stdlib_stats_distribution diff --git a/src/stdlib_stats_distribution_implementation.fypp b/src/stdlib_stats_distribution_implementation.fypp new file mode 100644 index 000000000..63dd28f0f --- /dev/null +++ b/src/stdlib_stats_distribution_implementation.fypp @@ -0,0 +1,149 @@ +#:include "common.fypp" +submodule (stdlib_stats_distribution) stdlib_stats_distribution_imp + + contains + + #:for k1, t1 in REAL_KINDS_TYPES + module procedure unif_dist_rvs_${t1[0]}$${k1}$ + if(scale == 0.0_${k1}$) call error_stop("Error: scale parameter is zero") + res = loc + scale * uni(scale) + end procedure unif_dist_rvs_${t1[0]}$${k1}$ + + #:endfor + #:for k1, t1 in REAL_KINDS_TYPES + #:for k2, t2 in INT_KINDS_TYPES + module procedure unif_dist_rvs_array_${t1[0]}$${k1}$_${k2}$ + ${t2}$ :: i + + if(scale == 0.0_${k1}$) call error_stop("Error: scale parameter is zero") + allocate(res(array_size)) + do i = 1_${k2}$, array_size + res(i) = loc + scale * uni(scale) + end do + end procedure unif_dist_rvs_array_${t1[0]}$${k1}$_${k2}$ + + #:endfor + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + module procedure unif_dist_pdf_${t1[0]}$${k1}$ + if(scale == 0.0_${k1}$) then + res = 1.0_${k1}$ + elseif(x < min(loc, loc + scale) .or. x > max(loc, loc + scale)) then + res = 0.0_${k1}$ + else + res = 1.0_${k1}$ / scale + end if + end procedure unif_dist_pdf_${t1[0]}$${k1}$ + + #:endfor + #:for k1, t1 in REAL_KINDS_TYPES + module procedure unif_dist_cdf_${t1[0]}$${k1}$ + if(scale == 0.0_${k1}$) then + res = 1.0_${k1}$ + elseif(x < min(loc, loc + scale) ) then + res = 0.0_${k1}$ + elseif(x >= min(loc, loc + scale) .and. x <= max(loc, loc + scale)) then + res = (x - loc) / scale + else + res = 1.0_${k1}$ + end if + end procedure unif_dist_cdf_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + module procedure norm_dist_rvs_${t1[0]}$${k1}$ + res = rnor(scale) * scale + loc + end procedure norm_dist_rvs_${t1[0]}$${k1}$ + + #:endfor + #:for k1, t1 in REAL_KINDS_TYPES + #:for k2, t2 in INT_KINDS_TYPES + module procedure norm_dist_rvs_array_${t1[0]}$${k1}$_${k2}$ + ${t2}$ :: i + + allocate(res(array_size)) + do i = 1_${k2}$, array_size + res(i) = rnor(scale) * scale + loc + end do + end procedure norm_dist_rvs_array_${t1[0]}$${k1}$_${k2}$ + + #:endfor + #:endfor + #:for k1, t1 in REAL_KINDS_TYPES + module procedure norm_dist_pdf_${t1[0]}$${k1}$ + res = exp(- 0.5_${k1}$ * (x - loc) * (x - loc) / (scale * scale)) / & + (sqrt_2_Pi_${t1[0]}$${k1}$ * scale) + end procedure norm_dist_pdf_${t1[0]}$${k1}$ + + #:endfor + #:for k1, t1 in REAL_KINDS_TYPES + module procedure norm_dist_cdf_${t1[0]}$${k1}$ + res = (1.0_${k1}$ + erf((x - loc) / (scale * sqrt_2_${t1[0]}$${k1}$))) & + / 2.0_${k1}$ + end procedure norm_dist_cdf_${t1[0]}$${k1}$ + + #:endfor + + + #:for k1, t1 in INT_KINDS_TYPES + module procedure binom_dist_rvs_${t1[0]}$${k1}$ + + if(n * p < 10.0) then + res = random_binomial_num1(n, p) + else + res = random_binomial_num2(n, p) + endif + end procedure binom_dist_rvs_${t1[0]}$${k1}$ + + #:endfor + #:for k1, t1 in INT_KINDS_TYPES + module procedure binom_dist_rvs_array_${t1[0]}$${k1}$ + ${t1}$ :: i + + allocate(res(array_size)) + if(n * p < 10.0) then + do i =1_${k1}$, array_size + res(i) = random_binomial_num1(n, p) + end do + else + do i = 1_${k1}$, array_size + res(i) = random_binomial_num2(n, p) + end do + endif + end procedure binom_dist_rvs_array_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in INT_KINDS_TYPES + module procedure binom_dist_pmf_${t1[0]}$${k1}$ + res = log_factorial(n) - & + log_factorial(k) - & + log_factorial(n - k) + & + k * log(p) + (n - k) * log(1.0 - p) + res = exp(res) + end procedure binom_dist_pmf_${t1[0]}$${k1}$ + + #:endfor + #:for k1, t1 in INT_KINDS_TYPES + module procedure binom_dist_cdf_${t1[0]}$${k1}$ + ${t1}$ :: i + real(dp) :: coeff, q, logpmf, sum, lq, lpq + + q = 1.0_dp - p + lq = log(q) + lpq = log(p / q) + coeff = 0.0_dp + sum = exp(n * lq) + do i = 1, k + coeff = coeff + log(real(n - i + 1, kind=dp)) - log(real(i, kind=dp)) + logpmf = coeff + i * lpq + n * lq + sum = sum + exp(logpmf) + end do + res = sum + end procedure binom_dist_cdf_${t1[0]}$${k1}$ + + #:endfor + +end submodule stdlib_stats_distribution_imp \ No newline at end of file diff --git a/src/stdlib_stats_distribution_rvs.fypp b/src/stdlib_stats_distribution_rvs.fypp new file mode 100644 index 000000000..9a4d25fb1 --- /dev/null +++ b/src/stdlib_stats_distribution_rvs.fypp @@ -0,0 +1,791 @@ +#:include "common.fypp" +Module stdlib_stats_distribution_rvs + use stdlib_kinds + implicit none + private + type :: binom_table + character(100) :: name + real(dp), allocatable :: tb(:) + end type binom_table + real(dp), parameter :: MESENNE_NUMBER = 1.0_dp / (2.0_dp ** 53 - 1.0_dp) + real(dp), parameter :: HALF = 0.5_dp, ONE = 1.0_dp, TWO = 2.0_dp, & + TWO_PI = TWO * acos(-one) + integer(int32), save :: kn(0:127), ke(0:255) + real(dp), save :: wn(0:127), fn(0:127), we(0:255), fe(0:255) + integer(int64), save :: st(4), si = 614872703977525537_int64 + logical, save :: zig_initialized = .false. + integer, save :: num_binom_tb = 0 + type(binom_table), allocatable, save :: binom_tb(:) + + public :: uni, rnor + public :: random_binomial_num1, random_binomial_num2 + public :: log_factorial, random_distribution_seed + + + interface uni + ! Generation of uniform distributed variates. + ! Based on the paper by Frederic Goualard, "Generating Random Floating- + ! Point Numbers By Dividing Integers: a Case Study", Proceedings of + ! ICCS 2020, June 20202, Amsterdam, Netherlands + ! + ! Fortran 90 by Jim-215-Fisher + ! + #:for k1, t1 in REAL_KINDS_TYPES + module procedure uni_${t1[0]}$${k1}$ + #:endfor + end interface uni + + interface rnor + ! Generation of normally distributed variates. + ! Marsaglia & Tsang generator for random normals & random exponentials. + ! Translated from C by Alan Miller (amiller@bigpond.net.au) + ! + ! Marsaglia, G. & Tsang, W.W. (2000) `The ziggurat method for generating + ! random variables', J. Statist. Software, v5(8). + ! + ! Latest version - 1 January 2001 + ! + #:for k1, t1 in REAL_KINDS_TYPES + module procedure rnor_${t1[0]}$${k1}$ + #:endfor + end interface rnor + + interface rexp + ! Generation of exponentially distributed variates. + ! Marsaglia & Tsang generator for random normals & random exponentials. + ! Translated from C by Alan Miller (amiller@bigpond.net.au) + ! + ! Latest version - 1 January 2001 + ! + #:for k1, t1 in REAL_KINDS_TYPES + module procedure rexp_${t1[0]}$${k1}$ + #:endfor + end interface rexp + + interface random_binomial_num1 + ! Generation of binomially distributed variates for np < 10 + ! Inverse transformation of binomial distribution function + ! Algorithm BINV + ! + #:for k1, t1 in INT_KINDS_TYPES + module procedure random_binomial_num1_${t1[0]}$${k1}$ + #:endfor + end interface random_binomial_num1 + + interface random_binomial_num2 + ! Generation of binomially distributed variates for np>=10 + #:for k1, t1 in INT_KINDS_TYPES + module procedure random_binomial_num2_${t1[0]}$${k1}$ + #:endfor + end interface random_binomial_num2 + + interface fc_stirling + ! Table and function for log factorial stirling approximation correction + ! term. + ! + #:for k1, t1 in INT_KINDS_TYPES + module procedure fc_stirling_${t1[0]}$${k1}$ + #:endfor + end interface fc_stirling + + interface log_factorial + ! Table and function for log factorial. + ! + #:for k1, t1 in INT_KINDS_TYPES + module procedure log_factorial_${t1[0]}$${k1}$ + #:endfor + end interface log_factorial + + + contains + + subroutine zigset + ! Marsaglia & Tsang generator for random normals & random exponentials. + ! Translated from C by Alan Miller (amiller@bigpond.net.au) + ! + ! 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 + ! + ! N.B. It is assumed that all integers are 32-bit. + ! N.B. The value of M2 has been halved to compensate for the lack of + ! unsigned integers in Fortran. + ! + ! Latest version - 1 January 2001 + ! + real(dp), parameter :: M1 = 2147483648.0_dp, M2 = 2147483648.0_dp + real(dp) :: dn = 3.442619855899_dp, tn, & + vn = 0.00991256303526217_dp, & + de = 7.697117470131487_dp, te, & + ve = 0.003949659822581572_dp, q + integer(int32) :: i + + tn = dn + te = de + ! 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) = 1.0 + 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 + + ! tables for random exponetials + q = ve * exp( de ) + ke(0) = int((de / q) * M2, kind = int32) + ke(1) = 0 + we(0) = q / M2 + we(255) = de / M2 + fe(0) = 1.0 + fe(255) = exp( -de ) + do i = 254, 1, -1 + de = -log( ve / de + exp( -de ) ) + ke(i+1) = int(M2 * (de / te), kind = int32) + te = de + fe(i) = exp( -de ) + we(i) = de / M2 + end do + zig_initialized = .true. + return + end subroutine zigset + + +! function shr3( ) result( ival ) +! ! generate random 32-bit integers +! integer(int32) :: ival +! integer(int32) :: jz +! +! jz = jsr +! jsr = ieor( jsr, ishft( jsr, 13 ) ) +! jsr = ieor( jsr, ishft( jsr, -17 ) ) +! jsr = ieor( jsr, ishft( jsr, 5 ) ) +! ival = jz + jsr +! return +! end function shr3 + + function xoshiro256ss( ) result (res) + ! Generate random 64-bit integers + ! Written in 2018 by David Blackman and Sebastiano Vigna (vigna@acm.org) + ! + ! This is xoshiro256** 1.0, one of our all-purpose, rock-solid + ! generators. It has excellent (sub-ns) speed, a state (256 bits) that is + ! large enough for any parallel application, and it passes all tests we + ! are aware of. + ! + ! The state must be seeded so that it is not everywhere zero. If you have + ! a 64-bit seed, we suggest to seed a splitmix64 generator and use its + ! output to fill st. + ! + ! Fortran 90 version translated from C by Jim-215-Fisher + ! + integer(int64) :: res, t + + res = rol64(st(2) * 5 , 7) * 9 + t = shiftl(st(2), 17) + st(3) = ieor(st(3), st(1)) + st(4) = ieor(st(4), st(2)) + st(2) = ieor(st(2), st(3)) + st(1) = ieor(st(1), st(4)) + st(3) = ieor(st(3), t) + st(4) = rol64(st(4), 45) + return + + contains + + function rol64(x, k) result(res) + integer(int64), intent(in) :: x + integer(int32), intent(in) :: k + integer(int64) :: t1, t2, res + + t1 = shiftr(x, (64 - k)) + t2 = shiftl(x, k) + res = ior(t1, t2) + end function rol64 + end function xoshiro256ss + + function splitmix64(s) result(res) + ! Written in 2015 by Sebastiano Vigna (vigna@acm.org) + ! This is a fixed-increment version of Java 8's SplittableRandom + ! generator. + ! See http://dx.doi.org/10.1145/2714064.2660195 and + ! http://docs.oracle.com/javase/8/docs/api/java/util/SplittableRandom.html + ! + ! It is a very fast generator passing BigCrush, and it can be useful if + ! for some reason you absolutely want 64 bits of state. + ! + ! Fortran 90 translated from C by Jim-215-Fisher + ! + integer(int64) :: res, int01, int02, int03 + integer(int64), intent(in), optional :: s + data int01, int02, int03/-7046029254386353131_int64, & + -4658895280553007687_int64, & + -7723592293110705685_int64/ + + if(present(s)) si = s + res = si + si = res + int01 + res = ieor(res, shiftr(res, 30)) * int02 + res = ieor(res, shiftr(res, 27)) * int03 + res = ieor(res, shiftr(res, 31)) + return + end function splitmix64 + + subroutine random_distribution_seed(put, get) + ! Random seed for distribution + integer(int32), intent(in) :: put + integer(int32), intent(out) :: get + integer(int64) :: tmp + integer(int32) :: i + + tmp = splitmix64(int(put, kind = int64)) + do i = 1, 10 + tmp = splitmix64( ) + end do + do i = 1, 4 + tmp = splitmix64( ) + st(i) = tmp + end do + get = int(tmp, kind = int32) + end subroutine random_distribution_seed + + #:for k1, t1 in REAL_KINDS_TYPES + function uni_${t1[0]}$${k1}$( x ) result( fn_val ) + ! Based on the paper by Frederic Goualard, "Generating Random Floating- + ! Point Numbers By Dividing Integers: a Case Study", Proceedings of + ! ICCS 2020, June 20202, Amsterdam, Netherlands + ! + ${t1}$, intent(in) :: x + ${t1}$ :: fn_val + integer(int64) :: tmp + + tmp = shiftr(xoshiro256ss( ), 11) + fn_val = real(tmp * MESENNE_NUMBER, kind = ${k1}$) + return + end function uni_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + function rnor_${t1[0]}$${k1}$( z ) result( fn_val ) + ! Normal random variate + ${t1}$, intent(in) :: z + ${t1}$ :: fn_val + ${t1}$, parameter :: r = 3.442619855899_${k1}$ + ${t1}$ :: x, y + integer(int32) :: hz, iz + + if( .not. zig_initialized ) call zigset + + ! original algorithm use 32bit + hz = shiftr(xoshiro256ss( ), 32) + + iz = iand( hz, 127 ) + if( abs( hz ) < kn(iz) ) then + fn_val = hz * wn(iz) + else + do + if( iz == 0 ) then + do + x = -log( uni(z) ) / r + y = -log( uni(z) ) + if( y + y >= x * x ) exit + end do + fn_val = r + x + if( hz <= 0 ) fn_val = -fn_val + return + end if + x = hz * wn(iz) + if( fn(iz) + uni(z) * (fn(iz-1) - fn(iz)) < exp(-HALF * & + x * x) ) then + fn_val = x + return + end if + + !original algorithm use 32bit + hz = shiftr(xoshiro256ss( ), 32) + + iz = iand( hz, 127 ) + if( abs( hz ) < kn(iz) ) then + fn_val = hz * wn(iz) + return + end if + end do + end if + return + end function rnor_${t1[0]}$${k1}$ + + #:endfor + + + #:for k1, t1 in REAL_KINDS_TYPES + function rexp_${t1[0]}$${k1}$( z ) result( fn_val ) + ! exponential variate algorithm + ${t1}$, intent(in) :: z + ${t1}$ :: fn_val, x + ${t1}$ :: r = 7.69711747013104972_${k1}$ + integer(int32) :: jz, iz + + if( .not. zig_initialized ) call zigset + + ! Original algorithm use 32bit + jz = shiftr(xoshiro256ss( ), 32) + + iz = iand( jz, 255 ) + if( abs( jz ) < ke(iz) ) then + fn_val = abs(jz) * we(iz) + return + end if + do + if( iz == 0 ) then + fn_val = r - log( uni(z) ) + return + end if + x = abs( jz ) * we(iz) + if( fe(iz) + uni(z) * (fe(iz-1) - fe(iz)) < exp( -x ) ) then + fn_val = x + return + end if + + !original algorithm use 32bit + jz = shiftr(xoshiro256ss( ), 32) + + iz = iand( jz, 255 ) + if( abs( jz ) < ke(iz) ) then + fn_val = abs( jz ) * we(iz) + return + end if + end do + return + end function rexp_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in INT_KINDS_TYPES + function random_binomial_num1_${t1[0]}$${k1}$(n, p) result(res) + ! BINV algorithm for n*p < 10 + ${t1}$, intent(in) :: n + real, intent(in) :: p + ${t1}$ :: res, k + real :: q, s, a, r, u, pq + + pq = min(p, 1.0 - p) + q = 1.0 - pq + s = pq / q + a = (n + 1) * s + r = q ** n + k = 0 + u = uni(q) + do + if(u <= r) exit + u = u - r + k = k + 1_${k1}$ + r = (a / k - s) * r + end do + if(p > 0.5) then + res = n - k + else + res = k + endif + return + end function random_binomial_num1_${t1[0]}$${k1}$ + + #:endfor + + #:set INVF=1 + #:set BTRD=0 + #:set BTPE=0 + #:if BTRD + + #:for k1, t1 in INT_KINDS_TYPES + function random_binomial_num2_${t1[0]}$${k1}$(n, p) result(res) + ! BTRD algorithm, may not pass chi^2 test + ! Generation of binomial variates for np >=10 and p<=0.5 + ! Wolfgang Hormann generator for binomial random variates. + ! + ! Fortran 90 written by Jim-215-Fisher + ! + ! Hormann Wolfgang (1993) `The generation of binomial random variates', + ! Journal of Statistical Computation and Simulation, v46(1-2). + ! + ! Article can be found at http://core.ac.uk/download/pdf/11007254.pdf + ! + + real, intent(in) :: p + ${t1}$, intent(in) :: n + ${t1}$ :: res, m, i, k, km, nm, nk + real :: r, nr, npq, a, b, c, alpha, vr, urvr, u, v, us, rho, f, t, h + + m = floor((n + 1) * p) + r = p / (1.0 - p) + nr = (n + 1) * r + npq = n * p * (1.0 - p) + b = 1.15 + 2.53 * sqrt(npq) + a = -0.0873 + 0.0248 * b + 0.01 * p + c = n * p + 0.5 + alpha = (2.83 + 5.1 / b) * sqrt(npq) + vr = 0.92 - 4.2 / b + urvr = 0.86 * vr + + do + ! Steps 1, 2 + v = uni(p) + if(v <= urvr) then + u = v / vr - 0.43 + res = floor(((2 * a) / (0.5 - abs(u)) + b) * u + c) - 1 + return + endif + if(v >= vr) then + u = uni(p) - 0.5 + else + u = v / vr - 0.93 + u = sign(1.0, u) * 0.5 - u + v = uni(p) * vr + endif + !step 3.0 + us = 0.5 - abs(u) + k = floor(((2 * a) / us + b) * u + c) + if(k < 0 .or. k > n) cycle + v = v * alpha / (b + a / (us * us)) + km = abs(k - m) + if(km <= 15_${k1}$) then + !Step 3.1 + f = 1.0 + if(m < k) then + do i = m, k + f = f * (nr / i - r) + end do + elseif(m > k) then + do i = k, m + v = v * (nr / i -r) + end do + endif + if(v <= f) then + res = k - 1 + return + endif + else + !Step 3.2 + v = log(v) + rho = (km / npq) * (((km / 3.0 + 0.625) * km + 1.0 / 6.0) & + / npq + 0.5) + t = - (km * km) / (2.0 * npq) + if(v < t - rho) then + res = k - 1 + return + endif + if(v <= t + rho) then + !Step 3.3 + nm = n - m + 1_${k1}$ + h = (m + 0.5) * log((m + 1.0) / (r * nm)) + & + fc_stirling(m) + fc_stirling(n - m) + !Step 3.4 + nk = n - k + 1_${k1}$ + h = h + (n + 1) * (log(real(nm)) - log(real(nk))) + & + (k + 0.5) * log(nk * r / (k + 1)) - & + fc_stirling(k) - fc_stirling(n - k) + if(v <= h) then + res = k - 1 + return + endif + endif + endif + end do + end function random_binomial_num2_${t1[0]}$${k1}$ + + #:endfor + + #:elif BTPE + + #:for k1, t1 in INT_KINDS_TYPES + function random_binomial_num2_${t1[0]}$${k1}$(n, p) result(res) + ! BTPE algorithm, may not pass chi^2 test. + real, intent(in) :: p + ${t1}$, intent(in) :: n + ${t1}$ :: res, m, y, i, k + real :: fm, q, r, nrq, p1, p2, p3, p4, xm, xl, xr, c, a, lamdal, lamdar + real :: u, v, x, fy, s, rho, t, aa, x1, x2, f1, f2,z, w, z2, w2, cc + + q = 1.0 - p + r = min(p, q) + fm = (n + 1) * r + nrq = n * r * q + m = int(fm, kind = ${k1}$) + p1 = int(2.195 * sqrt(nrq) - 4.6 * q) + 0.5 + xm = m + 0.5 + xl = xm - p1 + xr = xm + p1 + c = 0.134 + 20.5 / (15.3 + m) + a = (fm - xl)/ (fm - xl * r) + lamdal = a * (1 + a / 2) + a = (xr - fm) / (xr * q) + lamdar = a * ( 1 + a / 2) + p2 = p1 * (1 + 2 * c) + p3 = p2 + c / lamdal + p4 = p3 + c / lamdar + do +!step 1. + u = p4 * uni(q) + v = uni(q) + if(u <= p1) then + y=int(xm - p1 * v + u, kind = ${k1}$) + else +!step 2. + if(U <= p2) then + x = xl + (u - p1) / c + v = v * c + 1 - abs(m - x + 0.5) / p1 + if(v > 1) cycle + y = int(x, kind = ${k1}$) + else +!step 3. + if(u <= p3) then + y = int(xl + log(v) / lamdal, kind = ${k1}$) + if(y < 0) cycle + v = v * (u - p2) * lamdal + else +!step 4. + y = int(xr - log(v) / lamdar, kind = ${k1}$) + if(y > n) cycle + v = v * (u - p3) * lamdar + endif + endif +!step 5.0 + k = abs(y - m) + if(.not.(k > 20 .and. k < nrq / 2 - 1)) then +!step 5.1 + s = r / q + a = (n + 1) * s + fy = 1.0 + if(m < y) then + do i = m, y + fy = fy * (a / i - s) + end do + elseif(m > y) then + do i = y, m + fy = fy / (a / i - s) + end do + endif + if(v > fy) cycle + if(p > 0.5) then + res = n - y -1 + else + res = y - 1 + endif + return + else +!step 5.2. + rho = k / nrq + rho = rho * ((k * (k / 3.0 + 0.625) + 1.0 / 6) / nrq + 0.5) + t = - 0.5 * K / nrq * k + aa = log(v) + if(aa < t - rho) then + if(p > 0.5) then + res = n - y - 1 + else + res = y - 1 + endif + return + endif + if(aa > t + rho) cycle + endif +!step 5.3 + x1 = y + 1.0 + f1 = m + 1.0 + z = n + 1.0 - m + w = n - y + 1.0 + x2 = x1 * x1 + f2 = f1 * f1 + z2 = z * z + w2 = w * w + cc = xm * log(f1 / x1) + (n - m + 0.5) *log(z / w) + cc = cc + (y - m) * log(w * r / x1 / q) + cc = cc + (13860.0 - (462.0 - (132.0 - (99.0 - 140.0 / f2) & + / f2) / f2) / f2) / f1 / 166320.0 + cc = cc + (13860.0 - (462.0 - (132.0 - (99.0 - 140.0 / z2) & + / z2) / z2) / z2) / z / 166320.0 + cc = cc + (13860.0 - (462.0 - (132.0 - (99.0 - 140.0 / x2) & + / x2) / x2) / x2) / x1 / 166320.0 + cc = cc + (13860.0 - (462.0 - (132.0 - (99.0 - 140.0 / w2) & + / w2) / w2) / w2) / w / 166320.0 + if(aa > cc) cycle + endif +!step 6. + if(p > 0.5) then + res = n - y - 1 + else + res = y - 1 + endif + return + enddo + end function random_binomial_num2_${t1[0]}$${k1}$ + #:endfor + + #:elif INVF + + #:for k1, t1 in INT_KINDS_TYPES + function random_binomial_num2_${t1[0]}$${k1}$(n, p) result(res) + ${t1}$, intent(in) :: n + real, intent(in) :: p + real(dp) :: q, lq, lpq, logpmf, coeff, sum, u + ${t1}$ :: res, k + integer :: i, iz, ix, ik + character(100) :: buffer + type(binom_table), allocatable :: temp(:) + real(dp), allocatable :: tb(:) + + + if(.not. allocated(binom_tb)) allocate(binom_tb(5)) + ! convert input n and p into string as index for binomial table + write(buffer, '(i0,f6.4)') n, p + ik = 0 + ! Search binomial table, get matched number + do i = 1, num_binom_tb + if(binom_tb(i) % name == buffer) then + ik = i + exit + endif + end do + ! If no match, create binomial table + if(ik == 0) then + num_binom_tb = num_binom_tb + 1 + ! If already have 5 more tables, extend the table by 1 + if(num_binom_tb > 5) then + allocate(temp(num_binom_tb)) + temp(1:size(binom_tb))=binom_tb + call move_alloc(temp, binom_tb) + endif + ! Set the table index by string + binom_tb(num_binom_tb) % name = buffer + allocate(binom_tb(num_binom_tb) % tb(0:n), tb(0:n)) + ! Calculate the binomial cumulative distribution function + q = 1.0 - p + lq = log(q) + lpq = log(p / q) + coeff = 0.0_dp + tb(0) = exp(n * lq) + sum = tb(0) + do i = 1, n + coeff = coeff + log(real(n - i + 1, kind=dp)) - log(real(i, kind=dp)) + logpmf = coeff + i * lpq + n * lq + sum = sum + exp(logpmf) + tb(i) = sum + end do + binom_tb(num_binom_tb) % tb = tb + ik = num_binom_tb + deallocate(tb) + endif + ! Access possible binomial variate by inverse transformation + u = uni(p) + ! Approximate position of u in cdf + iz = floor(u * n) + ix = 0 + allocate(tb, source = binom_tb(ik) % tb) + ! Exact value among n which has closeset cdf value to u + if(u < binom_tb(ik) % tb(iz)) then + do + ix = ix +1 + if(iz - ix <= 0) then + res = 0 + exit + elseif(u >= tb(iz - ix)) then + res = iz - ix + 1 + exit + end if + end do + else + do + ix = ix + 1 + if(iz + ix >= n) then + res = n + exit + elseif(u < tb(iz + ix)) then + res = iz + ix + exit + endif + end do + endif + return + end function random_binomial_num2_${t1[0]}$${k1}$ + #:endfor + #:endif + + #:for k1, t1 in INT_KINDS_TYPES + elemental module function fc_stirling_${t1[0]}$${k1}$(k) result(x) + real(dp) :: x, k1_2 + ${t1}$, intent(in) :: k + ${t1}$ :: k_1 + + select case(k) + case (0) + x = 0.08106146679532726_dp + case (1) + x = 0.04134069595540929_dp + case (2) + x = 0.02767792568499834_dp + case (3) + x = 0.02079067210376509_dp + case (4) + x = 0.01664469118982119_dp + case (5) + x = 0.01387612882307075_dp + case (6) + x = 0.01189670994589177_dp + case (7) + x = 0.01041126526197209_dp + case (8) + x = 0.009255462182712733_dp + case (9) + x = 0.008330563433362871_dp + case (10:) + k_1 = k + 1 + k1_2 = real(k_1, kind=dp) * k_1 + x = ONE / (1260.0_dp * k1_2) + x = (ONE / 360.0_dp - x) / k1_2 + x = (ONE / 12.0_dp - x) / k_1 + end select + end function fc_stirling_${t1[0]}$${k1}$ + #:endfor + + + #:for k1, t1 in INT_KINDS_TYPES + elemental module function log_factorial_${t1[0]}$${k1}$(k) result(res) + ${t1}$, intent(in) :: k + real(dp) :: res + + select case (k) + case (0) + res = 0.0_dp + case (1) + res = 0.0_dp + case (2) + res = 0.6931471805599453_dp + case (3) + res = 1.791759469228055_dp + case (4) + res = 3.178053830347946_dp + case (5) + res = 4.787491742782046_dp + case (6) + res = 6.579251212010101_dp + case (7) + res = 8.525161361065415_dp + case (8) + res = 10.604602902745251_dp + case (9) + res = 12.801827480081469_dp + case (10:) + res = HALF * log(TWO_PI) + (k + HALF) * log(k + 1.0_dp) - (k + 1) + & + fc_stirling(k) + end select + end function log_factorial_${t1[0]}$${k1}$ + #:endfor +end module stdlib_stats_distribution_rvs \ No newline at end of file diff --git a/src/tests/stats/CMakeLists.txt b/src/tests/stats/CMakeLists.txt index 36ffc7aeb..974afa526 100644 --- a/src/tests/stats/CMakeLists.txt +++ b/src/tests/stats/CMakeLists.txt @@ -5,6 +5,7 @@ ADDTEST(moment) ADDTEST(rawmoment) ADDTEST(var) ADDTEST(varn) +ADDTEST(distribution) 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 aacaf98ab..0e4193668 100644 --- a/src/tests/stats/Makefile.manual +++ b/src/tests/stats/Makefile.manual @@ -1,3 +1,4 @@ -PROGS_SRC = test_mean.f90 test_moment.f90 test_var.f90 -include ../Makefile.manual.test.mk +PROGS_SRC = test_mean.f90 test_moment.f90 test_var.f90 test_distribution.f90 + +include ../Makefile.manual.test.mk \ No newline at end of file diff --git a/src/tests/stats/test_distribution.f90 b/src/tests/stats/test_distribution.f90 new file mode 100644 index 000000000..4f99da0f4 --- /dev/null +++ b/src/tests/stats/test_distribution.f90 @@ -0,0 +1,1861 @@ +program test_distribution + use stdlib_stats_distribution, & + uni_rvs => uniform_distribution_rvs, & + uni_pdf => uniform_distribution_pdf, & + uni_cdf => uniform_distribution_cdf, & + nor_rvs => normal_distribution_rvs, & + nor_pdf => normal_distribution_pdf, & + nor_cdf => normal_distribution_cdf, & + binom_rvs => binomial_distribution_rvs, & + binom_pmf => binomial_distribution_pmf, & + binom_cdf => binomial_distribution_cdf + use stdlib_error, only : check + use stdlib_kinds + + implicit none + real(sp), parameter :: sptol = 1000 * epsilon(1.0_sp) + real(dp), parameter :: dptol = 1000 * epsilon(1.0_dp) + real(qp), parameter :: qptol = 1000 * epsilon(1.0_qp) + logical :: warn = .true. + + call test_random_seed + call test_uniform_random_generator + call test_normal_random_generator + call test_binomial_random_generator + + call test_uni_rvs_sp_int8 + call test_uni_rvs_sp_int16 + call test_uni_rvs_sp_int32 + call test_uni_rvs_sp_int64 + call test_uni_rvs_dp_int8 + call test_uni_rvs_dp_int16 + call test_uni_rvs_dp_int32 + call test_uni_rvs_dp_int64 + call test_uni_rvs_qp_int8 + call test_uni_rvs_qp_int16 + call test_uni_rvs_qp_int32 + call test_uni_rvs_qp_int64 + + call test_uni_pdf_sp + call test_uni_pdf_dp + call test_uni_pdf_qp + call test_uni_cdf_sp + call test_uni_cdf_dp + call test_uni_cdf_qp + + call test_nor_rvs_sp_int8 + call test_nor_rvs_sp_int16 + call test_nor_rvs_sp_int32 + call test_nor_rvs_sp_int64 + call test_nor_rvs_dp_int8 + call test_nor_rvs_dp_int16 + call test_nor_rvs_dp_int32 + call test_nor_rvs_dp_int64 + call test_nor_rvs_qp_int8 + call test_nor_rvs_qp_int16 + call test_nor_rvs_qp_int32 + call test_nor_rvs_qp_int64 + + call test_nor_pdf_sp + call test_nor_pdf_dp + call test_nor_pdf_qp + call test_nor_cdf_sp + call test_nor_cdf_dp + call test_nor_cdf_qp + + call test_binom_rvs_int8 + call test_binom_rvs_int16 + call test_binom_rvs_int32 + call test_binom_rvs_int64 + call test_binom_pmf_int8 + call test_binom_pmf_int16 + call test_binom_pmf_int32 + call test_binom_pmf_int64 + call test_binom_cdf_int8 + call test_binom_cdf_int16 + call test_binom_cdf_int32 + call test_binom_cdf_int64 + + + contains + + subroutine test_random_seed + integer(int32) :: put, get, res(5) + integer(int32) :: ans(5) = [-1859553078, -1933696596, -642834430, & + 1711399314, 1548311463] + integer :: i + + print *, "" + print *, "Test random_seed" + put = 135792468 + do i = 1, 5 + call random_seed(put,get) + res(i) = get + put = get + end do + call check(all(res == ans), & + msg="random seed test failed.",warn=warn) + end subroutine test_random_seed + + subroutine test_uniform_random_generator + integer :: i, j, freq(0:1000), num=10000000 + real(dp) :: chisq, expct + + print *,"" + print *, "Test uniform random generator with chi-squared" + freq = 0 + do i = 1, num + j = 1000 * uni_rvs(0.0, 1.0) + freq(j) = freq(j) + 1 + end do + chisq = 0.0_dp + expct = num / 1000 + do i = 0, 999 + chisq = chisq + (freq(i) - expct) ** 2 / expct + end do + write(*,*) "The critical values for chi-squared with 1000 d. of f. are" + write(*,*) "867.48 and 1143.92" + write(*,*) "Chi-squared for uniform random generator is : ", chisq + call check((chisq < 1143.9 .and. chisq > 867.48) , & + msg="uniform randomness failed chi-squared test", warn=warn) + end subroutine test_uniform_random_generator + + subroutine test_normal_random_generator + integer :: i, j, freq(0:1000), num=10000000 + real(dp) :: chisq, expct + + print *, "" + print *, "Test normal random generator with chi-squared" + freq = 0 + do i = 1, num + j = 1000 * (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 / 1000 + do i = 0, 999 + chisq = chisq + (freq(i) - expct) ** 2 / expct + end do + write(*,*) "The critical values for chi-squared with 1000 d. of f. are" + write(*,*) "867.48 and 1143.92" + write(*,*) "Chi-squared for normal random generator is : ", chisq + call check((chisq < 1143.9 .and. chisq > 867.48), & + msg="normal randomness failed chi-squared test", warn=warn) + end subroutine test_normal_random_generator + + subroutine test_binomial_random_generator + integer :: i, j, n, num=1000000 + real(dp) :: chisq, expct + real :: p + integer, allocatable :: freq(:) + + print *, "" + print *, "Test binomial1 random generator with chi-squared" + n = 10 + p = 0.45 + allocate(freq(0:n)) + freq = 0 + do i = 1, num + j = binom_rvs(n, p) + freq(j) = freq(j) + 1 + end do + chisq = 0.0_dp + do i = 0, n + expct = num * binom_pmf(i, n, p) + if(expct < 1.0e-5) cycle + chisq = chisq + (freq(i) - expct) ** 2 / expct + end do + write(*,*) "The critical values for chi-squared with 9 d. of f. are" + write(*,*) "1.15 and 27.87" + write(*,*) "Chi-squared for binomial random generator is : ", chisq + call check((chisq < 27.88 .and. chisq > 1.15), & + msg="binomial1 randomness failed chi-squared test", warn=warn) + n = 140 + p = 0.34 + deallocate(freq) + allocate(freq(0:n)) + freq = 0 + do i = 1, num + j = binom_rvs(n, p) + freq(j) = freq(j) + 1 + end do + chisq = 0.0_dp + do i = 0, n + expct = num * binom_pmf(i, n, p) + if(expct < 1.0e-5) cycle + chisq = chisq + (freq(i) - expct) ** 2 / expct + end do + write(*,*) "" + write(*,*) "The critical values for chi-squared with 49 d. of f. are" + write(*,*) "23.98 and 85.35" + write(*,*) "Chi-squared for binomial random generator is : ", chisq + call check((chisq < 85.35 .and. chisq > 23.98), & + msg="binomial1 randomness failed chi-squared test", warn=warn) + end subroutine test_binomial_random_generator + + subroutine test_uni_rvs_sp_int8 + real(sp) :: res(20), loc, scale + integer(int8) :: k + integer :: i, n, seed, get + + real(sp) :: ans(20) = & + [0.457413316_sp, 0.183665052_sp, 0.887956202_sp, & + 0.442960650_sp, 0.475367814_sp, 0.170218706_sp, & + 0.441669136_sp, 0.918695927_sp, 0.148022801_sp, & + 0.691296339_sp,-0.132472515_sp,-0.878723383_sp, & + -0.901660025_sp,-0.164090633_sp,-0.333886743_sp, & + 0.119948030_sp,-0.859304368_sp,-0.952883482_sp, & + 0.854362130_sp, 0.692578673_sp] + + print *, "Test uniform_distribution_rvs_sp_int8" + seed = 258147369 + call random_seed(seed, get) + loc = 0.0_sp + scale = 1.0_sp + k = 5_int8 + do i = 1, 5 + res(i) = uni_rvs(loc, scale) + end do + res(6:10) = uni_rvs(loc, scale, k) + loc = -1.0 + scale = 2.0 + do i = 11, 15 + res(i) = uni_rvs(loc, scale) + end do + res(16:20) = uni_rvs(loc, scale, k) + call check(all(abs(res - ans) < sptol), & + msg="uniform_distribution_rvs_sp_int8 failed", warn=warn) + end subroutine test_uni_rvs_sp_int8 + + subroutine test_uni_rvs_sp_int16 + real(sp) :: res(20), loc, scale + integer(int16) :: k + integer :: i, n, seed, get + + real(sp) :: ans(20) = & + [0.457413316_sp, 0.183665052_sp, 0.887956202_sp, & + 0.442960650_sp, 0.475367814_sp, 0.170218706_sp, & + 0.441669136_sp, 0.918695927_sp, 0.148022801_sp, & + 0.691296339_sp,-0.132472515_sp,-0.878723383_sp, & + -0.901660025_sp,-0.164090633_sp,-0.333886743_sp, & + 0.119948030_sp,-0.859304368_sp,-0.952883482_sp, & + 0.854362130_sp, 0.692578673_sp] + + print *, "Test uniform_distribution_rvs_sp_int16" + seed = 258147369 + call random_seed(seed, get) + loc = 0.0_sp + scale = 1.0_sp + k = 5_int16 + do i = 1, 5 + res(i) = uni_rvs(loc, scale) + end do + res(6:10) = uni_rvs(loc, scale, k) + loc = -1.0 + scale = 2.0 + do i = 11, 15 + res(i) = uni_rvs(loc, scale) + end do + res(16:20) = uni_rvs(loc, scale, k) + call check(all(abs(res - ans) < sptol), & + msg="uniform_distribution_rvs_sp_int16 failed", warn=warn) + end subroutine test_uni_rvs_sp_int16 + + subroutine test_uni_rvs_sp_int32 + real(sp) :: res(20), loc, scale + integer(int32) :: k + integer :: i, n, seed, get + + real(sp) :: ans(20) = & + [0.457413316_sp, 0.183665052_sp, 0.887956202_sp, & + 0.442960650_sp, 0.475367814_sp, 0.170218706_sp, & + 0.441669136_sp, 0.918695927_sp, 0.148022801_sp, & + 0.691296339_sp,-0.132472515_sp,-0.878723383_sp, & + -0.901660025_sp,-0.164090633_sp,-0.333886743_sp, & + 0.119948030_sp,-0.859304368_sp,-0.952883482_sp, & + 0.854362130_sp, 0.692578673_sp] + + print *, "Test uniform_distribution_rvs_sp_int32" + seed = 258147369 + call random_seed(seed, get) + loc = 0.0_sp + scale = 1.0_sp + k = 5_int32 + do i = 1, 5 + res(i) = uni_rvs(loc, scale) + end do + res(6:10) = uni_rvs(loc, scale, k) + loc = -1.0 + scale = 2.0 + do i = 11, 15 + res(i) = uni_rvs(loc, scale) + end do + res(16:20) = uni_rvs(loc, scale, k) + call check(all(abs(res - ans) < sptol), & + msg="uniform_distribution_rvs_sp_int32 failed", warn=warn) + end subroutine test_uni_rvs_sp_int32 + + subroutine test_uni_rvs_sp_int64 + real(sp) :: res(20), loc, scale + integer(int64) :: k + integer :: i, n, seed, get + + real(sp) :: ans(20) = & + [0.457413316_sp, 0.183665052_sp, 0.887956202_sp, & + 0.442960650_sp, 0.475367814_sp, 0.170218706_sp, & + 0.441669136_sp, 0.918695927_sp, 0.148022801_sp, & + 0.691296339_sp,-0.132472515_sp,-0.878723383_sp, & + -0.901660025_sp,-0.164090633_sp,-0.333886743_sp, & + 0.119948030_sp,-0.859304368_sp,-0.952883482_sp, & + 0.854362130_sp, 0.692578673_sp] + + print *, "Test uniform_distribution_rvs_sp_int64" + seed = 258147369 + call random_seed(seed, get) + loc = 0.0_sp + scale = 1.0_sp + k = 5_int64 + do i = 1, 5 + res(i) = uni_rvs(loc, scale) + end do + res(6:10) = uni_rvs(loc, scale, k) + loc = -1.0 + scale = 2.0 + do i = 11, 15 + res(i) = uni_rvs(loc, scale) + end do + res(16:20) = uni_rvs(loc, scale, k) + call check(all(abs(res - ans) < sptol), & + msg="uniform_distribution_rvs_sp_int64 failed", warn=warn) + end subroutine test_uni_rvs_sp_int64 + + subroutine test_uni_rvs_dp_int8 + real(dp) :: res(20), loc, scale + integer(int8) :: k + integer :: i, n, seed, get + + real(dp) :: ans(20) = & + [0.45741331426937459_dp, 0.18366504933248320_dp, & + 0.88795621528854640_dp, 0.44296065449379507_dp, & + 0.47536782827149393_dp, 0.17021871307147243_dp, & + 0.44166914074652641_dp, 0.91869594694692958_dp, & + 0.14802280170069973_dp, 0.69129636492557078_dp, & + -0.13247249397818517_dp,-0.87872336629421621_dp, & + -0.90166004614151585_dp,-0.16409061414773740_dp, & + -0.33388671819038429_dp, 0.11994803941701271_dp, & + -0.85930438039845658_dp,-0.95288345131793517_dp, & + 0.85436218295364763_dp, 0.69257869952149043_dp] + + + print *, "Test uniform_distribution_rvs_dp_int8" + seed = 258147369 + call random_seed(seed, get) + loc = 0.0_dp + scale = 1.0_dp + k = 5_int8 + do i = 1, 5 + res(i) = uni_rvs(loc, scale) + end do + res(6:10) = uni_rvs(loc, scale, k) + loc = -1.0 + scale = 2.0 + do i = 11, 15 + res(i) = uni_rvs(loc, scale) + end do + res(16:20) = uni_rvs(loc, scale, k) + call check(all(abs(res - ans) < dptol), & + msg="uniform_distribution_rvs_dp_int8 failed", warn=warn) + end subroutine test_uni_rvs_dp_int8 + + subroutine test_uni_rvs_dp_int16 + real(dp) :: res(20), loc, scale + integer(int16) :: k + integer :: i, n, seed, get + + real(dp) :: ans(20) = & + [0.45741331426937459_dp, 0.18366504933248320_dp, & + 0.88795621528854640_dp, 0.44296065449379507_dp, & + 0.47536782827149393_dp, 0.17021871307147243_dp, & + 0.44166914074652641_dp, 0.91869594694692958_dp, & + 0.14802280170069973_dp, 0.69129636492557078_dp, & + -0.13247249397818517_dp,-0.87872336629421621_dp, & + -0.90166004614151585_dp,-0.16409061414773740_dp, & + -0.33388671819038429_dp, 0.11994803941701271_dp, & + -0.85930438039845658_dp,-0.95288345131793517_dp, & + 0.85436218295364763_dp, 0.69257869952149043_dp] + + + print *, "Test uniform_distribution_rvs_dp_int16" + seed = 258147369 + call random_seed(seed, get) + loc = 0.0_dp + scale = 1.0_dp + k = 5_int16 + do i = 1, 5 + res(i) = uni_rvs(loc, scale) + end do + res(6:10) = uni_rvs(loc, scale, k) + loc = -1.0 + scale = 2.0 + do i = 11, 15 + res(i) = uni_rvs(loc, scale) + end do + res(16:20) = uni_rvs(loc, scale, k) + call check(all(abs(res - ans) < dptol), & + msg="uniform_distribution_rvs_dp_int16 failed", warn=warn) + end subroutine test_uni_rvs_dp_int16 + + subroutine test_uni_rvs_dp_int32 + real(dp) :: res(20), loc, scale + integer(int32) :: k + integer :: i, n, seed, get + + real(dp) :: ans(20) = & + [0.45741331426937459_dp, 0.18366504933248320_dp, & + 0.88795621528854640_dp, 0.44296065449379507_dp, & + 0.47536782827149393_dp, 0.17021871307147243_dp, & + 0.44166914074652641_dp, 0.91869594694692958_dp, & + 0.14802280170069973_dp, 0.69129636492557078_dp, & + -0.13247249397818517_dp,-0.87872336629421621_dp, & + -0.90166004614151585_dp,-0.16409061414773740_dp, & + -0.33388671819038429_dp, 0.11994803941701271_dp, & + -0.85930438039845658_dp,-0.95288345131793517_dp, & + 0.85436218295364763_dp, 0.69257869952149043_dp] + + + print *, "Test uniform_distribution_rvs_dp_int32" + seed = 258147369 + call random_seed(seed, get) + loc = 0.0_dp + scale = 1.0_dp + k = 5_int32 + do i = 1, 5 + res(i) = uni_rvs(loc, scale) + end do + res(6:10) = uni_rvs(loc, scale, k) + loc = -1.0 + scale = 2.0 + do i = 11, 15 + res(i) = uni_rvs(loc, scale) + end do + res(16:20) = uni_rvs(loc, scale, k) + call check(all(abs(res - ans) < dptol), & + msg="uniform_distribution_rvs_dp_int32 failed", warn=warn) + end subroutine test_uni_rvs_dp_int32 + + subroutine test_uni_rvs_dp_int64 + real(dp) :: res(20), loc, scale + integer(int64) :: k + integer :: i, n, seed, get + + real(dp) :: ans(20) = & + [0.45741331426937459_dp, 0.18366504933248320_dp, & + 0.88795621528854640_dp, 0.44296065449379507_dp, & + 0.47536782827149393_dp, 0.17021871307147243_dp, & + 0.44166914074652641_dp, 0.91869594694692958_dp, & + 0.14802280170069973_dp, 0.69129636492557078_dp, & + -0.13247249397818517_dp,-0.87872336629421621_dp, & + -0.90166004614151585_dp,-0.16409061414773740_dp, & + -0.33388671819038429_dp, 0.11994803941701271_dp, & + -0.85930438039845658_dp,-0.95288345131793517_dp, & + 0.85436218295364763_dp, 0.69257869952149043_dp] + + + print *, "Test uniform_distribution_rvs_dp_int64" + seed = 258147369 + call random_seed(seed, get) + loc = 0.0_dp + scale = 1.0_dp + k = 5_int64 + do i = 1, 5 + res(i) = uni_rvs(loc, scale) + end do + res(6:10) = uni_rvs(loc, scale, k) + loc = -1.0 + scale = 2.0 + do i = 11, 15 + res(i) = uni_rvs(loc, scale) + end do + res(16:20) = uni_rvs(loc, scale, k) + call check(all(abs(res - ans) < dptol), & + msg="uniform_distribution_rvs_dp_int64 failed", warn=warn) + end subroutine test_uni_rvs_dp_int64 + + subroutine test_uni_rvs_qp_int8 + real(qp) :: res(20), loc, scale + integer(int8) :: k + integer :: i, n, seed, get + + real(qp) :: ans(20) = & + [0.457413314269374593479255963757168502_qp, & + 0.183665049332483204524990583195176441_qp, & + 0.887956215288546402142344504682114348_qp, & + 0.442960654493795069619466175936395302_qp, & + 0.475367828271493930714086673106066883_qp, & + 0.170218713071472432796227280960010830_qp, & + 0.441669140746526411867023398372111842_qp, & + 0.918695946946929575815943280758801848_qp, & + 0.148022801700699729865462472844228614_qp, & + 0.691296364925570783199759716808330268_qp, & + -0.132472493978185168472805344208609313_qp, & + -0.878723366294216184924081858298450243_qp, & + -0.901660046141515819639877804547722917_qp, & + -0.164090614147737401395943379611708224_qp, & + -0.333886718190384290672056977200554684_qp, & + 0.119948039417012708440779533702880144_qp, & + -0.859304380398456552070385328079282772_qp, & + -0.952883451317935156743565983106236672_qp, & + 0.854362182953647630867521911568474025_qp, & + 0.692578699521490426249670235847588629_qp] + + + print *, "Test uniform_distribution_rvs_qp_int8" + seed = 258147369 + call random_seed(seed, get) + loc = 0.0_qp + scale = 1.0_qp + k = 5_int8 + do i = 1, 5 + res(i) = uni_rvs(loc, scale) + end do + res(6:10) = uni_rvs(loc, scale, k) + loc = -1.0 + scale = 2.0 + do i = 11, 15 + res(i) = uni_rvs(loc, scale) + end do + res(16:20) = uni_rvs(loc, scale, k) + call check(all(abs(res - ans) < qptol), & + msg="uniform_distribution_rvs_qp_int8 failed", warn=warn) + end subroutine test_uni_rvs_qp_int8 + + subroutine test_uni_rvs_qp_int16 + real(qp) :: res(20), loc, scale + integer(int16) :: k + integer :: i, n, seed, get + + real(qp) :: ans(20) = & + [0.457413314269374593479255963757168502_qp, & + 0.183665049332483204524990583195176441_qp, & + 0.887956215288546402142344504682114348_qp, & + 0.442960654493795069619466175936395302_qp, & + 0.475367828271493930714086673106066883_qp, & + 0.170218713071472432796227280960010830_qp, & + 0.441669140746526411867023398372111842_qp, & + 0.918695946946929575815943280758801848_qp, & + 0.148022801700699729865462472844228614_qp, & + 0.691296364925570783199759716808330268_qp, & + -0.132472493978185168472805344208609313_qp, & + -0.878723366294216184924081858298450243_qp, & + -0.901660046141515819639877804547722917_qp, & + -0.164090614147737401395943379611708224_qp, & + -0.333886718190384290672056977200554684_qp, & + 0.119948039417012708440779533702880144_qp, & + -0.859304380398456552070385328079282772_qp, & + -0.952883451317935156743565983106236672_qp, & + 0.854362182953647630867521911568474025_qp, & + 0.692578699521490426249670235847588629_qp] + + + print *, "Test uniform_distribution_rvs_qp_int16" + seed = 258147369 + call random_seed(seed, get) + loc = 0.0_qp + scale = 1.0_qp + k = 5_int16 + do i = 1, 5 + res(i) = uni_rvs(loc, scale) + end do + res(6:10) = uni_rvs(loc, scale, k) + loc = -1.0 + scale = 2.0 + do i = 11, 15 + res(i) = uni_rvs(loc, scale) + end do + res(16:20) = uni_rvs(loc, scale, k) + call check(all(abs(res - ans) < qptol), & + msg="uniform_distribution_rvs_qp_int16 failed", warn=warn) + end subroutine test_uni_rvs_qp_int16 + + subroutine test_uni_rvs_qp_int32 + real(qp) :: res(20), loc, scale + integer(int32) :: k + integer :: i, n, seed, get + + real(qp) :: ans(20) = & + [0.457413314269374593479255963757168502_qp, & + 0.183665049332483204524990583195176441_qp, & + 0.887956215288546402142344504682114348_qp, & + 0.442960654493795069619466175936395302_qp, & + 0.475367828271493930714086673106066883_qp, & + 0.170218713071472432796227280960010830_qp, & + 0.441669140746526411867023398372111842_qp, & + 0.918695946946929575815943280758801848_qp, & + 0.148022801700699729865462472844228614_qp, & + 0.691296364925570783199759716808330268_qp, & + -0.132472493978185168472805344208609313_qp, & + -0.878723366294216184924081858298450243_qp, & + -0.901660046141515819639877804547722917_qp, & + -0.164090614147737401395943379611708224_qp, & + -0.333886718190384290672056977200554684_qp, & + 0.119948039417012708440779533702880144_qp, & + -0.859304380398456552070385328079282772_qp, & + -0.952883451317935156743565983106236672_qp, & + 0.854362182953647630867521911568474025_qp, & + 0.692578699521490426249670235847588629_qp] + + + print *, "Test uniform_distribution_rvs_qp_int32" + seed = 258147369 + call random_seed(seed, get) + loc = 0.0_qp + scale = 1.0_qp + k = 5_int32 + do i = 1, 5 + res(i) = uni_rvs(loc, scale) + end do + res(6:10) = uni_rvs(loc, scale, k) + loc = -1.0 + scale = 2.0 + do i = 11, 15 + res(i) = uni_rvs(loc, scale) + end do + res(16:20) = uni_rvs(loc, scale, k) + call check(all(abs(res - ans) < qptol), & + msg="uniform_distribution_rvs_qp_int32 failed", warn=warn) + end subroutine test_uni_rvs_qp_int32 + + subroutine test_uni_rvs_qp_int64 + real(qp) :: res(20), loc, scale + integer(int64) :: k + integer :: i, n, seed, get + + real(qp) :: ans(20) = & + [0.457413314269374593479255963757168502_qp, & + 0.183665049332483204524990583195176441_qp, & + 0.887956215288546402142344504682114348_qp, & + 0.442960654493795069619466175936395302_qp, & + 0.475367828271493930714086673106066883_qp, & + 0.170218713071472432796227280960010830_qp, & + 0.441669140746526411867023398372111842_qp, & + 0.918695946946929575815943280758801848_qp, & + 0.148022801700699729865462472844228614_qp, & + 0.691296364925570783199759716808330268_qp, & + -0.132472493978185168472805344208609313_qp, & + -0.878723366294216184924081858298450243_qp, & + -0.901660046141515819639877804547722917_qp, & + -0.164090614147737401395943379611708224_qp, & + -0.333886718190384290672056977200554684_qp, & + 0.119948039417012708440779533702880144_qp, & + -0.859304380398456552070385328079282772_qp, & + -0.952883451317935156743565983106236672_qp, & + 0.854362182953647630867521911568474025_qp, & + 0.692578699521490426249670235847588629_qp] + + + print *, "Test uniform_distribution_rvs_qp_int64" + seed = 258147369 + call random_seed(seed, get) + loc = 0.0_qp + scale = 1.0_qp + k = 5_int64 + do i = 1, 5 + res(i) = uni_rvs(loc, scale) + end do + res(6:10) = uni_rvs(loc, scale, k) + loc = -1.0 + scale = 2.0 + do i = 11, 15 + res(i) = uni_rvs(loc, scale) + end do + res(16:20) = uni_rvs(loc, scale, k) + call check(all(abs(res - ans) < qptol), & + msg="uniform_distribution_rvs_qp_int64 failed", warn=warn) + end subroutine test_uni_rvs_qp_int64 + + + subroutine test_uni_pdf_sp + real(sp) :: x1, x2(3,4), res(3,5), loc, scale + integer :: i, n + integer :: seed, get + real(sp) :: ans(15) = [(0.2_sp, i=1,15)] + + + print *, "Test uniform_distribution_pdf_sp" + seed = 147258639 + call random_seed(seed, get) + loc = 0.0_sp + scale = 5.0_sp + x1 = uni_rvs(loc, scale) + x2 = reshape(uni_rvs(loc, scale, 12), [3,4]) + res(:,1) = uni_pdf(x1, loc, scale) + res(:, 2:5) = uni_pdf(x2, loc, scale) + call check(all(abs(res - reshape(ans,[3,5])) < sptol), & + msg="uniform_distribution_pdf_sp failed", warn=warn) + end subroutine test_uni_pdf_sp + + subroutine test_uni_pdf_dp + real(dp) :: x1, x2(3,4), res(3,5), loc, scale + integer :: i, n + integer :: seed, get + real(dp) :: ans(15) = [(0.2_dp, i=1,15)] + + + print *, "Test uniform_distribution_pdf_dp" + seed = 147258639 + call random_seed(seed, get) + loc = 0.0_dp + scale = 5.0_dp + x1 = uni_rvs(loc, scale) + x2 = reshape(uni_rvs(loc, scale, 12), [3,4]) + res(:,1) = uni_pdf(x1, loc, scale) + res(:, 2:5) = uni_pdf(x2, loc, scale) + call check(all(abs(res - reshape(ans,[3,5])) < dptol), & + msg="uniform_distribution_pdf_dp failed", warn=warn) + end subroutine test_uni_pdf_dp + + subroutine test_uni_pdf_qp + real(qp) :: x1, x2(3,4), res(3,5), loc, scale + integer :: i, n + integer :: seed, get + real(qp) :: ans(15) = [(0.2_qp, i=1,15)] + + + print *, "Test uniform_distribution_pdf_qp" + seed = 147258639 + call random_seed(seed, get) + loc = 0.0_qp + scale = 5.0_qp + x1 = uni_rvs(loc, scale) + x2 = reshape(uni_rvs(loc, scale, 12), [3,4]) + res(:,1) = uni_pdf(x1, loc, scale) + res(:, 2:5) = uni_pdf(x2, loc, scale) + call check(all(abs(res - reshape(ans,[3,5])) < qptol), & + msg="uniform_distribution_pdf_qp failed", warn=warn) + end subroutine test_uni_pdf_qp + + + subroutine test_uni_cdf_sp + real(sp) :: x1, x2(3,4), res(3,5), loc, scale + integer :: i, n + integer :: seed, get + real(sp) :: ans(15) = & + [0.170192957_sp, 0.170192957_sp, 0.170192957_sp, & + 0.276106149_sp, 0.754930079_sp, 0.406620681_sp, & + 0.187742829_sp, 0.651605546_sp, 0.934733927_sp, & + 0.151271492_sp, 0.987674534_sp, 0.130533904_sp, & + 0.106271923_sp, 9.27578807E-02_sp, 0.203399420_sp] + + print *, "Test uniform_distribution_cdf_sp" + seed = 369147258 + call random_seed(seed, get) + loc = -1.0_sp + scale = 2.0_sp + x1 = uni_rvs(loc, scale) + x2 = reshape(uni_rvs(loc, scale, 12), [3,4]) + res(:,1) = uni_cdf(x1, loc, scale) + res(:, 2:5) = uni_cdf(x2, loc, scale) + call check(all(abs(res - reshape(ans,[3,5])) < sptol), & + msg="uniform_distribution_cdf_sp failed", warn=warn) + end subroutine test_uni_cdf_sp + + subroutine test_uni_cdf_dp + real(dp) :: x1, x2(3,4), res(3,5), loc, scale + integer :: i, n + integer :: seed, get + real(dp) :: ans(15) = & + [0.17019294429755738_dp, 0.17019294429755738_dp, & + 0.17019294429755738_dp, 0.27610614627464619_dp, & + 0.75493009747387507_dp, 0.40662068257311801_dp, & + 0.18774281919180108_dp, 0.65160552609050759_dp, & + 0.93473394973210489_dp, 0.15127149185161326_dp, & + 0.98767452279771961_dp, 0.13053389946340466_dp, & + 0.10627190592187685_dp, 9.2757865224011304E-002_dp,& + 0.20339942685342044_dp] + + print *, "Test uniform_distribution_cdf_dp" + seed = 369147258 + call random_seed(seed, get) + loc = -1.0_dp + scale = 2.0_dp + x1 = uni_rvs(loc, scale) + x2 = reshape(uni_rvs(loc, scale, 12), [3,4]) + res(:,1) = uni_cdf(x1, loc, scale) + res(:, 2:5) = uni_cdf(x2, loc, scale) + call check(all(abs(res - reshape(ans,[3,5])) < dptol), & + msg="uniform_distribution_cdf_dp failed", warn=warn) + end subroutine test_uni_cdf_dp + + subroutine test_uni_cdf_qp + real(qp) :: x1, x2(3,4), res(3,5), loc, scale + integer :: i, n + integer :: seed, get + real(qp) :: ans(15) = & + [0.170192944297557408050991512027394492_qp, & + 0.170192944297557408050991512027394492_qp, & + 0.170192944297557408050991512027394492_qp, & + 0.276106146274646191418611351764411665_qp, & + 0.754930097473875072466853453079238534_qp, & + 0.406620682573118008562573777453508228_qp, & + 0.187742819191801080247472555129206739_qp, & + 0.651605526090507591874256831943057477_qp, & + 0.934733949732104885121941606485052034_qp, & + 0.151271491851613287815681019310432021_qp, & + 0.987674522797719611766353864368284121_qp, & + 0.130533899463404684526679488953959662_qp, & + 0.106271905921876880229959283497009892_qp, & + 9.27578652240113182836367400341259781E-0002_qp,& + 0.203399426853420439709196898547816090_qp] + + print *, "Test uniform_distribution_cdf_qp" + seed = 369147258 + call random_seed(seed, get) + loc = -1.0_qp + scale = 2.0_qp + x1 = uni_rvs(loc, scale) + x2 = reshape(uni_rvs(loc, scale, 12), [3,4]) + res(:,1) = uni_cdf(x1, loc, scale) + res(:, 2:5) = uni_cdf(x2, loc, scale) + call check(all(abs(res - reshape(ans,[3,5])) < qptol), & + msg="uniform_distribution_cdf_qp failed", warn=warn) + end subroutine test_uni_cdf_qp + + + subroutine test_nor_rvs_sp_int8 + real(sp) :: res(20), loc, scale + integer(int8) :: k + integer :: i, n + integer :: seed, get + real(sp) :: ans(20) = & + [1.08354020_sp, 0.930153966_sp, 0.388561100_sp, & + -1.44566274_sp, 0.612832963_sp, -1.00310886_sp, & + 0.817594171_sp, -0.568394125_sp, 0.993003964_sp, & + -1.76855731_sp, -3.22572017_sp, -3.99702096_sp, & + 0.428807259_sp, -2.44686961_sp, -3.48141479_sp, & + -4.33496284_sp, -0.154625356_sp,-0.830695271_sp, & + -3.90960717_sp, 1.53445792_sp] + + print *, "Test normal_distribution_rvs_sp_int8" + seed = 25836914 + call random_seed(seed, get) + loc = 0.0_sp + scale = 1.0_sp + k = 5_int8 + do i = 1, 5 + res(i) = nor_rvs(loc, scale) + end do + res(6:10) = nor_rvs(loc, scale, k) + loc = -1.0 + scale = 2.0 + do i = 11, 15 + res(i) = nor_rvs(loc, scale) + end do + res(16:20) = nor_rvs(loc, scale, k) + call check(all(abs(res - ans) < sptol), & + msg="normal_distribution_rvs_sp_int8 failed", warn=warn) + end subroutine test_nor_rvs_sp_int8 + + subroutine test_nor_rvs_sp_int16 + real(sp) :: res(20), loc, scale + integer(int16) :: k + integer :: i, n + integer :: seed, get + real(sp) :: ans(20) = & + [1.08354020_sp, 0.930153966_sp, 0.388561100_sp, & + -1.44566274_sp, 0.612832963_sp, -1.00310886_sp, & + 0.817594171_sp, -0.568394125_sp, 0.993003964_sp, & + -1.76855731_sp, -3.22572017_sp, -3.99702096_sp, & + 0.428807259_sp, -2.44686961_sp, -3.48141479_sp, & + -4.33496284_sp, -0.154625356_sp,-0.830695271_sp, & + -3.90960717_sp, 1.53445792_sp] + + print *, "Test normal_distribution_rvs_sp_int16" + seed = 25836914 + call random_seed(seed, get) + loc = 0.0_sp + scale = 1.0_sp + k = 5_int16 + do i = 1, 5 + res(i) = nor_rvs(loc, scale) + end do + res(6:10) = nor_rvs(loc, scale, k) + loc = -1.0 + scale = 2.0 + do i = 11, 15 + res(i) = nor_rvs(loc, scale) + end do + res(16:20) = nor_rvs(loc, scale, k) + call check(all(abs(res - ans) < sptol), & + msg="normal_distribution_rvs_sp_int16 failed", warn=warn) + end subroutine test_nor_rvs_sp_int16 + + subroutine test_nor_rvs_sp_int32 + real(sp) :: res(20), loc, scale + integer(int32) :: k + integer :: i, n + integer :: seed, get + real(sp) :: ans(20) = & + [1.08354020_sp, 0.930153966_sp, 0.388561100_sp, & + -1.44566274_sp, 0.612832963_sp, -1.00310886_sp, & + 0.817594171_sp, -0.568394125_sp, 0.993003964_sp, & + -1.76855731_sp, -3.22572017_sp, -3.99702096_sp, & + 0.428807259_sp, -2.44686961_sp, -3.48141479_sp, & + -4.33496284_sp, -0.154625356_sp,-0.830695271_sp, & + -3.90960717_sp, 1.53445792_sp] + + print *, "Test normal_distribution_rvs_sp_int32" + seed = 25836914 + call random_seed(seed, get) + loc = 0.0_sp + scale = 1.0_sp + k = 5_int32 + do i = 1, 5 + res(i) = nor_rvs(loc, scale) + end do + res(6:10) = nor_rvs(loc, scale, k) + loc = -1.0 + scale = 2.0 + do i = 11, 15 + res(i) = nor_rvs(loc, scale) + end do + res(16:20) = nor_rvs(loc, scale, k) + call check(all(abs(res - ans) < sptol), & + msg="normal_distribution_rvs_sp_int32 failed", warn=warn) + end subroutine test_nor_rvs_sp_int32 + + subroutine test_nor_rvs_sp_int64 + real(sp) :: res(20), loc, scale + integer(int64) :: k + integer :: i, n + integer :: seed, get + real(sp) :: ans(20) = & + [1.08354020_sp, 0.930153966_sp, 0.388561100_sp, & + -1.44566274_sp, 0.612832963_sp, -1.00310886_sp, & + 0.817594171_sp, -0.568394125_sp, 0.993003964_sp, & + -1.76855731_sp, -3.22572017_sp, -3.99702096_sp, & + 0.428807259_sp, -2.44686961_sp, -3.48141479_sp, & + -4.33496284_sp, -0.154625356_sp,-0.830695271_sp, & + -3.90960717_sp, 1.53445792_sp] + + print *, "Test normal_distribution_rvs_sp_int64" + seed = 25836914 + call random_seed(seed, get) + loc = 0.0_sp + scale = 1.0_sp + k = 5_int64 + do i = 1, 5 + res(i) = nor_rvs(loc, scale) + end do + res(6:10) = nor_rvs(loc, scale, k) + loc = -1.0 + scale = 2.0 + do i = 11, 15 + res(i) = nor_rvs(loc, scale) + end do + res(16:20) = nor_rvs(loc, scale, k) + call check(all(abs(res - ans) < sptol), & + msg="normal_distribution_rvs_sp_int64 failed", warn=warn) + end subroutine test_nor_rvs_sp_int64 + + subroutine test_nor_rvs_dp_int8 + real(dp) :: res(20), loc, scale + integer(int8) :: k + integer :: i, n + integer :: seed, get + real(dp) :: ans(20) = & + [1.0835401965902034_dp, 0.93015397468064165_dp, & + 0.38856109396542121_dp, -1.4456627206540740_dp, & + 0.61283297553014326_dp, -1.0031088776888382_dp, & + 0.81759417579176041_dp, -0.56839412687107116_dp, & + 0.99300393889422900_dp, -1.7685573691899061_dp, & + -3.2257201976639149_dp, -3.9970210500520191_dp, & + 0.42880723775826013_dp, -2.4468696145147510_dp, & + -3.4814147687231882_dp, -4.3349631940225235_dp, & + -0.15462537592816106_dp, -0.83069527405234211_dp, & + -3.9096071456063441_dp, 1.5344578674016103_dp] + + print *, "Test normal_distribution_rvs_dp_int8" + seed = 25836914 + call random_seed(seed, get) + loc = 0.0_dp + scale = 1.0_dp + k = 5_int8 + do i = 1, 5 + res(i) = nor_rvs(loc, scale) + end do + res(6:10) = nor_rvs(loc, scale, k) + loc = -1.0 + scale = 2.0 + do i = 11, 15 + res(i) = nor_rvs(loc, scale) + end do + res(16:20) = nor_rvs(loc, scale, k) + call check(all(abs(res - ans) < dptol), & + msg="normal_distribution_rvs_dp_int8 failed", warn=warn) + end subroutine test_nor_rvs_dp_int8 + + subroutine test_nor_rvs_dp_int16 + real(dp) :: res(20), loc, scale + integer(int16) :: k + integer :: i, n + integer :: seed, get + real(dp) :: ans(20) = & + [1.0835401965902034_dp, 0.93015397468064165_dp, & + 0.38856109396542121_dp, -1.4456627206540740_dp, & + 0.61283297553014326_dp, -1.0031088776888382_dp, & + 0.81759417579176041_dp, -0.56839412687107116_dp, & + 0.99300393889422900_dp, -1.7685573691899061_dp, & + -3.2257201976639149_dp, -3.9970210500520191_dp, & + 0.42880723775826013_dp, -2.4468696145147510_dp, & + -3.4814147687231882_dp, -4.3349631940225235_dp, & + -0.15462537592816106_dp, -0.83069527405234211_dp, & + -3.9096071456063441_dp, 1.5344578674016103_dp] + + print *, "Test normal_distribution_rvs_dp_int16" + seed = 25836914 + call random_seed(seed, get) + loc = 0.0_dp + scale = 1.0_dp + k = 5_int16 + do i = 1, 5 + res(i) = nor_rvs(loc, scale) + end do + res(6:10) = nor_rvs(loc, scale, k) + loc = -1.0 + scale = 2.0 + do i = 11, 15 + res(i) = nor_rvs(loc, scale) + end do + res(16:20) = nor_rvs(loc, scale, k) + call check(all(abs(res - ans) < dptol), & + msg="normal_distribution_rvs_dp_int16 failed", warn=warn) + end subroutine test_nor_rvs_dp_int16 + + subroutine test_nor_rvs_dp_int32 + real(dp) :: res(20), loc, scale + integer(int32) :: k + integer :: i, n + integer :: seed, get + real(dp) :: ans(20) = & + [1.0835401965902034_dp, 0.93015397468064165_dp, & + 0.38856109396542121_dp, -1.4456627206540740_dp, & + 0.61283297553014326_dp, -1.0031088776888382_dp, & + 0.81759417579176041_dp, -0.56839412687107116_dp, & + 0.99300393889422900_dp, -1.7685573691899061_dp, & + -3.2257201976639149_dp, -3.9970210500520191_dp, & + 0.42880723775826013_dp, -2.4468696145147510_dp, & + -3.4814147687231882_dp, -4.3349631940225235_dp, & + -0.15462537592816106_dp, -0.83069527405234211_dp, & + -3.9096071456063441_dp, 1.5344578674016103_dp] + + print *, "Test normal_distribution_rvs_dp_int32" + seed = 25836914 + call random_seed(seed, get) + loc = 0.0_dp + scale = 1.0_dp + k = 5_int32 + do i = 1, 5 + res(i) = nor_rvs(loc, scale) + end do + res(6:10) = nor_rvs(loc, scale, k) + loc = -1.0 + scale = 2.0 + do i = 11, 15 + res(i) = nor_rvs(loc, scale) + end do + res(16:20) = nor_rvs(loc, scale, k) + call check(all(abs(res - ans) < dptol), & + msg="normal_distribution_rvs_dp_int32 failed", warn=warn) + end subroutine test_nor_rvs_dp_int32 + + subroutine test_nor_rvs_dp_int64 + real(dp) :: res(20), loc, scale + integer(int64) :: k + integer :: i, n + integer :: seed, get + real(dp) :: ans(20) = & + [1.0835401965902034_dp, 0.93015397468064165_dp, & + 0.38856109396542121_dp, -1.4456627206540740_dp, & + 0.61283297553014326_dp, -1.0031088776888382_dp, & + 0.81759417579176041_dp, -0.56839412687107116_dp, & + 0.99300393889422900_dp, -1.7685573691899061_dp, & + -3.2257201976639149_dp, -3.9970210500520191_dp, & + 0.42880723775826013_dp, -2.4468696145147510_dp, & + -3.4814147687231882_dp, -4.3349631940225235_dp, & + -0.15462537592816106_dp, -0.83069527405234211_dp, & + -3.9096071456063441_dp, 1.5344578674016103_dp] + + print *, "Test normal_distribution_rvs_dp_int64" + seed = 25836914 + call random_seed(seed, get) + loc = 0.0_dp + scale = 1.0_dp + k = 5_int64 + do i = 1, 5 + res(i) = nor_rvs(loc, scale) + end do + res(6:10) = nor_rvs(loc, scale, k) + loc = -1.0 + scale = 2.0 + do i = 11, 15 + res(i) = nor_rvs(loc, scale) + end do + res(16:20) = nor_rvs(loc, scale, k) + call check(all(abs(res - ans) < dptol), & + msg="normal_distribution_rvs_dp_int64 failed", warn=warn) + end subroutine test_nor_rvs_dp_int64 + + subroutine test_nor_rvs_qp_int8 + real(qp) :: res(20), loc, scale + integer(int8) :: k + integer :: i, n + integer :: seed, get + real(qp) :: ans(20) = & + [1.08354019659020339716448688704986125_qp, & + 0.930153974680641648653534048207802698_qp, & + 0.388561093965421211482436092410353012_qp, & + -1.44566272065407397384717569366330281_qp, & + 0.612832975530143264641935729741817340_qp, & + -1.00310887768883816306697553955018520_qp, & + 0.817594175791760413574138510739430785_qp, & + -0.568394126871071159179393816884839907_qp, & + 0.993003938894228999068047869513975456_qp, & + -1.76855736918990613659730115614365786_qp, & + -3.22572019766391493433843606908340007_qp, & + -3.99702105005201913101586796983610839_qp, & + 0.428807237758260129112386493943631649_qp, & + -2.44686961451475104567521157150622457_qp, & + -3.48141476872318822444185570930130780_qp, & + -4.33496319402252350272419789689593017_qp, & + -0.154625375928161057359488950169179589_qp, & + -0.830695274052342114146085805259644985_qp, & + -3.90960714560634414738160558044910431_qp, & + 1.53445786740161027594808729190845042_qp] + + print *, "Test normal_distribution_rvs_qp_int8" + seed = 25836914 + call random_seed(seed, get) + loc = 0.0_qp + scale = 1.0_qp + k = 5_int8 + do i = 1, 5 + res(i) = nor_rvs(loc, scale) + end do + res(6:10) = nor_rvs(loc, scale, k) + loc = -1.0 + scale = 2.0 + do i = 11, 15 + res(i) = nor_rvs(loc, scale) + end do + res(16:20) = nor_rvs(loc, scale, k) + call check(all(abs(res - ans) < qptol), & + msg="normal_distribution_rvs_qp_int8 failed", warn=warn) + end subroutine test_nor_rvs_qp_int8 + + subroutine test_nor_rvs_qp_int16 + real(qp) :: res(20), loc, scale + integer(int16) :: k + integer :: i, n + integer :: seed, get + real(qp) :: ans(20) = & + [1.08354019659020339716448688704986125_qp, & + 0.930153974680641648653534048207802698_qp, & + 0.388561093965421211482436092410353012_qp, & + -1.44566272065407397384717569366330281_qp, & + 0.612832975530143264641935729741817340_qp, & + -1.00310887768883816306697553955018520_qp, & + 0.817594175791760413574138510739430785_qp, & + -0.568394126871071159179393816884839907_qp, & + 0.993003938894228999068047869513975456_qp, & + -1.76855736918990613659730115614365786_qp, & + -3.22572019766391493433843606908340007_qp, & + -3.99702105005201913101586796983610839_qp, & + 0.428807237758260129112386493943631649_qp, & + -2.44686961451475104567521157150622457_qp, & + -3.48141476872318822444185570930130780_qp, & + -4.33496319402252350272419789689593017_qp, & + -0.154625375928161057359488950169179589_qp, & + -0.830695274052342114146085805259644985_qp, & + -3.90960714560634414738160558044910431_qp, & + 1.53445786740161027594808729190845042_qp] + + print *, "Test normal_distribution_rvs_qp_int16" + seed = 25836914 + call random_seed(seed, get) + loc = 0.0_qp + scale = 1.0_qp + k = 5_int16 + do i = 1, 5 + res(i) = nor_rvs(loc, scale) + end do + res(6:10) = nor_rvs(loc, scale, k) + loc = -1.0 + scale = 2.0 + do i = 11, 15 + res(i) = nor_rvs(loc, scale) + end do + res(16:20) = nor_rvs(loc, scale, k) + call check(all(abs(res - ans) < qptol), & + msg="normal_distribution_rvs_qp_int16 failed", warn=warn) + end subroutine test_nor_rvs_qp_int16 + + subroutine test_nor_rvs_qp_int32 + real(qp) :: res(20), loc, scale + integer(int32) :: k + integer :: i, n + integer :: seed, get + real(qp) :: ans(20) = & + [1.08354019659020339716448688704986125_qp, & + 0.930153974680641648653534048207802698_qp, & + 0.388561093965421211482436092410353012_qp, & + -1.44566272065407397384717569366330281_qp, & + 0.612832975530143264641935729741817340_qp, & + -1.00310887768883816306697553955018520_qp, & + 0.817594175791760413574138510739430785_qp, & + -0.568394126871071159179393816884839907_qp, & + 0.993003938894228999068047869513975456_qp, & + -1.76855736918990613659730115614365786_qp, & + -3.22572019766391493433843606908340007_qp, & + -3.99702105005201913101586796983610839_qp, & + 0.428807237758260129112386493943631649_qp, & + -2.44686961451475104567521157150622457_qp, & + -3.48141476872318822444185570930130780_qp, & + -4.33496319402252350272419789689593017_qp, & + -0.154625375928161057359488950169179589_qp, & + -0.830695274052342114146085805259644985_qp, & + -3.90960714560634414738160558044910431_qp, & + 1.53445786740161027594808729190845042_qp] + + print *, "Test normal_distribution_rvs_qp_int32" + seed = 25836914 + call random_seed(seed, get) + loc = 0.0_qp + scale = 1.0_qp + k = 5_int32 + do i = 1, 5 + res(i) = nor_rvs(loc, scale) + end do + res(6:10) = nor_rvs(loc, scale, k) + loc = -1.0 + scale = 2.0 + do i = 11, 15 + res(i) = nor_rvs(loc, scale) + end do + res(16:20) = nor_rvs(loc, scale, k) + call check(all(abs(res - ans) < qptol), & + msg="normal_distribution_rvs_qp_int32 failed", warn=warn) + end subroutine test_nor_rvs_qp_int32 + + subroutine test_nor_rvs_qp_int64 + real(qp) :: res(20), loc, scale + integer(int64) :: k + integer :: i, n + integer :: seed, get + real(qp) :: ans(20) = & + [1.08354019659020339716448688704986125_qp, & + 0.930153974680641648653534048207802698_qp, & + 0.388561093965421211482436092410353012_qp, & + -1.44566272065407397384717569366330281_qp, & + 0.612832975530143264641935729741817340_qp, & + -1.00310887768883816306697553955018520_qp, & + 0.817594175791760413574138510739430785_qp, & + -0.568394126871071159179393816884839907_qp, & + 0.993003938894228999068047869513975456_qp, & + -1.76855736918990613659730115614365786_qp, & + -3.22572019766391493433843606908340007_qp, & + -3.99702105005201913101586796983610839_qp, & + 0.428807237758260129112386493943631649_qp, & + -2.44686961451475104567521157150622457_qp, & + -3.48141476872318822444185570930130780_qp, & + -4.33496319402252350272419789689593017_qp, & + -0.154625375928161057359488950169179589_qp, & + -0.830695274052342114146085805259644985_qp, & + -3.90960714560634414738160558044910431_qp, & + 1.53445786740161027594808729190845042_qp] + + print *, "Test normal_distribution_rvs_qp_int64" + seed = 25836914 + call random_seed(seed, get) + loc = 0.0_qp + scale = 1.0_qp + k = 5_int64 + do i = 1, 5 + res(i) = nor_rvs(loc, scale) + end do + res(6:10) = nor_rvs(loc, scale, k) + loc = -1.0 + scale = 2.0 + do i = 11, 15 + res(i) = nor_rvs(loc, scale) + end do + res(16:20) = nor_rvs(loc, scale, k) + call check(all(abs(res - ans) < qptol), & + msg="normal_distribution_rvs_qp_int64 failed", warn=warn) + end subroutine test_nor_rvs_qp_int64 + + + + subroutine test_nor_pdf_sp + real(sp) :: x1, x2(3,4), res(3,5), loc, scale + integer :: i, n + integer :: seed, get + real(sp) :: ans(15) = & + [0.322576135_sp, 0.322576135_sp, 0.322576135_sp, & + 0.300806433_sp, 8.49242210E-02_sp, 0.358480453_sp, & + 0.398903936_sp, 0.393221349_sp, 0.374799609_sp, & + 5.98081723E-02_sp, 0.398853570_sp, 0.241967395_sp, & + 0.373766601_sp, 0.356140822_sp, 0.233544141_sp] + + print *, "Test normal_distribution_pdf_sp" + seed = 741852963 + call random_seed(seed, get) + loc = 0.0_sp + scale = 1.0_sp + 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])) < sptol), & + msg="normal_distribution_pdf_sp failed", warn=warn) + end subroutine test_nor_pdf_sp + + subroutine test_nor_pdf_dp + real(dp) :: x1, x2(3,4), res(3,5), loc, scale + integer :: i, n + integer :: seed, get + real(dp) :: ans(15) = & + [0.32257615048492366_dp, 0.32257615048492366_dp, & + 0.32257615048492366_dp, 0.30080644003932094_dp, & + 8.4924229110490274E-002_dp, 0.35848043641803234_dp, & + 0.39890395411791441_dp, 0.39322133798111997_dp, & + 0.37479961337242840_dp, 5.9808167624805800E-002_dp, & + 0.39885355470530021_dp, 0.24196740475597517_dp, & + 0.37376661053141419_dp, 0.35614082586331985_dp, & + 0.23354412957618306_dp] + + print *, "Test normal_distribution_pdf_dp" + seed = 741852963 + call random_seed(seed, get) + loc = 0.0_dp + scale = 1.0_dp + 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])) < dptol), & + msg="normal_distribution_pdf_dp failed", warn=warn) + end subroutine test_nor_pdf_dp + + subroutine test_nor_pdf_qp + real(qp) :: x1, x2(3,4), res(3,5), loc, scale + integer :: i, n + integer :: seed, get + real(qp) :: ans(15) = & + [0.322576150484923624816177827114417878_qp, & + 0.322576150484923624816177827114417878_qp, & + 0.322576150484923624816177827114417878_qp, & + 0.300806440039320895258949727681658841_qp, & + 8.49242291104902651552862536156033805E-0002_qp, & + 0.358480436418032272301373539277868731_qp, & + 0.398903954117914366544457439174698638_qp, & + 0.393221337981119941663547835562354415_qp, & + 0.374799613372428368299981031485342040_qp, & + 5.98081676248057976816689135970544796E-0002_qp, & + 0.398853554705300200346860042267642541_qp, & + 0.241967404755975138058416435199686964_qp, & + 0.373766610531414167998075638655556567_qp, & + 0.356140825863319809711905710918457326_qp, & + 0.233544129576183026277279390942135717_qp] + + print *, "Test normal_distribution_pdf_qp" + seed = 741852963 + call random_seed(seed, get) + loc = 0.0_qp + scale = 1.0_qp + 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])) < qptol), & + msg="normal_distribution_pdf_qp failed", warn=warn) + end subroutine test_nor_pdf_qp + + + subroutine test_nor_cdf_sp + real(sp) :: x1, x2(3,4), res(3,5), loc, scale + integer :: i, n + integer :: seed, get + real(sp) :: ans(15) = & + [7.50826299E-02_sp, 7.50826299E-02_sp, 7.50826299E-02_sp, & + 0.143119842_sp, 0.241425395_sp, 0.284345925_sp, & + 0.233239830_sp, 0.341059506_sp, 0.353156865_sp, & + 6.81066811E-02_sp, 4.38792408E-02_sp, 0.763679624_sp, & + 0.363722205_sp, 0.868187129_sp, 0.626506805_sp] + + print *, "Test normal_distribution_cdf_sp" + seed = 369147582 + call random_seed(seed, get) + loc = -1.0_sp + scale = 2.0_sp + 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])) < sptol), & + msg="normal_distribution_cdf_sp failed", warn=warn) + end subroutine test_nor_cdf_sp + + subroutine test_nor_cdf_dp + real(dp) :: x1, x2(3,4), res(3,5), loc, scale + integer :: i, n + integer :: seed, get + real(dp) :: ans(15) = & + [7.5082630503844117E-002_dp, 7.5082630503844117E-002_dp, & + 7.5082630503844117E-002_dp, 0.14311983410871804_dp, & + 0.24142542152570318_dp, 0.28434587862603933_dp, & + 0.23323983636601592_dp, 0.34105950613721914_dp, & + 0.35315685019983512_dp, 6.8106676639663855E-002_dp, & + 4.3879233144168306E-002_dp, 0.76367963788286075_dp, & + 0.36372218758735508_dp, 0.86818711488498046_dp, & + 0.62650679980965285_dp] + + print *, "Test normal_distribution_cdf_dp" + seed = 369147582 + call random_seed(seed, get) + loc = -1.0_dp + scale = 2.0_dp + 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])) < dptol), & + msg="normal_distribution_cdf_dp failed", warn=warn) + end subroutine test_nor_cdf_dp + + subroutine test_nor_cdf_qp + real(qp) :: x1, x2(3,4), res(3,5), loc, scale + integer :: i, n + integer :: seed, get + real(qp) :: ans(15) = & + [7.50826305038441048487991102776954069E-0002_qp, & + 7.50826305038441048487991102776954069E-0002_qp, & + 7.50826305038441048487991102776954069E-0002_qp, & + 0.143119834108717983250834016885129863_qp, & + 0.241425421525703182028420560765471735_qp, & + 0.284345878626039240974266199229875972_qp, & + 0.233239836366015928845367994433532733_qp, & + 0.341059506137219171082517155967522896_qp, & + 0.353156850199835111081038166086606192_qp, & + 6.81066766396638231790017005897813364E-0002_qp, & + 4.38792331441682923984716366123285768E-0002_qp, & + 0.763679637882860826030745070304416929_qp, & + 0.363722187587355040667876190724308059_qp, & + 0.868187114884980488672309198087692444_qp, & + 0.626506799809652872401992867475200722_qp] + + print *, "Test normal_distribution_cdf_qp" + seed = 369147582 + call random_seed(seed, get) + loc = -1.0_qp + scale = 2.0_qp + 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])) < qptol), & + msg="normal_distribution_cdf_qp failed", warn=warn) + end subroutine test_nor_cdf_qp + + + subroutine test_binom_rvs_int8 + integer(int8) :: res(40), k, n + integer :: i, m + real :: p + integer :: seed, get + integer(int8) :: ans(40) = & + [71_int8, 81_int8, 75_int8, 78_int8, 74_int8, & + 77_int8, 81_int8, 80_int8, 81_int8, 78_int8, & + 10_int8, 5_int8, 6_int8, 12_int8, 10_int8, & + 10_int8, 10_int8, 11_int8, 11_int8, 9_int8, & + 28_int8, 27_int8, 26_int8, 28_int8, 19_int8, & + 20_int8, 29_int8, 25_int8, 23_int8, 27_int8, & + 6_int8, 4_int8, 5_int8, 5_int8, 7_int8, & + 8_int8, 6_int8, 7_int8, 8_int8, 7_int8] + + print *, "Test binomial_distribution_rvs_int8" + seed = 852693417 + call random_seed(seed, get) + n = 100_int8 + p = 0.76 + k = 5_int8 + do i = 1, 5 + res(i) = binom_rvs(n, p) + end do + res(6:10) = binom_rvs(n, p, k) + n = 20_int8 + p = 0.4 + do i = 11, 15 + res(i) = binom_rvs(n, P) + end do + res(16:20) = binom_rvs(n, p, k) + + n = 100_int8 + p = 0.26 + k = 5_int8 + do i = 21, 25 + res(i) = binom_rvs(n, p) + end do + res(26:30) = binom_rvs(n, p, k) + n = 10_int8 + p = 0.7 + do i = 31, 35 + res(i) = binom_rvs(n, P) + end do + res(36:40) = binom_rvs(n, p, k) + call check(all(res == ans), & + msg="binomial_distribution_rvs_int8_int64 failed", warn=warn) + end subroutine test_binom_rvs_int8 + + subroutine test_binom_rvs_int16 + integer(int16) :: res(40), k, n + integer :: i, m + real :: p + integer :: seed, get + integer(int16) :: ans(40) = & + [71_int16, 81_int16, 75_int16, 78_int16, 74_int16, & + 77_int16, 81_int16, 80_int16, 81_int16, 78_int16, & + 10_int16, 5_int16, 6_int16, 12_int16, 10_int16, & + 10_int16, 10_int16, 11_int16, 11_int16, 9_int16, & + 28_int16, 27_int16, 26_int16, 28_int16, 19_int16, & + 20_int16, 29_int16, 25_int16, 23_int16, 27_int16, & + 6_int16, 4_int16, 5_int16, 5_int16, 7_int16, & + 8_int16, 6_int16, 7_int16, 8_int16, 7_int16] + + print *, "Test binomial_distribution_rvs_int16" + seed = 852693417 + call random_seed(seed, get) + n = 100_int16 + p = 0.76 + k = 5_int16 + do i = 1, 5 + res(i) = binom_rvs(n, p) + end do + res(6:10) = binom_rvs(n, p, k) + n = 20_int16 + p = 0.4 + do i = 11, 15 + res(i) = binom_rvs(n, P) + end do + res(16:20) = binom_rvs(n, p, k) + + n = 100_int16 + p = 0.26 + k = 5_int16 + do i = 21, 25 + res(i) = binom_rvs(n, p) + end do + res(26:30) = binom_rvs(n, p, k) + n = 10_int16 + p = 0.7 + do i = 31, 35 + res(i) = binom_rvs(n, P) + end do + res(36:40) = binom_rvs(n, p, k) + call check(all(res == ans), & + msg="binomial_distribution_rvs_int16_int64 failed", warn=warn) + end subroutine test_binom_rvs_int16 + + subroutine test_binom_rvs_int32 + integer(int32) :: res(40), k, n + integer :: i, m + real :: p + integer :: seed, get + integer(int32) :: ans(40) = & + [71_int32, 81_int32, 75_int32, 78_int32, 74_int32, & + 77_int32, 81_int32, 80_int32, 81_int32, 78_int32, & + 10_int32, 5_int32, 6_int32, 12_int32, 10_int32, & + 10_int32, 10_int32, 11_int32, 11_int32, 9_int32, & + 28_int32, 27_int32, 26_int32, 28_int32, 19_int32, & + 20_int32, 29_int32, 25_int32, 23_int32, 27_int32, & + 6_int32, 4_int32, 5_int32, 5_int32, 7_int32, & + 8_int32, 6_int32, 7_int32, 8_int32, 7_int32] + + print *, "Test binomial_distribution_rvs_int32" + seed = 852693417 + call random_seed(seed, get) + n = 100_int32 + p = 0.76 + k = 5_int32 + do i = 1, 5 + res(i) = binom_rvs(n, p) + end do + res(6:10) = binom_rvs(n, p, k) + n = 20_int32 + p = 0.4 + do i = 11, 15 + res(i) = binom_rvs(n, P) + end do + res(16:20) = binom_rvs(n, p, k) + + n = 100_int32 + p = 0.26 + k = 5_int32 + do i = 21, 25 + res(i) = binom_rvs(n, p) + end do + res(26:30) = binom_rvs(n, p, k) + n = 10_int32 + p = 0.7 + do i = 31, 35 + res(i) = binom_rvs(n, P) + end do + res(36:40) = binom_rvs(n, p, k) + call check(all(res == ans), & + msg="binomial_distribution_rvs_int32_int64 failed", warn=warn) + end subroutine test_binom_rvs_int32 + + subroutine test_binom_rvs_int64 + integer(int64) :: res(40), k, n + integer :: i, m + real :: p + integer :: seed, get + integer(int64) :: ans(40) = & + [71_int64, 81_int64, 75_int64, 78_int64, 74_int64, & + 77_int64, 81_int64, 80_int64, 81_int64, 78_int64, & + 10_int64, 5_int64, 6_int64, 12_int64, 10_int64, & + 10_int64, 10_int64, 11_int64, 11_int64, 9_int64, & + 28_int64, 27_int64, 26_int64, 28_int64, 19_int64, & + 20_int64, 29_int64, 25_int64, 23_int64, 27_int64, & + 6_int64, 4_int64, 5_int64, 5_int64, 7_int64, & + 8_int64, 6_int64, 7_int64, 8_int64, 7_int64] + + print *, "Test binomial_distribution_rvs_int64" + seed = 852693417 + call random_seed(seed, get) + n = 100_int64 + p = 0.76 + k = 5_int64 + do i = 1, 5 + res(i) = binom_rvs(n, p) + end do + res(6:10) = binom_rvs(n, p, k) + n = 20_int64 + p = 0.4 + do i = 11, 15 + res(i) = binom_rvs(n, P) + end do + res(16:20) = binom_rvs(n, p, k) + + n = 100_int64 + p = 0.26 + k = 5_int64 + do i = 21, 25 + res(i) = binom_rvs(n, p) + end do + res(26:30) = binom_rvs(n, p, k) + n = 10_int64 + p = 0.7 + do i = 31, 35 + res(i) = binom_rvs(n, P) + end do + res(36:40) = binom_rvs(n, p, k) + call check(all(res == ans), & + msg="binomial_distribution_rvs_int64_int64 failed", warn=warn) + end subroutine test_binom_rvs_int64 + + + subroutine test_binom_pmf_int8 + integer(int8) :: x1, x2(3,4), n + integer :: i, m + integer :: seed, get + real :: p, res(3,5) + real :: ans(15) = [7.78146312E-02, 7.78146312E-02, 7.78146312E-02, & + 0.109103285, 3.38223181E-04, 7.78146312E-02, & + 0.114558451, 1.53707610E-02, 8.07851329E-02, & + 9.58787426E-02, 6.05889633E-02, 0.114558451, & + 4.04635631E-02, 5.83609045E-02, 5.83609045E-02] + + print *, "Test binomial_distribution_pmf_int8" + seed = 630852741 + call random_seed(seed, get) + n = 50_int8 + p = 0.6 + x1 = binom_rvs(n, p) + x2 = reshape(binom_rvs(n, p, 12_int8), [3,4]) + res(:,1) = binom_pmf(x1, n, p) + res(:, 2:5) = binom_pmf(x2, n, p) + call check(all(abs(res - reshape(ans, [3,5])) < sptol), & + msg="binomial_distribution_pmf_int8 failed", warn=warn) + end subroutine test_binom_pmf_int8 + + subroutine test_binom_pmf_int16 + integer(int16) :: x1, x2(3,4), n + integer :: i, m + integer :: seed, get + real :: p, res(3,5) + real :: ans(15) = [7.78146312E-02, 7.78146312E-02, 7.78146312E-02, & + 0.109103285, 3.38223181E-04, 7.78146312E-02, & + 0.114558451, 1.53707610E-02, 8.07851329E-02, & + 9.58787426E-02, 6.05889633E-02, 0.114558451, & + 4.04635631E-02, 5.83609045E-02, 5.83609045E-02] + + print *, "Test binomial_distribution_pmf_int16" + seed = 630852741 + call random_seed(seed, get) + n = 50_int16 + p = 0.6 + x1 = binom_rvs(n, p) + x2 = reshape(binom_rvs(n, p, 12_int16), [3,4]) + res(:,1) = binom_pmf(x1, n, p) + res(:, 2:5) = binom_pmf(x2, n, p) + call check(all(abs(res - reshape(ans, [3,5])) < sptol), & + msg="binomial_distribution_pmf_int16 failed", warn=warn) + end subroutine test_binom_pmf_int16 + + subroutine test_binom_pmf_int32 + integer(int32) :: x1, x2(3,4), n + integer :: i, m + integer :: seed, get + real :: p, res(3,5) + real :: ans(15) = [7.78146312E-02, 7.78146312E-02, 7.78146312E-02, & + 0.109103285, 3.38223181E-04, 7.78146312E-02, & + 0.114558451, 1.53707610E-02, 8.07851329E-02, & + 9.58787426E-02, 6.05889633E-02, 0.114558451, & + 4.04635631E-02, 5.83609045E-02, 5.83609045E-02] + + print *, "Test binomial_distribution_pmf_int32" + seed = 630852741 + call random_seed(seed, get) + n = 50_int32 + p = 0.6 + x1 = binom_rvs(n, p) + x2 = reshape(binom_rvs(n, p, 12_int32), [3,4]) + res(:,1) = binom_pmf(x1, n, p) + res(:, 2:5) = binom_pmf(x2, n, p) + call check(all(abs(res - reshape(ans, [3,5])) < sptol), & + msg="binomial_distribution_pmf_int32 failed", warn=warn) + end subroutine test_binom_pmf_int32 + + subroutine test_binom_pmf_int64 + integer(int64) :: x1, x2(3,4), n + integer :: i, m + integer :: seed, get + real :: p, res(3,5) + real :: ans(15) = [7.78146312E-02, 7.78146312E-02, 7.78146312E-02, & + 0.109103285, 3.38223181E-04, 7.78146312E-02, & + 0.114558451, 1.53707610E-02, 8.07851329E-02, & + 9.58787426E-02, 6.05889633E-02, 0.114558451, & + 4.04635631E-02, 5.83609045E-02, 5.83609045E-02] + + print *, "Test binomial_distribution_pmf_int64" + seed = 630852741 + call random_seed(seed, get) + n = 50_int64 + p = 0.6 + x1 = binom_rvs(n, p) + x2 = reshape(binom_rvs(n, p, 12_int64), [3,4]) + res(:,1) = binom_pmf(x1, n, p) + res(:, 2:5) = binom_pmf(x2, n, p) + call check(all(abs(res - reshape(ans, [3,5])) < sptol), & + msg="binomial_distribution_pmf_int64 failed", warn=warn) + end subroutine test_binom_pmf_int64 + + + subroutine test_binom_cdf_int8 + integer(int8) :: x1, x2(3,4), n + integer :: i, m + integer :: seed, get + real :: p, res(3,5) + real :: ans(15) = [0.978971064, 0.978971064, 0.978971064, & + 0.993534148, 0.872478724, 0.250010669, & + 5.09519465E-02, 0.125598967, 0.943473637, & + 0.872478724, 0.978971064, 0.595598698, & + 0.250010669, 0.125598967, 0.872478724] + + print *, "Test binomial_distribution_cdf_int8" + seed = 17428396 + call random_seed(seed, get) + n = 20_int8 + p = 0.4 + x1 = binom_rvs(n, p) + x2 = reshape(binom_rvs(n, p, 12_int8), [3,4]) + res(:,1) = binom_cdf(x1, n, p) + res(:, 2:5) = binom_cdf(x2, n, p) + call check(all(abs(res - reshape(ans, [3,5])) < sptol), & + msg="binomial_distribution_cdf_int8 failed", warn=warn) + end subroutine test_binom_cdf_int8 + + subroutine test_binom_cdf_int16 + integer(int16) :: x1, x2(3,4), n + integer :: i, m + integer :: seed, get + real :: p, res(3,5) + real :: ans(15) = [0.978971064, 0.978971064, 0.978971064, & + 0.993534148, 0.872478724, 0.250010669, & + 5.09519465E-02, 0.125598967, 0.943473637, & + 0.872478724, 0.978971064, 0.595598698, & + 0.250010669, 0.125598967, 0.872478724] + + print *, "Test binomial_distribution_cdf_int16" + seed = 17428396 + call random_seed(seed, get) + n = 20_int16 + p = 0.4 + x1 = binom_rvs(n, p) + x2 = reshape(binom_rvs(n, p, 12_int16), [3,4]) + res(:,1) = binom_cdf(x1, n, p) + res(:, 2:5) = binom_cdf(x2, n, p) + call check(all(abs(res - reshape(ans, [3,5])) < sptol), & + msg="binomial_distribution_cdf_int16 failed", warn=warn) + end subroutine test_binom_cdf_int16 + + subroutine test_binom_cdf_int32 + integer(int32) :: x1, x2(3,4), n + integer :: i, m + integer :: seed, get + real :: p, res(3,5) + real :: ans(15) = [0.978971064, 0.978971064, 0.978971064, & + 0.993534148, 0.872478724, 0.250010669, & + 5.09519465E-02, 0.125598967, 0.943473637, & + 0.872478724, 0.978971064, 0.595598698, & + 0.250010669, 0.125598967, 0.872478724] + + print *, "Test binomial_distribution_cdf_int32" + seed = 17428396 + call random_seed(seed, get) + n = 20_int32 + p = 0.4 + x1 = binom_rvs(n, p) + x2 = reshape(binom_rvs(n, p, 12_int32), [3,4]) + res(:,1) = binom_cdf(x1, n, p) + res(:, 2:5) = binom_cdf(x2, n, p) + call check(all(abs(res - reshape(ans, [3,5])) < sptol), & + msg="binomial_distribution_cdf_int32 failed", warn=warn) + end subroutine test_binom_cdf_int32 + + subroutine test_binom_cdf_int64 + integer(int64) :: x1, x2(3,4), n + integer :: i, m + integer :: seed, get + real :: p, res(3,5) + real :: ans(15) = [0.978971064, 0.978971064, 0.978971064, & + 0.993534148, 0.872478724, 0.250010669, & + 5.09519465E-02, 0.125598967, 0.943473637, & + 0.872478724, 0.978971064, 0.595598698, & + 0.250010669, 0.125598967, 0.872478724] + + print *, "Test binomial_distribution_cdf_int64" + seed = 17428396 + call random_seed(seed, get) + n = 20_int64 + p = 0.4 + x1 = binom_rvs(n, p) + x2 = reshape(binom_rvs(n, p, 12_int64), [3,4]) + res(:,1) = binom_cdf(x1, n, p) + res(:, 2:5) = binom_cdf(x2, n, p) + call check(all(abs(res - reshape(ans, [3,5])) < sptol), & + msg="binomial_distribution_cdf_int64 failed", warn=warn) + end subroutine test_binom_cdf_int64 + + +end program test_distribution \ No newline at end of file