From c9682d0a14cc516c6349ace6db9a1b0b3a2050b2 Mon Sep 17 00:00:00 2001 From: HugoMVale <57530119+HugoMVale@users.noreply.github.com> Date: Thu, 8 Jun 2023 08:40:01 +0200 Subject: [PATCH 1/5] replace calls to error_stop by NaN, convert procedures to pure when possible, otherwise explicit impure --- ...stdlib_stats_distribution_exponential.fypp | 65 ++++++++++--------- 1 file changed, 34 insertions(+), 31 deletions(-) diff --git a/src/stdlib_stats_distribution_exponential.fypp b/src/stdlib_stats_distribution_exponential.fypp index 01d4e8eb8..498758975 100644 --- a/src/stdlib_stats_distribution_exponential.fypp +++ b/src/stdlib_stats_distribution_exponential.fypp @@ -1,8 +1,8 @@ #:include "common.fypp" #:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES module stdlib_stats_distribution_exponential + use ieee_arithmetic, only: ieee_value, ieee_quiet_nan use stdlib_kinds, only : sp, dp, xdp, qp, int32 - use stdlib_error, only : error_stop use stdlib_random, only : dist_rand use stdlib_stats_distribution_uniform, only : uni=>rvs_uniform @@ -71,7 +71,7 @@ module stdlib_stats_distribution_exponential contains - subroutine zigset + impure subroutine zigset ! Marsaglia & Tsang generator for random normals & random exponentials. ! Translated from C by Alan Miller (amiller@bigpond.net.au) ! @@ -90,7 +90,7 @@ contains de = 7.697117470131487_dp te = de - !tables for random exponetials + ! tables for random exponentials q = ve * exp(de) ke(0) = int((de / q) * M2, kind = int32) ke(1) = 0 @@ -112,7 +112,7 @@ contains #:for k1, t1 in REAL_KINDS_TYPES - function rvs_exp_0_${t1[0]}$${k1}$( ) result(res) + impure function rvs_exp_0_${t1[0]}$${k1}$( ) result(res) ! ! Standard exponential random variate (lambda=1) ! @@ -122,8 +122,8 @@ contains if(.not. zig_exp_initialized ) call zigset iz = 0 - jz = dist_rand(1_int32) !32bit random integer - iz = iand( jz, 255 ) !random integer in [0, 255] + jz = dist_rand(1_int32) ! 32bit random integer + iz = iand( jz, 255 ) ! random integer in [0, 255] if( abs( jz ) < ke(iz) ) then res = abs(jz) * we(iz) else @@ -153,18 +153,19 @@ contains #:for k1, t1 in REAL_KINDS_TYPES - function rvs_exp_${t1[0]}$${k1}$(lambda) result(res) + impure elemental function rvs_exp_${t1[0]}$${k1}$(lambda) result(res) ! ! Exponential distributed random variate ! ${t1}$, intent(in) :: lambda ${t1}$ :: res - - if(lambda <= 0.0_${k1}$) call error_stop("Error(rvs_exp): Exponen" & - //"tial distribution lambda parameter must be greater than zero") - res = rvs_exp_0_${t1[0]}$${k1}$( ) - res = res / lambda + if (lambda <= 0._${k1}$) then + res = ieee_value(1._${k1}$, ieee_quiet_nan) + else + res = rvs_exp_0_${t1[0]}$${k1}$( ) + res = res / lambda + end if end function rvs_exp_${t1[0]}$${k1}$ #:endfor @@ -173,7 +174,7 @@ contains #:for k1, t1 in CMPLX_KINDS_TYPES - function rvs_exp_${t1[0]}$${k1}$(lambda) result(res) + impure elemental function rvs_exp_${t1[0]}$${k1}$(lambda) result(res) ${t1}$, intent(in) :: lambda ${t1}$ :: res real(${k1}$) :: tr, ti @@ -189,15 +190,17 @@ contains #:for k1, t1 in REAL_KINDS_TYPES - function rvs_exp_array_${t1[0]}$${k1}$(lambda, array_size) result(res) + impure function rvs_exp_array_${t1[0]}$${k1}$(lambda, array_size) result(res) ${t1}$, intent(in) :: lambda integer, intent(in) :: array_size ${t1}$ :: res(array_size), x, re ${t1}$, parameter :: r = 7.69711747013104972_${k1}$ integer :: jz, iz, i - if(lambda <= 0.0_${k1}$) call error_stop("Error(rvs_exp_array): Exp" & - //"onential distribution lambda parameter must be greater than zero") + if (lambda <= 0._${k1}$) then + res = ieee_value(1._${k1}$, ieee_quiet_nan) + return + end if if(.not. zig_exp_initialized) call zigset do i = 1, array_size @@ -235,7 +238,7 @@ contains #:for k1, t1 in CMPLX_KINDS_TYPES - function rvs_exp_array_${t1[0]}$${k1}$(lambda, array_size) result(res) + impure function rvs_exp_array_${t1[0]}$${k1}$(lambda, array_size) result(res) ${t1}$, intent(in) :: lambda integer, intent(in) :: array_size ${t1}$ :: res(array_size) @@ -255,18 +258,18 @@ contains #:for k1, t1 in REAL_KINDS_TYPES - impure elemental function pdf_exp_${t1[0]}$${k1}$(x, lambda) result(res) + elemental function pdf_exp_${t1[0]}$${k1}$(x, lambda) result(res) ! ! Exponential Distribution Probability Density Function ! ${t1}$, intent(in) :: x, lambda real(${k1}$) :: res - if(lambda <= 0.0_${k1}$) call error_stop("Error(pdf_exp): Expon" & - //"ential distribution lambda parameter must be greater than zero") - if(x < 0.0_${k1}$) call error_stop("Error(pdf_exp): Exponential" & - //" distribution variate x must be non-negative") - res = exp(- x * lambda) * lambda + if ((lambda <= 0._${k1}$) .or. (x < 0._${k1}$)) then + res = ieee_value(1._${k1}$, ieee_quiet_nan) + else + res = exp(- x * lambda) * lambda + end if end function pdf_exp_${t1[0]}$${k1}$ #:endfor @@ -275,7 +278,7 @@ contains #:for k1, t1 in CMPLX_KINDS_TYPES - impure elemental function pdf_exp_${t1[0]}$${k1}$(x, lambda) result(res) + elemental function pdf_exp_${t1[0]}$${k1}$(x, lambda) result(res) ${t1}$, intent(in) :: x, lambda real(${k1}$) :: res @@ -289,18 +292,18 @@ contains #:for k1, t1 in REAL_KINDS_TYPES - impure elemental function cdf_exp_${t1[0]}$${k1}$(x, lambda) result(res) + elemental function cdf_exp_${t1[0]}$${k1}$(x, lambda) result(res) ! ! Exponential Distribution Cumulative Distribution Function ! ${t1}$, intent(in) :: x, lambda real(${k1}$) :: res - if(lambda <= 0.0_${k1}$) call error_stop("Error(cdf_exp): Expon" & - //"ential distribution lambda parameter must be greater than zero") - if(x < 0.0_${k1}$) call error_stop("Error(cdf_exp): Exponential" & - //" distribution variate x must be non-negative") - res = 1.0_${k1}$ - exp(- x * lambda) + if ((lambda <= 0._${k1}$) .or. (x < 0._${k1}$)) then + res = ieee_value(1._${k1}$, ieee_quiet_nan) + else + res = 1.0_${k1}$ - exp(- x * lambda) + end if end function cdf_exp_${t1[0]}$${k1}$ #:endfor @@ -309,7 +312,7 @@ contains #:for k1, t1 in CMPLX_KINDS_TYPES - impure elemental function cdf_exp_${t1[0]}$${k1}$(x, lambda) result(res) + elemental function cdf_exp_${t1[0]}$${k1}$(x, lambda) result(res) ${t1}$, intent(in) :: x, lambda real(${k1}$) :: res From f6b2b320121c602ef77997b2983c2caab237e55b Mon Sep 17 00:00:00 2001 From: HugoMVale <57530119+HugoMVale@users.noreply.github.com> Date: Thu, 8 Jun 2023 17:18:30 +0200 Subject: [PATCH 2/5] improve readibility, fix equation rendering, update specs --- .../stdlib_stats_distribution_exponential.md | 47 ++++++++++--------- 1 file changed, 26 insertions(+), 21 deletions(-) diff --git a/doc/specs/stdlib_stats_distribution_exponential.md b/doc/specs/stdlib_stats_distribution_exponential.md index 4be87df65..f94fac5d9 100644 --- a/doc/specs/stdlib_stats_distribution_exponential.md +++ b/doc/specs/stdlib_stats_distribution_exponential.md @@ -14,16 +14,17 @@ Experimental ### Description -An exponential distribution is the distribution of time between events in a Poisson point process. The inverse scale parameter `lambda` specifies the average time between events, also called the rate of events. +An exponential distribution is the distribution of time between events in a Poisson point process. The inverse scale parameter `lambda` specifies the average time between events ($\lambda$), also called the rate of events. -Without argument the function returns a random sample from the standard exponential distribution `E(1)` with `lambda = 1`. +Without argument, the function returns a random sample from the standard exponential distribution $E(\lambda=1)$. -With a single argument, the function returns a random sample from the exponential distribution `E(lambda)`. +With a single argument, the function returns a random sample from the exponential distribution $E(\lambda=\text{lambda})$. For complex arguments, the real and imaginary parts are sampled independently of each other. -With two arguments the function returns a rank one array of exponentially distributed random variates. +With two arguments, the function returns a rank-1 array of exponentially distributed random variates. -Note: the algorithm used for generating exponetial random variates is fundamentally limited to double precision. Ref.: Marsaglia, G. & Tsang, W.W. (2000) `The ziggurat method for generating random variables', J. Statist. Software, v5(8). +@note +The algorithm used for generating exponential random variates is fundamentally limited to double precision.[^1] ### Syntax @@ -31,17 +32,19 @@ Note: the algorithm used for generating exponetial random variates is fundamenta ### Class -Function +Elemental function ### Arguments -`lambda`: optional argument has `intent(in)` and is a scalar of type `real` or `complex`. The value of `lambda` has to be non-negative. +`lambda`: optional argument has `intent(in)` and is a scalar of type `real` or `complex`. +If `lambda` is `real`, its value must be positive. If `lambda` is `complex`, both the real and imaginary components must be positive. `array_size`: optional argument has `intent(in)` and is a scalar of type `integer` with default kind. ### Return value -The result is a scalar or rank one array with a size of `array_size`, and has the same type of `lambda`. +The result is a scalar or rank-1 array with a size of `array_size`, and the same type as `lambda`. +If `lambda` is non-positive, the result is `NaN`. ### Example @@ -57,15 +60,13 @@ Experimental ### Description -The probability density function (pdf) of the single real variable exponential distribution: +The probability density function (pdf) of the single real variable exponential distribution is: -$$f(x)=\begin{cases} \lambda e^{-\lambda x} &x\geqslant 0 \\\\ 0 &x< 0\end{}$$ +$$f(x)=\begin{cases} \lambda e^{-\lambda x} &x\geqslant 0 \\\\ 0 &x< 0\end{cases}$$ -For a complex variable (x + y i) with independent real x and imaginary y parts, the joint probability density function -is the product of the corresponding marginal pdf of real and imaginary pdf (for more details, see -"Probability and Random Processes with Applications to Signal Processing and Communications", 2nd ed., Scott L. Miller and Donald Childers, 2012, p.197): +For a complex variable $z=(x + y i)$ with independent real $x$ and imaginary $y$ parts, the joint probability density function is the product of the corresponding real and imaginary marginal pdfs:[^2] -$$f(x+\mathit{i}y)=f(x)f(y)=\begin{cases} \lambda_{x} \lambda_{y} e^{-(\lambda_{x} x + \lambda_{y} y)} &x\geqslant 0, y\geqslant 0 \\\\ 0 &otherwise\end{}$$ +$$f(x+\mathit{i}y)=f(x)f(y)=\begin{cases} \lambda_{x} \lambda_{y} e^{-(\lambda_{x} x + \lambda_{y} y)} &x\geqslant 0, y\geqslant 0 \\\\ 0 &\text{otherwise}\end{cases}$$ ### Syntax @@ -80,12 +81,13 @@ Elemental function `x`: has `intent(in)` and is a scalar of type `real` or `complex`. `lambda`: has `intent(in)` and is a scalar of type `real` or `complex`. +If `lambda` is `real`, its value must be positive. If `lambda` is `complex`, both the real and imaginary components must be positive. All arguments must have the same type. ### Return value -The result is a scalar or an array, with a shape conformable to arguments, and has the same type of input arguments. +The result is a scalar or an array, with a shape conformable to the arguments, and the same type as the input arguments. If `lambda` is non-positive, the result is `NaN`. ### Example @@ -103,13 +105,11 @@ Experimental Cumulative distribution function (cdf) of the single real variable exponential distribution: -$$F(x)=\begin{cases}1 - e^{-\lambda x} &x\geqslant 0 \\\\ 0 &x< 0\end{}$$ +$$F(x)=\begin{cases}1 - e^{-\lambda x} &x\geqslant 0 \\\\ 0 &x< 0\end{cases}$$ -For a complex variable (x + y i) with independent real x and imaginary y parts, the joint cumulative distribution -function is the product of corresponding marginal cdf of real and imaginary cdf (for more details, see -"Probability and Random Processes with Applications to Signal Processing and Communications", 2nd ed., Scott L. Miller and Donald Childers, 2012, p.197): +For a complex variable $z=(x + y i)$ with independent real $x$ and imaginary $y$ parts, the joint cumulative distribution function is the product of corresponding real and imaginary marginal cdfs:[^2] -$$F(x+\mathit{i}y)=F(x)F(y)=\begin{cases} (1 - e^{-\lambda_{x} x})(1 - e^{-\lambda_{y} y}) &x\geqslant 0, \;\; y\geqslant 0 \\\\ 0 &otherwise \end{}$$ +$$F(x+\mathit{i}y)=F(x)F(y)=\begin{cases} (1 - e^{-\lambda_{x} x})(1 - e^{-\lambda_{y} y}) &x\geqslant 0, \;\; y\geqslant 0 \\\\ 0 & \text{otherwise} \end{cases}$$ ### Syntax @@ -124,15 +124,20 @@ Elemental function `x`: has `intent(in)` and is a scalar of type `real` or `complex`. `lambda`: has `intent(in)` and is a scalar of type `real` or `complex`. +If `lambda` is `real`, its value must be positive. If `lambda` is `complex`, both the real and imaginary components must be positive. All arguments must have the same type. ### Return value -The result is a scalar or an array, with a shape conformable to arguments, and has the same type of input arguments. +The result is a scalar or an array, with a shape conformable to the arguments, and the same type as the input arguments. If `lambda` is non-positive, the result is `NaN`. ### Example ```fortran {!example/stats_distribution_exponential/example_exponential_cdf.f90!} ``` + +[^1] Marsaglia, George, and Wai Wan Tsang. "The ziggurat method for generating random variables." _Journal of statistical software_ 5 (2000): 1-7. + +[^2] Miller, Scott, and Donald Childers. _Probability and random processes: With applications to signal processing and communications_. Academic Press, 2012 (p. 197). From 47dfe53262fab2d5c8885b7460c5b74310f4805b Mon Sep 17 00:00:00 2001 From: HugoMVale <57530119+HugoMVale@users.noreply.github.com> Date: Thu, 8 Jun 2023 17:56:19 +0200 Subject: [PATCH 3/5] clean up + add examples with NaN and concurrent call --- .../example_exponential_pdf.f90 | 71 ++++++++++++------- 1 file changed, 47 insertions(+), 24 deletions(-) diff --git a/example/stats_distribution_exponential/example_exponential_pdf.f90 b/example/stats_distribution_exponential/example_exponential_pdf.f90 index f8eb7a636..77c534f7f 100644 --- a/example/stats_distribution_exponential/example_exponential_pdf.f90 +++ b/example/stats_distribution_exponential/example_exponential_pdf.f90 @@ -1,39 +1,62 @@ program example_exponential_pdf use stdlib_random, only: random_seed use stdlib_stats_distribution_exponential, only: exp_pdf => pdf_exp, & - rexp => rvs_exp + rexp => rvs_exp implicit none - real :: x(2, 3, 4), a(2, 3, 4) + real, dimension(2, 3, 4) :: x, lambda + real :: xsum complex :: scale - integer :: seed_put, seed_get + integer :: seed_put, seed_get, i seed_put = 1234567 call random_seed(seed_put, seed_get) - print *, exp_pdf(1.0, 1.0) !a probability density at 1.0 in standard expon - -! 0.367879450 - - print *, exp_pdf(2.0, 2.0) !a probability density at 2.0 with lambda=2.0 - -! 3.66312787E-02 - - x = reshape(rexp(0.5, 24), [2, 3, 4]) ! standard expon random variates array - a(:, :, :) = 0.5 - print *, exp_pdf(x, a) ! a rank 3 standard expon probability density - -! 0.457115263 0.451488823 0.492391467 0.485233188 0.446215510 -! 0.401670188 0.485127628 0.316924453 0.418474048 0.483173639 -! 0.307366133 0.285812140 0.448017836 0.426440030 0.403896868 -! 0.334653258 0.410376132 0.485370994 0.333617479 0.263791025 -! 0.249779820 0.457159877 0.495636940 0.482243657 - + ! probability density at x=1.0 in standard exponential + print *, exp_pdf(1.0, 1.0) + ! 0.367879450 + + ! probability density at x=2.0 with lambda=2.0 + print *, exp_pdf(2.0, 2.0) + ! 3.66312787E-02 + + ! probability density at x=2.0 with lambda=-1.0 (out of range) + print *, exp_pdf(2.0, -1.0) + ! NaN + + ! standard exponential random variates array + x = reshape(rexp(0.5, 24), [2, 3, 4]) + + ! a rank-3 exponential probability density + lambda(:, :, :) = 0.5 + print *, exp_pdf(x, lambda) + ! 0.349295378 0.332413018 0.470253497 0.443498343 0.317152828 + ! 0.208242029 0.443112582 8.07073265E-02 0.245337561 0.436016470 + ! 7.14025944E-02 5.33841923E-02 0.322308093 0.264558554 0.212898195 + ! 0.100339092 0.226891592 0.444002301 9.91026312E-02 3.87373678E-02 + ! 3.11400592E-02 0.349431813 0.482774824 0.432669312 + + ! probability density array where lambda<=0.0 for certain elements + print *, exp_pdf([1.0, 1.0, 1.0], [1.0, 0.0, -1.0]) + ! 0.367879450 NaN NaN + + ! `pdf_exp` is pure and, thus, can be called concurrently + xsum = 0.0 + do concurrent (i=1:size(x,3)) + xsum = xsum + sum(exp_pdf(x(:,:,i), lambda(:,:,i))) + end do + print *, xsum + ! 6.45566940 + + ! complex exponential probability density function at (1.5,1.0) with real part + ! of lambda=1.0 and imaginary part of lambda=2.0 scale = (1.0, 2.) print *, exp_pdf((1.5, 1.0), scale) -! a complex expon probability density function at (1.5,1.0) with real part -!of lambda=1.0 and imaginary part of lambda=2.0 + ! 6.03947677E-02 -! 6.03947677E-02 + ! As above, but with lambda%re < 0 + scale = (-1.0, 2.) + print *, exp_pdf((1.5, 1.0), scale) + ! NaN end program example_exponential_pdf From 94ed13fe2ae1be8115101c47db5d2f19d3324a65 Mon Sep 17 00:00:00 2001 From: HugoMVale <57530119+HugoMVale@users.noreply.github.com> Date: Thu, 8 Jun 2023 18:05:36 +0200 Subject: [PATCH 4/5] pdf and cdf should be zero if x<0 --- src/stdlib_stats_distribution_exponential.fypp | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/stdlib_stats_distribution_exponential.fypp b/src/stdlib_stats_distribution_exponential.fypp index 498758975..7e3a43ae4 100644 --- a/src/stdlib_stats_distribution_exponential.fypp +++ b/src/stdlib_stats_distribution_exponential.fypp @@ -265,8 +265,10 @@ contains ${t1}$, intent(in) :: x, lambda real(${k1}$) :: res - if ((lambda <= 0._${k1}$) .or. (x < 0._${k1}$)) then + if (lambda <= 0._${k1}$) then res = ieee_value(1._${k1}$, ieee_quiet_nan) + else if (x < 0._${k1}$) then + res = 0._${k1}$ else res = exp(- x * lambda) * lambda end if @@ -299,8 +301,10 @@ contains ${t1}$, intent(in) :: x, lambda real(${k1}$) :: res - if ((lambda <= 0._${k1}$) .or. (x < 0._${k1}$)) then + if (lambda <= 0._${k1}$) then res = ieee_value(1._${k1}$, ieee_quiet_nan) + else if (x < 0._${k1}$) then + res = 0._${k1}$ else res = 1.0_${k1}$ - exp(- x * lambda) end if From afc1844607616b307b9a9a73d37971011b13e6e2 Mon Sep 17 00:00:00 2001 From: HugoMVale <57530119+HugoMVale@users.noreply.github.com> Date: Thu, 8 Jun 2023 22:28:45 +0200 Subject: [PATCH 5/5] clean up + add examples with NaN and concurrent call --- .../example_exponential_cdf.f90 | 59 +++++++++++++------ 1 file changed, 41 insertions(+), 18 deletions(-) diff --git a/example/stats_distribution_exponential/example_exponential_cdf.f90 b/example/stats_distribution_exponential/example_exponential_cdf.f90 index cc6f2c74f..9f108a686 100644 --- a/example/stats_distribution_exponential/example_exponential_cdf.f90 +++ b/example/stats_distribution_exponential/example_exponential_cdf.f90 @@ -4,38 +4,61 @@ program example_exponential_cdf rexp => rvs_exp implicit none - real :: x(2, 3, 4), a(2, 3, 4) + real, dimension(2, 3, 4) :: x, lambda + real :: xsum complex :: scale - integer :: seed_put, seed_get + integer :: seed_put, seed_get, i seed_put = 1234567 call random_seed(seed_put, seed_get) - print *, exp_cdf(1.0, 1.0) ! a standard exponential cumulative at 1.0 + ! standard exponential cumulative distribution at x=1.0 + print *, exp_cdf(1.0, 1.0) + ! 0.632120550 -! 0.632120550 + ! cumulative distribution at x=2.0 with lambda=2 + print *, exp_cdf(2.0, 2.0) + ! 0.981684387 - print *, exp_cdf(2.0, 2.0) ! a cumulative at 2.0 with lambda=2 - -! 0.981684387 + ! cumulative distribution at x=2.0 with lambda=-1.0 (out of range) + print *, exp_cdf(2.0, -1.0) + ! NaN + ! standard exponential random variates array x = reshape(rexp(0.5, 24), [2, 3, 4]) -! standard exponential random variates array - a(:, :, :) = 0.5 - print *, exp_cdf(x, a) ! a rank 3 array of standard exponential cumulative -! 8.57694745E-02 9.70223546E-02 1.52170658E-02 2.95336246E-02 -! 0.107568979 0.196659625 2.97447443E-02 0.366151094 0.163051903 -! 3.36527228E-02 0.385267735 0.428375721 0.103964329 0.147119939 -! 0.192206264 0.330693483 0.179247737 2.92580128E-02 0.332765043 -! 0.472417951 0.500440359 8.56802464E-02 8.72612000E-03 3.55126858E-02 + ! a rank-3 exponential cumulative distribution + lambda(:, :, :) = 0.5 + print *, exp_cdf(x, lambda) + ! 0.301409245 0.335173965 5.94930053E-02 0.113003314 + ! 0.365694344 0.583515942 0.113774836 0.838585377 + ! 0.509324908 0.127967060 0.857194781 0.893231630 + ! 0.355383813 0.470882893 0.574203610 0.799321830 + ! 0.546216846 0.111995399 0.801794767 0.922525287 + ! 0.937719882 0.301136374 3.44503522E-02 0.134661376 + + ! cumulative distribution array where lambda<=0.0 for certain elements + print *, exp_cdf([1.0, 1.0, 1.0], [1.0, 0.0, -1.0]) + ! 0.632120550 NaN NaN + + ! `cdf_exp` is pure and, thus, can be called concurrently + xsum = 0.0 + do concurrent (i=1:size(x,3)) + xsum = xsum + sum(exp_cdf(x(:,:,i), lambda(:,:,i))) + end do + print *, xsum + ! 11.0886612 + ! complex exponential cumulative distribution at (0.5,0.5) with real part of + ! lambda=0.5 and imaginary part of lambda=1.0 scale = (0.5, 1.0) print *, exp_cdf((0.5, 0.5), scale) -!complex exponential cumulative distribution at (0.5,0.5) with real part of -!lambda=0.5 and imaginary part of lambda=1.0 + ! 8.70351046E-02 -! 8.70351046E-02 + ! As above, but with lambda%im < 0 + scale = (1.0, -2.0) + print *, exp_cdf((1.5, 1.0), scale) + ! NaN end program example_exponential_cdf