rvs_normal - normal distribution random variatesExperimental
A normal continuous random variate distribution, also known as Gaussian, or Gauss or Laplace-Gauss distribution.
The location loc specifies the mean or expectation (). The scale specifies the standard deviation ().
Without argument, the function returns a standard normal distributed random variate .
With two arguments, the function returns a normal distributed random variate . For complex arguments, the real and imaginary parts are independent of each other.
With three arguments, the function returns a rank-1 array of normal distributed random variates.
With one or two arguments where the first is array_size, the function returns a rank-1 array of standard normal distributed random variates . The mold argument determines the output type and kind; it is optional only for real(dp) (and defaults to real(dp) when omitted), but required for all other types.
Note
The algorithm used for generating exponential random variates is fundamentally limited to double precision.1
result = rvs_normal ([loc, scale] [[, array_size]])
result = rvs_normal (array_size [, mold])
Elemental function (passing both loc and scale).
loc: optional argument has intent(in) and is a scalar of type real or complex.
scale: optional argument has intent(in) and is a positive scalar of type real or complex.
array_size: optional argument has intent(in) and is a scalar of type integer. When used with loc and scale, specifies the size of the output array. When used alone or with mold, must be provided as the first argument.
mold: optional argument (only for real(dp); required for other types) has intent(in) and is a scalar of type real or complex. Used only to determine the type and kind of the output; its value is not referenced. When omitted (only allowed for real(dp)), defaults to real(dp). When provided, generates standard normal variates of the specified type and kind.
loc and scale arguments must be of the same type.
The result is a scalar or rank-1 array, with a size of array_size, and the same type as scale and loc (or same type and kind as mold when using the array_size [, mold] form; defaults to real(dp) when mold is omitted). If scale is non-positive, the result is NaN.
program example_normal_rvs
use stdlib_random, only: random_seed
use stdlib_stats_distribution_normal, only: norm => rvs_normal
implicit none
real :: a(2, 3, 4), b(2, 3, 4)
complex :: loc, scale
integer :: seed_put, seed_get
seed_put = 1234567
call random_seed(seed_put, seed_get)
print *, norm() !single standard normal random variate
! 0.563655198
print *, norm(1.0, 2.0)
!normal random variate mu=1.0, sigma=2.0
! -0.633261681
print *, norm(0.0, 1.0, 10) !an array of 10 standard norml random variates
! -3.38123664E-02 -0.190365672 0.220678389 -0.424612164 -0.249541596
! 0.865260184 1.11086845 -0.328349441 1.10873628 1.27049923
a(:, :, :) = 1.0
b(:, :, :) = 1.0
print *, norm(a, b) ! a rank 3 random variates array
!0.152776539 -7.51764774E-02 1.47208166 0.180561781 1.32407105
! 1.20383692 0.123445868 -0.455737948 -0.469808221 1.60750175
! 1.05748117 0.720934749 0.407810807 1.48165631 2.31749439
! 0.414566994 3.06084275 1.86505437 1.36338580 7.26878643E-02
! 0.178585172 1.39557445 0.828021586 0.872084975
loc = (-1.0, 2.0)
scale = (2.0, 1.0)
print *, norm(loc, scale)
!single complex normal random variate with real part of mu=-1, sigma=2;
!imagainary part of mu=2.0 and sigma=1.0
! (1.22566295,2.12518454)
end program example_normal_rvs
pdf_normal - normal distribution probability density functionExperimental
The probability density function (pdf) of the single real variable normal distribution:
For a complex varible with independent real and imaginary parts, the joint probability density function is the product of the the corresponding real and imaginary marginal pdfs:2
result = pdf_normal (x, loc, scale)
Elemental function
x: has intent(in) and is a scalar of type real or complex.
loc: has intent(in) and is a scalar of type real or complex.
scale: has intent(in) and is a positive scalar of type real or complex.
All three arguments must have the same type.
The result is a scalar or an array, with a shape conformable to the arguments, and the same type as the input arguments. If scale is non-positive, the result is NaN.
program example_normal_pdf
use stdlib_random, only: random_seed
use stdlib_stats_distribution_normal, only: norm_pdf => pdf_normal, &
norm => rvs_normal
implicit none
real, dimension(3, 4, 5) :: x, mu, sigma
real :: xsum
complex :: loc, scale
integer :: seed_put, seed_get, i
seed_put = 1234567
call random_seed(seed_put, seed_get)
! probability density at x=1.0 in standard normal
print *, norm_pdf(1.0, 0., 1.)
! 0.241970733
! probability density at x=2.0 with mu=-1.0 and sigma=2.0
print *, norm_pdf(2.0, -1.0, 2.0)
! 6.47588000E-02
! probability density at x=1.0 with mu=1.0 and sigma=-1.0 (out of range)
print *, norm_pdf(1.0, 1.0, -1.0)
! NaN
! standard normal random variates array
x = reshape(norm(0.0, 1.0, 60), [3, 4, 5])
! standard normal probability density array
mu(:, :, :) = 0.0
sigma(:, :, :) = 1.0
print *, norm_pdf(x, mu, sigma)
! 0.340346158 0.285823315 0.398714304 0.391778737 0.389345556
! 0.364551932 0.386712372 0.274370432 0.215250477 0.378006011
! 0.215760440 0.177990928 0.278640658 0.223813817 0.356875211
! 0.285167664 0.378533930 0.390739858 0.271684974 0.138273031
! 0.135456234 0.331718773 0.398283750 0.383706540
! probability density array where sigma<=0.0 for certain elements
print *, norm_pdf([1.0, 1.0, 1.0], [1.0, 1.0, 1.0], [1.0, 0.0, -1.0])
! 0.398942292 NaN NaN
! `pdf_normal` is pure and, thus, can be called concurrently
xsum = 0.0
do concurrent (i=1:size(x,3))
xsum = xsum + sum(norm_pdf(x(:,:,i), mu(:,:,i), sigma(:,:,i)))
end do
print *, xsum
! 18.0754433
! complex normal probability density function at x=(1.5,1.0) with real part
! of mu=1.0, sigma=1.0 and imaginary part of mu=-0.5, sigma=2.0
loc = (1.0, -0.5)
scale = (1.0, 2.)
print *, norm_pdf((1.5, 1.0), loc, scale)
! 5.30100204E-02
end program example_normal_pdf
cdf_normal - normal distribution cumulative distribution functionExperimental
Cumulative distribution function of the single real variable normal distribution:
For the complex variable with independent real and imaginary parts, the joint cumulative distribution function is the product of the corresponding real and imaginary marginal cdfs:2
result = cdf_normal (x, loc, scale)
Elemental function
x: has intent(in) and is a scalar of type real or complex.
loc: has intent(in) and is a scalar of type real or complex.
scale: has intent(in) and is a positive scalar of type real or complex.
All three arguments must have the same type.
The result is a scalar or an array, with a shape conformable to the arguments, and the same type as the input arguments. If scale is non-positive, the result is NaN.
program example_normal_cdf
use stdlib_random, only: random_seed
use stdlib_stats_distribution_normal, only: norm_cdf => cdf_normal, &
norm => rvs_normal
implicit none
real, dimension(2, 3, 4) :: x, mu, sigma
real :: xsum
complex :: loc, scale
integer :: seed_put, seed_get, i
seed_put = 1234567
call random_seed(seed_put, seed_get)
! standard normal cumulative probability at x=0.0
print *, norm_cdf(0.0, 0.0, 1.0)
! 0.500000000
! cumulative probability at x=2.0 with mu=-1.0 sigma=2.0
print *, norm_cdf(2.0, -1.0, 2.0)
! 0.933192849
! cumulative probability at x=1.0 with mu=1.0 and sigma=-1.0 (out of range)
print *, norm_cdf(1.0, 1.0, -1.0)
! NaN
! standard normal random variates array
x = reshape(norm(0.0, 1.0, 24), [2, 3, 4])
! standard normal cumulative array
mu(:, :, :) = 0.0
sigma(:, :, :) = 1.0
print *, norm_cdf(x, mu, sigma)
! 0.713505626 0.207069695 0.486513376 0.424511284 0.587328553
! 0.335559726 0.401470929 0.806552052 0.866687536 0.371323735
! 0.866228044 0.898046613 0.198435277 0.141147852 0.681565762
! 0.206268221 0.627057910 0.580759525 0.190364420 7.27325380E-02
! 7.08068311E-02 0.728241026 0.522919059 0.390097380
! cumulative probability array where sigma<=0.0 for certain elements
print *, norm_cdf([1.0, 1.0, 1.0], [1.0, 1.0, 1.0], [1.0, 0.0, -1.0])
! 0.500000000 NaN NaN
! `cdf_normal` is pure and, thus, can be called concurrently
xsum = 0.0
do concurrent (i=1:size(x,3))
xsum = xsum + sum(norm_cdf(x(:,:,i), mu(:,:,i), sigma(:,:,i)))
end do
print *, xsum
! 11.3751936
! complex normal cumulative distribution at x=(0.5,-0.5) with real part of
! mu=1.0, sigma=0.5 and imaginary part of mu=0.0, sigma=1.0
loc = (1.0, 0.0)
scale = (0.5, 1.0)
print *, norm_cdf((0.5, -0.5), loc, scale)
! 4.89511043E-02
end program example_normal_cdf
Marsaglia, George, and Wai Wan Tsang. "The ziggurat method for generating random variables." Journal of statistical software 5 (2000): 1-7. ↩
Miller, Scott, and Donald Childers. Probability and random processes: With applications to signal processing and communications. Academic Press, 2012 (p. 197). ↩↩