From 43f88eb9fd7c2f143069ec3019e25e4bd9ddf023 Mon Sep 17 00:00:00 2001 From: "Sebastian G. Mutz" Date: Thu, 8 May 2025 11:51:50 +0200 Subject: [PATCH 1/7] modified parameterisation for stats exp dist. procedures --- .../example_exponential_cdf.f90 | 58 ++++--- .../example_exponential_pdf.f90 | 55 ++++--- .../example_exponential_rvs.f90 | 26 ++- ...stdlib_stats_distribution_exponential.fypp | 70 ++++---- test/stats/test_distribution_exponential.fypp | 154 +++++++++--------- 5 files changed, 193 insertions(+), 170 deletions(-) diff --git a/example/stats_distribution_exponential/example_exponential_cdf.f90 b/example/stats_distribution_exponential/example_exponential_cdf.f90 index 9f108a686..f319662c4 100644 --- a/example/stats_distribution_exponential/example_exponential_cdf.f90 +++ b/example/stats_distribution_exponential/example_exponential_cdf.f90 @@ -4,60 +4,72 @@ program example_exponential_cdf rexp => rvs_exp implicit none - real, dimension(2, 3, 4) :: x, lambda + real, dimension(2, 3, 4) :: x, loc, scale real :: xsum - complex :: scale + complex :: cloc, cscale integer :: seed_put, seed_get, i seed_put = 1234567 call random_seed(seed_put, seed_get) - ! standard exponential cumulative distribution at x=1.0 - print *, exp_cdf(1.0, 1.0) + ! standard exponential cumulative distribution at x=1.0 with loc=0.0, scale=1.0 + print *, exp_cdf(1.0, 0.0, 1.0) ! 0.632120550 - ! cumulative distribution at x=2.0 with lambda=2 - print *, exp_cdf(2.0, 2.0) + ! cumulative distribution at x=2.0 with loc=0.0 and scale=0.5 (equivalent of lambda=2) + print *, exp_cdf(2.0, 0.0, 0.5) ! 0.981684387 - ! cumulative distribution at x=2.0 with lambda=-1.0 (out of range) - print *, exp_cdf(2.0, -1.0) + ! cumulative distribution at x=2.5 with loc=0.5 and scale=0.5 (equivalent of lambda=2) + print *, exp_cdf(2.5, 0.5, 0.5) + ! 0.981684387 + + ! cumulative distribution at x=2.0 with loc=0.0 and scale=-1.0 (out of range) + print *, exp_cdf(2.0, 0.0, -1.0) ! NaN + ! cumulative distribution at x=0.5 with loc=1.0 and scale=1.0, putting x below the minimum + print *, exp_cdf(0.5, 1.0, 1.0) + ! 0.00000000 + ! standard exponential random variates array - x = reshape(rexp(0.5, 24), [2, 3, 4]) + x = reshape(rexp(0.0, 2.0, 24), [2, 3, 4]) ! a rank-3 exponential cumulative distribution - lambda(:, :, :) = 0.5 - print *, exp_cdf(x, lambda) + loc(:, :, :) = 0.0 + scale(:, :, :) = 2.0 + print *, exp_cdf(x, loc, scale) ! 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 + ! 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]) + ! cumulative distribution array where scale<=0.0 for certain elements + print *, exp_cdf([1.0, 1.0, 1.0], [0.0, 0.0, 0.0], [1.0, 0.0, -1.0]) ! 0.632120550 NaN NaN - ! `cdf_exp` is pure and, thus, can be called concurrently + ! `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))) + xsum = xsum + sum(exp_cdf(x(:,:,i), loc(:,:,i), scale(:,:,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.0, 2) with real part of + ! scale=2 and imaginary part of scale=1.0 + cloc = (0.0, 0.0) + cscale = (2, 1.0) + print *, exp_cdf((0.5, 0.5), cloc, cscale) ! 8.70351046E-02 - ! As above, but with lambda%im < 0 - scale = (1.0, -2.0) - print *, exp_cdf((1.5, 1.0), scale) + ! As above, but with scale%im < 0 + cloc = (0.0, 0.0) + cscale = (1.0, -2.0) + print *, exp_cdf((1.5, 1.0), cloc, cscale) ! NaN end program example_exponential_cdf diff --git a/example/stats_distribution_exponential/example_exponential_pdf.f90 b/example/stats_distribution_exponential/example_exponential_pdf.f90 index 77c534f7f..3f4ef6c2e 100644 --- a/example/stats_distribution_exponential/example_exponential_pdf.f90 +++ b/example/stats_distribution_exponential/example_exponential_pdf.f90 @@ -4,59 +4,66 @@ program example_exponential_pdf rexp => rvs_exp implicit none - real, dimension(2, 3, 4) :: x, lambda + real, dimension(2, 3, 4) :: x, loc, scale real :: xsum - complex :: scale + complex :: cloc, cscale integer :: seed_put, seed_get, i seed_put = 1234567 call random_seed(seed_put, seed_get) - ! probability density at x=1.0 in standard exponential - print *, exp_pdf(1.0, 1.0) + ! probability density at x=1.0 with loc=0 and in standard exponential + print *, exp_pdf(1.0, 0.0, 1.0) ! 0.367879450 - ! probability density at x=2.0 with lambda=2.0 - print *, exp_pdf(2.0, 2.0) + ! probability density at x=2.0 with loc=0.0 and scale=0.5 (lambda=2.0) + print *, exp_pdf(2.0, 0.0, 0.5) ! 3.66312787E-02 - ! probability density at x=2.0 with lambda=-1.0 (out of range) - print *, exp_pdf(2.0, -1.0) + ! probability density at x=1.5 with loc=0.5 and scale=0.5 (lambda=2.0) + print *, exp_pdf(2.5, 0.5, 0.5) + ! 3.66312787E-02 + + ! probability density at x=2.0 with loc=0.0 and scale=-1.0 (out of range) + print *, exp_pdf(2.0, 0.0, -1.0) ! NaN - ! standard exponential random variates array - x = reshape(rexp(0.5, 24), [2, 3, 4]) + ! standard exponential random variates array + x = reshape(rexp(0.0, 2.0, 24), [2, 3, 4]) ! a rank-3 exponential probability density - lambda(:, :, :) = 0.5 - print *, exp_pdf(x, lambda) + loc(:, :, :) = 0.0 + scale(:, :, :) = 2.0 + print *, exp_pdf(x, loc, scale) ! 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 + ! 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]) + ! probability density array where scale<=0.0 for certain elements (loc = 0.0) + print *, exp_pdf([1.0, 1.0, 1.0], [0.0, 0.0, 0.0], [1.0, 0.0, -1.0]) ! 0.367879450 NaN NaN - ! `pdf_exp` is pure and, thus, can be called concurrently + ! `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))) + xsum = xsum + sum(exp_pdf(x(:,:,i), loc(:,:,i), scale(:,:,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) + ! complex exponential probability density function at (1.5, 0.0, 1.0) with real part + ! of scale=1.0 and imaginary part of scale=0.5 + cloc = (0.0, 0.0) + cscale = (1.0, 0.5) + print *, exp_pdf((1.5, 1.0), cloc, cscale) ! 6.03947677E-02 - ! As above, but with lambda%re < 0 - scale = (-1.0, 2.) - print *, exp_pdf((1.5, 1.0), scale) + ! As above, but with scale%re < 0 + cloc = (0.0, 0.0) + cscale = (-1.0, 2.0) + print *, exp_pdf((1.5, 1.0), cloc, cscale) ! NaN end program example_exponential_pdf diff --git a/example/stats_distribution_exponential/example_exponential_rvs.f90 b/example/stats_distribution_exponential/example_exponential_rvs.f90 index 551d2fdce..af189f109 100644 --- a/example/stats_distribution_exponential/example_exponential_rvs.f90 +++ b/example/stats_distribution_exponential/example_exponential_rvs.f90 @@ -3,30 +3,28 @@ program example_exponential_rvs use stdlib_stats_distribution_exponential, only: rexp => rvs_exp implicit none - complex :: scale + complex :: cloc, cscale integer :: seed_put, seed_get seed_put = 1234567 call random_seed(seed_put, seed_get) print *, rexp() !single standard exponential random variate - ! 0.358690143 - print *, rexp(2.0) !exponential random variate with lambda=2.0 - -! 0.816459715 - - print *, rexp(0.3, 10) !an array of 10 variates with lambda=0.3 + print *, rexp(0.6, 0.2) !exponential random variate with loc=0.6 and scale=0.5 (lambda=2) +! 0.681645989 -! 1.84008647E-02 3.59742008E-02 0.136567295 0.262772143 3.62352766E-02 -! 0.547133625 0.213591918 4.10784185E-02 0.583882213 0.671128035 + print *, rexp(0.0, 3.0, 10) !an array of 10 variates with loc=0.0 and scale=3.0 (lambda=1/3) +! 0.184008643 0.359742016 1.36567295 2.62772131 0.362352759 +! 5.47133636 2.13591909 0.410784155 5.83882189 6.71128035 - scale = (2.0, 0.7) - print *, rexp(scale) -!single complex exponential random variate with real part of lambda=2.0; -!imagainary part of lambda=0.7 + cloc = (0.0, 0.0) + cscale = (0.5, 1.6) + print *, rexp(cloc, cscale) +!single complex exponential random variate with real part of scale=0.5 (lambda=2.0); +!imagainary part of scale=1.6 (lambda=0.625) -! (1.41435969,4.081114382E-02) +! (0.219550118,1.01847279) end program example_exponential_rvs diff --git a/src/stdlib_stats_distribution_exponential.fypp b/src/stdlib_stats_distribution_exponential.fypp index 7e3a43ae4..ffebb3299 100644 --- a/src/stdlib_stats_distribution_exponential.fypp +++ b/src/stdlib_stats_distribution_exponential.fypp @@ -114,7 +114,7 @@ contains #:for k1, t1 in REAL_KINDS_TYPES impure function rvs_exp_0_${t1[0]}$${k1}$( ) result(res) ! - ! Standard exponential random variate (lambda=1) + ! Standard exponential random variate (lambda=1; scale=1/lambda=1) ! ${t1}$ :: res, x ${t1}$, parameter :: r = 7.69711747013104972_${k1}$ @@ -153,18 +153,18 @@ contains #:for k1, t1 in REAL_KINDS_TYPES - impure elemental function rvs_exp_${t1[0]}$${k1}$(lambda) result(res) + impure elemental function rvs_exp_${t1[0]}$${k1}$(loc, scale) result(res) ! ! Exponential distributed random variate ! - ${t1}$, intent(in) :: lambda + ${t1}$, intent(in) :: loc, scale ${t1}$ :: res - if (lambda <= 0._${k1}$) then + if (scale <= 0._${k1}$) then res = ieee_value(1._${k1}$, ieee_quiet_nan) else res = rvs_exp_0_${t1[0]}$${k1}$( ) - res = res / lambda + res = res * scale + loc end if end function rvs_exp_${t1[0]}$${k1}$ @@ -174,13 +174,13 @@ contains #:for k1, t1 in CMPLX_KINDS_TYPES - impure elemental function rvs_exp_${t1[0]}$${k1}$(lambda) result(res) - ${t1}$, intent(in) :: lambda + impure elemental function rvs_exp_${t1[0]}$${k1}$(loc, scale) result(res) + ${t1}$, intent(in) :: loc, scale ${t1}$ :: res real(${k1}$) :: tr, ti - tr = rvs_exp_r${k1}$(lambda % re) - ti = rvs_exp_r${k1}$(lambda % im) + tr = rvs_exp_r${k1}$(loc%re, scale%re) + ti = rvs_exp_r${k1}$(loc%im, scale%im) res = cmplx(tr, ti, kind=${k1}$) end function rvs_exp_${t1[0]}$${k1}$ @@ -190,14 +190,14 @@ contains #:for k1, t1 in REAL_KINDS_TYPES - impure function rvs_exp_array_${t1[0]}$${k1}$(lambda, array_size) result(res) - ${t1}$, intent(in) :: lambda + impure function rvs_exp_array_${t1[0]}$${k1}$(loc, scale, array_size) result(res) + ${t1}$, intent(in) :: loc, scale integer, intent(in) :: array_size ${t1}$ :: res(array_size), x, re ${t1}$, parameter :: r = 7.69711747013104972_${k1}$ integer :: jz, iz, i - if (lambda <= 0._${k1}$) then + if (scale <= 0._${k1}$) then res = ieee_value(1._${k1}$, ieee_quiet_nan) return end if @@ -228,7 +228,7 @@ contains end if end do L1 endif - res(i) = re / lambda + res(i) = re * scale + loc end do end function rvs_exp_array_${t1[0]}$${k1}$ @@ -238,16 +238,16 @@ contains #:for k1, t1 in CMPLX_KINDS_TYPES - impure function rvs_exp_array_${t1[0]}$${k1}$(lambda, array_size) result(res) - ${t1}$, intent(in) :: lambda + impure function rvs_exp_array_${t1[0]}$${k1}$(loc, scale, array_size) result(res) + ${t1}$, intent(in) :: loc, scale integer, intent(in) :: array_size ${t1}$ :: res(array_size) integer :: i real(${k1}$) :: tr, ti do i = 1, array_size - tr = rvs_exp_r${k1}$(lambda % re) - ti = rvs_exp_r${k1}$(lambda % im) + tr = rvs_exp_r${k1}$(loc%re, scale%re) + ti = rvs_exp_r${k1}$(loc%im, scale%im) res(i) = cmplx(tr, ti, kind=${k1}$) end do end function rvs_exp_array_${t1[0]}$${k1}$ @@ -258,19 +258,19 @@ contains #:for k1, t1 in REAL_KINDS_TYPES - elemental function pdf_exp_${t1[0]}$${k1}$(x, lambda) result(res) + elemental function pdf_exp_${t1[0]}$${k1}$(x, loc, scale) result(res) ! ! Exponential Distribution Probability Density Function ! - ${t1}$, intent(in) :: x, lambda + ${t1}$, intent(in) :: x, loc, scale real(${k1}$) :: res - if (lambda <= 0._${k1}$) then + if (scale <= 0._${k1}$) then res = ieee_value(1._${k1}$, ieee_quiet_nan) - else if (x < 0._${k1}$) then + else if (x < loc) then res = 0._${k1}$ else - res = exp(- x * lambda) * lambda + res = exp(- (x - loc) / scale) / scale end if end function pdf_exp_${t1[0]}$${k1}$ @@ -280,12 +280,12 @@ contains #:for k1, t1 in CMPLX_KINDS_TYPES - elemental function pdf_exp_${t1[0]}$${k1}$(x, lambda) result(res) - ${t1}$, intent(in) :: x, lambda + elemental function pdf_exp_${t1[0]}$${k1}$(x, loc, scale) result(res) + ${t1}$, intent(in) :: x, loc, scale real(${k1}$) :: res - res = pdf_exp_r${k1}$(x % re, lambda % re) - res = res * pdf_exp_r${k1}$(x % im, lambda % im) + res = pdf_exp_r${k1}$(x%re, loc%re, scale%re) + res = res * pdf_exp_r${k1}$(x%im, loc%im, scale%im) end function pdf_exp_${t1[0]}$${k1}$ #:endfor @@ -294,19 +294,19 @@ contains #:for k1, t1 in REAL_KINDS_TYPES - elemental function cdf_exp_${t1[0]}$${k1}$(x, lambda) result(res) + elemental function cdf_exp_${t1[0]}$${k1}$(x, loc, scale) result(res) ! ! Exponential Distribution Cumulative Distribution Function ! - ${t1}$, intent(in) :: x, lambda + ${t1}$, intent(in) :: x, loc, scale real(${k1}$) :: res - if (lambda <= 0._${k1}$) then + if (scale <= 0._${k1}$) then res = ieee_value(1._${k1}$, ieee_quiet_nan) - else if (x < 0._${k1}$) then + else if (x < loc) then res = 0._${k1}$ else - res = 1.0_${k1}$ - exp(- x * lambda) + res = 1.0_sp - exp(- (x - loc) / scale) end if end function cdf_exp_${t1[0]}$${k1}$ @@ -316,12 +316,12 @@ contains #:for k1, t1 in CMPLX_KINDS_TYPES - elemental function cdf_exp_${t1[0]}$${k1}$(x, lambda) result(res) - ${t1}$, intent(in) :: x, lambda + elemental function cdf_exp_${t1[0]}$${k1}$(x, loc, scale) result(res) + ${t1}$, intent(in) :: x, loc, scale real(${k1}$) :: res - res = cdf_exp_r${k1}$(x % re, lambda % re) - res = res * cdf_exp_r${k1}$(x % im, lambda % im) + res = cdf_exp_r${k1}$(x%re, loc%re, scale%re) + res = res * cdf_exp_r${k1}$(x%im, loc%im, scale%im) end function cdf_exp_${t1[0]}$${k1}$ #:endfor diff --git a/test/stats/test_distribution_exponential.fypp b/test/stats/test_distribution_exponential.fypp index 01d3f5bb5..5df3b29fb 100644 --- a/test/stats/test_distribution_exponential.fypp +++ b/test/stats/test_distribution_exponential.fypp @@ -53,7 +53,7 @@ contains print *, "Test exponential random generator with chi-squared" freq = 0 do i = 1, num - j = 1000 * (1 - exp(- expon_rvs(1.0))) + j = 1000 * (1 - exp(- expon_rvs(0.0, 1.0))) freq(j) = freq(j) + 1 end do chisq = 0.0_dp @@ -72,46 +72,46 @@ contains #:for k1, t1 in RC_KINDS_TYPES subroutine test_expon_rvs_${t1[0]}$${k1}$ - ${t1}$ :: res(10), scale + ${t1}$ :: res(10), loc, scale integer, parameter :: k = 5 integer :: i integer :: seed, get #:if t1[0] == "r" #! for real type ${t1}$, parameter :: ans(10) = & - [0.609680481289574416337018192280083895_${k1}$, & - 0.137541023585612635452927558314210417_${k1}$, & - 0.134921508232253721063879462841820585_${k1}$, & - 1.33766060689229752493171569464417802_${k1}$, & - 0.111148487576340881943792737729381770_${k1}$, & - 0.533951653963536868966836361020492979_${k1}$, & - 1.96897428558727671799033487332053483_${k1}$, & - 0.371111977992924465160247867364281152_${k1}$, & - 0.811918715695663687862785688290993341_${k1}$, & - 0.404637854946697868759504975362991277_${k1}$] + [ 1.37178108290154243675829093263018876_${k1}$, & + 0.309467303067628429769087006206973456_${k1}$, & + 0.303573393522570872393728791394096334_${k1}$, & + 3.00973636550766943109636031294940040_${k1}$, & + 0.250084097046766984373533659891108982_${k1}$, & + 1.20139122141795795517538181229610927_${k1}$, & + 4.43019214257137261547825346497120336_${k1}$, & + 0.835001950484080046610557701569632627_${k1}$, & + 1.82681711031524329769126779865473509_${k1}$, & + 0.910435173630070204708886194566730410_${k1}$] #:else #! for complex type ${t1}$, parameter :: ans(10) = & - [(1.30645817419194517786503898345732266_${k1}$, & - 0.158701181060322271676454874977935106_${k1}$), & - (0.289117517640543687994027420375329869_${k1}$, & - 1.54345454641418945184428733997405138_${k1}$), & - (0.238175330520730461308127295134389521_${k1}$, & - 0.616098062265619464192503493485184250_${k1}$), & - (4.21923061197273582426500329997257485_${k1}$, & - 0.428206128453374382877209077728016710_${k1}$), & - (1.73982581934785075970596933205212874_${k1}$, & - 0.466889832630805233184044202341912994_${k1}$), & - (2.22649889847873832288931745486999202_${k1}$, & - 0.879109337848515628785697537331053851_${k1}$), & - (8.76802198822945553859296653951917464_${k1}$, & - 0.200128045239398311139211728004738688_${k1}$), & - (0.694821947760945587572020290930855262_${k1}$, & - 0.101964167346166995492113143812345625_${k1}$), & - (0.141476585024528208770330398432893829_${k1}$, & - 3.989655879458742013468417133900891716E-0002_${k1}$), & - (2.10676792861163792685325850990401309_${k1}$, & - 0.249356813451327473065187125310051027_${k1}$)] + [(0.640164505354053137153869101894088070_${k1}$, & + 0.268204995991944639133208738712710328_${k1}$), & + (0.141667583643866407117073435983911608_${k1}$, & + 2.60843818343998017361684560455614716_${k1}$), & + (0.116705911955157926040982374615850854_${k1}$, & + 1.04120572522889689448533090398996145_${k1}$), & + (2.06742299986664055388985161698656149_${k1}$, & + 0.723668357086202707062483341360348315_${k1}$), & + (0.852514651480446872255924972705542983_${k1}$, & + 0.789043817146060844081034701957833041_${k1}$), & + (1.09098446025458177821576555288629603_${k1}$, & + 1.48569478096399141264782883808948111_${k1}$), & + (4.29633077423243321391055360436439499_${k1}$, & + 0.338216396454583145825267820328008412_${k1}$), & + (0.340462754402863337910289942556119029_${k1}$, & + 0.172319442815022222381671213042864120_${k1}$), & + (6.932352666201882229746189523211795805E-0002_${k1}$, & + 6.742518436285274002761624956292507704E-0002_${k1}$), & + (1.03231628501970258415809666985296648_${k1}$, & + 0.421413014732743429480166241773986277_${k1}$)] #:endif print *, "Test exponential_distribution_rvs_${t1[0]}$${k1}$" @@ -120,16 +120,18 @@ contains #:if t1[0] == "r" #! for real type + loc = 0.0_${k1}$ scale = 1.5_${k1}$ #:else #! for complex type + loc = (0.0_${k1}$, 0.0_${k1}$) scale = (0.7_${k1}$, 1.3_${k1}$) #:endif do i = 1, k - res(i) = expon_rvs(scale) ! 1 dummy + res(i) = expon_rvs(loc, scale) ! 1 dummy end do - res(6:10) = expon_rvs(scale, k) ! 2 dummies + res(6:10) = expon_rvs(loc, scale, k) ! 2 dummies call check(all(abs(res - ans) < ${k1}$tol), & msg="exponential_distribution_rvs_${t1[0]}$${k1}$ failed", warn=warn) end subroutine test_expon_rvs_${t1[0]}$${k1}$ @@ -142,45 +144,45 @@ contains #:for k1, t1 in RC_KINDS_TYPES subroutine test_expon_pdf_${t1[0]}$${k1}$ - ${t1}$ :: x1, x2(3,4), scale + ${t1}$ :: x1, x2(3,4), loc, scale integer :: seed, get real(${k1}$) :: res(3,5) #:if t1[0] == "r" #! for real type real(${k1}$), parameter :: ans(15) = & - [0.362692289054629718342313806171796533_${k1}$, & - 0.362692289054629718342313806171796533_${k1}$, & - 0.362692289054629718342313806171796533_${k1}$, & - 1.44877092399186122284289290051705535_${k1}$, & - 1.08871761038277651996081144393335589_${k1}$, & - 0.203258408490339213767867275283195750_${k1}$, & - 0.730004225568590859263284124264208147_${k1}$, & - 0.237394827760488451509080387833146683_${k1}$, & - 0.301732182586179598102005265289645959_${k1}$, & - 1.35079274124711914255014934401469271_${k1}$, & - 0.416578245043239337295928202660090263_${k1}$, & - 1.44039177901335374382803898226703593_${k1}$, & - 0.196044829271295768265275728683411055_${k1}$, & - 0.271373826917613661285112379170965958_${k1}$, & - 1.00108987409617105109732206933052664_${k1}$] + [0.161196572913168763707695024965242901_${k1}$, & + 0.161196572913168763707695024965242901_${k1}$, & + 0.161196572913168763707695024965242901_${k1}$, & + 0.643898188440827210152396844674246781_${k1}$, & + 0.483874493503456231093693975081491447_${k1}$, & + 9.03370704401507616746076779036425611E-0002_${k1}$, & + 0.324446322474929270783681833006314743_${k1}$, & + 0.105508812337994867337369061259176307_${k1}$, & + 0.134103192260524265823113451239842656_${k1}$, & + 0.600352329443164063355621930673196730_${k1}$, & + 0.185145886685884149909301423404484569_${k1}$, & + 0.640174124005934997256906214340904849_${k1}$, & + 8.71310352316870081179003238592937984E-0002_${k1}$, & + 0.120610589741161627237827724075984870_${k1}$, & + 0.444928832931631578265476475258011855_${k1}$] #:else #! for complex type real(${k1}$), parameter :: ans(15) = & - [0.112097715784191810518066563334849515_${k1}$, & - 0.112097715784191810518066563334849515_${k1}$, & - 0.112097715784191810518066563334849515_${k1}$, & - 4.72087485401191174735651518020251204E-0002_${k1}$, & - 3.69705018439006691768174449531170720E-0002_${k1}$, & - 8.69498969681198520061798177185735738E-0002_${k1}$, & - 0.128007654288233028296342302153338001_${k1}$, & - 1.13496395875758374774198906169957218E-0002_${k1}$, & - 0.294260498264128747413785056084385424_${k1}$, & - 4.66169813179250908948018478030960097E-0002_${k1}$, & - 2.84438693906889813143446828488861951E-0002_${k1}$, & - 0.161859307815385236742977105439660254_${k1}$, & - 4.22904796362406579112752522035325397E-0002_${k1}$, & - 0.176117981883470250164040199296778089_${k1}$, & - 0.107352342201327219885025541854724060_${k1}$] + [0.486535224757776955373552792251950892_${k1}$, & + 0.486535224757776955373552792251950892_${k1}$, & + 0.486535224757776955373552792251950892_${k1}$, & + 0.204899082205378114034570971362956231_${k1}$, & + 0.160462247586374432191047938164570635_${k1}$, & + 0.377386705590797968776822125514642235_${k1}$, & + 0.555588777292678074202874575318307322_${k1}$, & + 4.92605884877423501624127196918217120E-0002_${k1}$, & + 1.27717230149361435509455319481070053_${k1}$, & + 0.202330648081272095897577464423159785_${k1}$, & + 0.123454294230420925843509908198290816_${k1}$, & + 0.702514356837609534474727020137414207_${k1}$, & + 0.183552428976738966628798837688943286_${k1}$, & + 0.764400963035895183003646698336710458_${k1}$, & + 0.465938985248816058528756692077795393_${k1}$] #:endif print *, "Test exponential_distribution_pdf_${t1[0]}$${k1}$" @@ -188,16 +190,18 @@ contains call random_seed(seed, get) #:if t1[0] == "r" #! for real type + loc = 0.0_${k1}$ scale = 1.5_${k1}$ #:else #! for complex type + loc = (0.0_${k1}$, 0.0_${k1}$) scale = (0.3_${k1}$, 1.6_${k1}$) #:endif - x1 = expon_rvs(scale) - x2 = reshape(expon_rvs(scale, 12), [3,4]) - res(:,1) = expon_pdf(x1, scale) - res(:, 2:5) = expon_pdf(x2, scale) + x1 = expon_rvs(loc, scale) + x2 = reshape(expon_rvs(loc, scale, 12), [3,4]) + res(:,1) = expon_pdf(x1, loc, scale) + res(:, 2:5) = expon_pdf(x2, loc, scale) call check(all(abs(res - reshape(ans, [3,5])) < ${k1}$tol), & msg="exponential_distribution_pdf_${t1[0]}$${k1}$ failed", warn=warn) end subroutine test_expon_pdf_${t1[0]}$${k1}$ @@ -210,7 +214,7 @@ contains #:for k1, t1 in RC_KINDS_TYPES subroutine test_expon_cdf_${t1[0]}$${k1}$ - ${t1}$ :: x1, x2(3,4), scale + ${t1}$ :: x1, x2(3,4), loc, scale integer :: seed, get real(${k1}$) :: res(3,5) #:if t1[0] == "r" @@ -256,15 +260,17 @@ contains #:if t1[0] == "r" #! for real type + loc = 0.0_${k1}$ scale = 2.0_${k1}$ #:else + loc = (0.0_${k1}$, 0.0_${k1}$) scale = (1.3_${k1}$, 2.1_${k1}$) #:endif - x1 = expon_rvs(scale) - x2 = reshape(expon_rvs(scale, 12), [3,4]) - res(:,1) = expon_cdf(x1, scale) - res(:, 2:5) = expon_cdf(x2, scale) + x1 = expon_rvs(loc, scale) + x2 = reshape(expon_rvs(loc, scale, 12), [3,4]) + res(:,1) = expon_cdf(x1, loc, scale) + res(:, 2:5) = expon_cdf(x2, loc, scale) call check(all(abs(res - reshape(ans,[3,5])) < ${k1}$tol), & msg="exponential_distribution_cdf_${t1[0]}$${k1}$ failed", warn=warn) end subroutine test_expon_cdf_${t1[0]}$${k1}$ From e7984cb93a97dbebce57f73d622afc3215e3aa4f Mon Sep 17 00:00:00 2001 From: Sebastian Date: Tue, 13 May 2025 19:12:16 +0100 Subject: [PATCH 2/7] reintroduced old interface through wrappers and updated specs --- .../stdlib_stats_distribution_exponential.md | 60 ++++- .../example_exponential_cdf.f90 | 8 + .../example_exponential_pdf.f90 | 10 +- .../example_exponential_rvs.f90 | 32 ++- ...stdlib_stats_distribution_exponential.fypp | 171 +++++++++++++- test/stats/test_distribution_exponential.fypp | 222 +++++++++++------- 6 files changed, 399 insertions(+), 104 deletions(-) diff --git a/doc/specs/stdlib_stats_distribution_exponential.md b/doc/specs/stdlib_stats_distribution_exponential.md index b11ca6a6c..27f5919a4 100644 --- a/doc/specs/stdlib_stats_distribution_exponential.md +++ b/doc/specs/stdlib_stats_distribution_exponential.md @@ -15,14 +15,19 @@ 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 (\(\lambda\)), also called the rate of events. +The inverse `scale` parameter `lambda` specifies the average time between events (\(\lambda\)), also called the rate of events. The location `loc` specifies the value by which the distribution is shifted. -Without argument, the function returns a random sample from the standard exponential distribution \(E(\lambda=1)\). +Without argument, the function returns a random sample from the unshifted standard exponential distribution \(E(\lambda=1)\) or \(E(loc=0, scale=1)\). -With a single argument, the function returns a random sample from the exponential distribution \(E(\lambda=\text{lambda})\). +With a single argument of type `real`, 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-1 array of exponentially distributed random variates. +With one argument of type `real` and one argument of type `integer`, the function returns a rank-1 array of exponentially distributed random variates for (E(\lambda=\text{lambda})\). + +With two arguments of type `real`, the function returns a random sample from the exponential distribution \(E(loc=loc, scale=scale)\). +For complex arguments, the real and imaginary parts are sampled independently of each other. + +With two arguments of type `real` and one argument of type `integer`, the function returns a rank-1 array of exponentially distributed random variates for \(E(loc=loc, scale=scale)\). @note The algorithm used for generating exponential random variates is fundamentally limited to double precision.[^1] @@ -31,6 +36,10 @@ The algorithm used for generating exponential random variates is fundamentally l `result = ` [[stdlib_stats_distribution_exponential(module):rvs_exp(interface)]] `([lambda] [[, array_size]])` +or + +`result = ` [[stdlib_stats_distribution_exponential(module):rvs_exp(interface)]] `([loc, scale] [[, array_size]])` + ### Class Elemental function @@ -40,13 +49,21 @@ Elemental function `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. +`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`. +If `scale` is `real`, its value must be positive. If `scale` 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-1 array with a size of `array_size`, and the same type as `lambda`. +If `lambda` is passed, 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`. +If `loc` and `scale` are passed, the result is a scalar or rank-1 array with a size of `array_size`, and the same type as `scale`. +If `scale` is non-positive, the result is `NaN`. + ### Example ```fortran @@ -69,10 +86,16 @@ For a complex variable \(z=(x + y i)\) with independent real \(x\) and imaginary $$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}$$ +Instead of of the inverse scale parameter `lambda`, it is possible to pass `loc` and `scale`, where \(scale = \frac{1}{\lambda}\) and `loc` specifies the value by which the distribution is shifted. + ### Syntax `result = ` [[stdlib_stats_distribution_exponential(module):pdf_exp(interface)]] `(x, lambda)` +or + +`result = ` [[stdlib_stats_distribution_exponential(module):pdf_exp(interface)]] `(x, loc, scale)` + ### Class Elemental function @@ -84,11 +107,20 @@ Elemental function `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. +`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`. +If `scale` is `real`, its value must be positive. If `scale` 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 the arguments, and the same type as the input arguments. If `lambda` is non-positive, the result is `NaN`. +If `lambda` is passed, 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`. + + +If `loc` and `scale` are passed, the result is a scalar or an array, with a shape conformable to the arguments, and the same type as the input arguments. If `scale` is non-positive, the result is `NaN`. + ### Example @@ -112,10 +144,16 @@ For a complex variable \(z=(x + y i)\) with independent real \(x\) and imaginar $$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}$$ +Instead of of the inverse scale parameter `lambda`, it is possible to pass `loc` and `scale`, where \(scale = \frac{1}{\lambda}\) and `loc` specifies the value by which the distribution is shifted. + ### Syntax `result = ` [[stdlib_stats_distribution_exponential(module):cdf_exp(interface)]] `(x, lambda)` +or + +`result = ` [[stdlib_stats_distribution_exponential(module):cdf_exp(interface)]] `(x, loc, scale)` + ### Class Elemental function @@ -127,11 +165,19 @@ Elemental function `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. +`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`. +If `scale` is `real`, its value must be positive. If `scale` 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 the arguments, and the same type as the input arguments. If `lambda` is non-positive, the result is `NaN`. +If `lamba` is passed, 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`. + + +If `loc` and `scale` are passed, the result is a scalar or an array, with a shape conformable to the arguments, and the same type as the input arguments. If `scale` is non-positive, the result is `NaN`. ### Example diff --git a/example/stats_distribution_exponential/example_exponential_cdf.f90 b/example/stats_distribution_exponential/example_exponential_cdf.f90 index f319662c4..e99edaed9 100644 --- a/example/stats_distribution_exponential/example_exponential_cdf.f90 +++ b/example/stats_distribution_exponential/example_exponential_cdf.f90 @@ -12,10 +12,18 @@ program example_exponential_cdf seed_put = 1234567 call random_seed(seed_put, seed_get) + ! standard exponential cumulative distribution at x=1.0 with lambda=1.0 + print *, exp_cdf(1.0, 1.0) + ! 0.632120550 + ! standard exponential cumulative distribution at x=1.0 with loc=0.0, scale=1.0 print *, exp_cdf(1.0, 0.0, 1.0) ! 0.632120550 + ! cumulative distribution at x=2.0 with lambda=2 + print *, exp_cdf(2.0, 2.0) + ! 0.981684387 + ! cumulative distribution at x=2.0 with loc=0.0 and scale=0.5 (equivalent of lambda=2) print *, exp_cdf(2.0, 0.0, 0.5) ! 0.981684387 diff --git a/example/stats_distribution_exponential/example_exponential_pdf.f90 b/example/stats_distribution_exponential/example_exponential_pdf.f90 index 3f4ef6c2e..c8f9aa379 100644 --- a/example/stats_distribution_exponential/example_exponential_pdf.f90 +++ b/example/stats_distribution_exponential/example_exponential_pdf.f90 @@ -12,10 +12,18 @@ program example_exponential_pdf seed_put = 1234567 call random_seed(seed_put, seed_get) - ! probability density at x=1.0 with loc=0 and in standard exponential + ! probability density at x=1.0 with lambda=1.0 + print *, exp_pdf(1.0, 1.0) + ! 0.367879450 + + ! probability density at x=1.0 with loc=0 and scale=1.0 print *, exp_pdf(1.0, 0.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 loc=0.0 and scale=0.5 (lambda=2.0) print *, exp_pdf(2.0, 0.0, 0.5) ! 3.66312787E-02 diff --git a/example/stats_distribution_exponential/example_exponential_rvs.f90 b/example/stats_distribution_exponential/example_exponential_rvs.f90 index af189f109..65fa84a65 100644 --- a/example/stats_distribution_exponential/example_exponential_rvs.f90 +++ b/example/stats_distribution_exponential/example_exponential_rvs.f90 @@ -9,22 +9,32 @@ program example_exponential_rvs seed_put = 1234567 call random_seed(seed_put, seed_get) - print *, rexp() !single standard exponential random variate -! 0.358690143 + ! single standard exponential random variate + print *, rexp() + ! 0.358690143 - print *, rexp(0.6, 0.2) !exponential random variate with loc=0.6 and scale=0.5 (lambda=2) -! 0.681645989 + ! exponential random variate with lambda=2 + print *, rexp(2.0) + ! 0.204114929 - print *, rexp(0.0, 3.0, 10) !an array of 10 variates with loc=0.0 and scale=3.0 (lambda=1/3) -! 0.184008643 0.359742016 1.36567295 2.62772131 0.362352759 -! 5.47133636 2.13591909 0.410784155 5.83882189 6.71128035 + ! exponential random variate with loc=0 and scale=0.5 (lambda=2) + print *, rexp(0.0, 0.5) + ! 0.122672431 + ! exponential random variate with loc=0.6 and scale=0.2 (lambda=5) + print *, rexp(0.6, 0.2) + ! 0.681645989 + + ! an array of 10 variates with loc=0.0 and scale=3.0 (lambda=1/3) + print *, rexp(0.0, 3.0, 10) + ! 1.36567295 2.62772131 0.362352759 5.47133636 2.13591909 + ! 0.410784155 5.83882189 6.71128035 1.31730068 1.90963650 + + ! single complex exponential random variate with real part of scale=0.5 (lambda=2.0); + ! imagainary part of scale=1.6 (lambda=0.625) cloc = (0.0, 0.0) cscale = (0.5, 1.6) print *, rexp(cloc, cscale) -!single complex exponential random variate with real part of scale=0.5 (lambda=2.0); -!imagainary part of scale=1.6 (lambda=0.625) - -! (0.219550118,1.01847279) + ! (0.426896989,2.56968451) end program example_exponential_rvs diff --git a/src/stdlib_stats_distribution_exponential.fypp b/src/stdlib_stats_distribution_exponential.fypp index ffebb3299..cf4b43ab3 100644 --- a/src/stdlib_stats_distribution_exponential.fypp +++ b/src/stdlib_stats_distribution_exponential.fypp @@ -9,7 +9,7 @@ module stdlib_stats_distribution_exponential implicit none private - integer :: ke(0:255) + integer :: ke(0:255) real(dp) :: we(0:255), fe(0:255) logical :: zig_exp_initialized = .false. @@ -26,14 +26,26 @@ module stdlib_stats_distribution_exponential !! ([Specification](../page/specs/stdlib_stats_distribution_exponential.html# !! rvs_exp-exponential-distribution-random-variates)) !! - module procedure rvs_exp_0_rsp !0 dummy variable + module procedure rvs_exp_0_rsp !0 dummy variable + + ! new interfaces using loc and scale + + #:for k1, t1 in RC_KINDS_TYPES + module procedure rvs_exp_${t1[0]}$${k1}$ !1 dummy variable + #:endfor + + #:for k1, t1 in RC_KINDS_TYPES + module procedure rvs_exp_array_${t1[0]}$${k1}$ !2 dummy variables + #:endfor + + ! original interfaces using lambda #:for k1, t1 in RC_KINDS_TYPES - module procedure rvs_exp_${t1[0]}$${k1}$ !1 dummy variable + module procedure rvs_exp_lambda_${t1[0]}$${k1}$ !1 dummy variable #:endfor #:for k1, t1 in RC_KINDS_TYPES - module procedure rvs_exp_array_${t1[0]}$${k1}$ !2 dummy variables + module procedure rvs_exp_array_lambda_${t1[0]}$${k1}$ !2 dummy variables #:endfor end interface rvs_exp @@ -46,9 +58,17 @@ module stdlib_stats_distribution_exponential !! ([Specification](../page/specs/stdlib_stats_distribution_exponential.html# !! pdf_exp-exponential-distribution-probability-density-function)) !! + ! new interfaces using loc and scale + #:for k1, t1 in RC_KINDS_TYPES module procedure pdf_exp_${t1[0]}$${k1}$ #:endfor + + ! original interfaces using lambda + + #:for k1, t1 in RC_KINDS_TYPES + module procedure pdf_exp_lambda_${t1[0]}$${k1}$ + #:endfor end interface pdf_exp @@ -60,9 +80,17 @@ module stdlib_stats_distribution_exponential !! ([Specification](../page/specs/stdlib_stats_distribution_exponential.html# !! cdf_exp-exponential-distribution-cumulative-distribution-function)) !! + ! new interfaces using loc and scale + #:for k1, t1 in RC_KINDS_TYPES module procedure cdf_exp_${t1[0]}$${k1}$ #:endfor + + ! original interfaces using lambda + + #:for k1, t1 in RC_KINDS_TYPES + module procedure cdf_exp_lambda_${t1[0]}$${k1}$ + #:endfor end interface cdf_exp @@ -306,7 +334,7 @@ contains else if (x < loc) then res = 0._${k1}$ else - res = 1.0_sp - exp(- (x - loc) / scale) + res = 1.0_${k1}$ - exp(- (x - loc) / scale) end if end function cdf_exp_${t1[0]}$${k1}$ @@ -326,4 +354,137 @@ contains #:endfor + + + + ! + ! below: wrapper functions for interfaces using lambda (for backward compatibility) + ! + + + + + #:for k1, t1 in REAL_KINDS_TYPES + impure elemental function rvs_exp_lambda_${t1[0]}$${k1}$(lambda) result(res) + ! + ! Exponential distributed random variate + ! + ${t1}$, intent(in) :: lambda + ${t1}$ :: res + + res = rvs_exp_${t1[0]}$${k1}$(0._${k1}$, 1.0_${k1}$/lambda) + end function rvs_exp_lambda_${t1[0]}$${k1}$ + + #:endfor + + + + + #:for k1, t1 in CMPLX_KINDS_TYPES + impure elemental function rvs_exp_lambda_${t1[0]}$${k1}$(lambda) result(res) + ${t1}$, intent(in) :: lambda + ${t1}$ :: res + real(${k1}$) :: tr, ti + + tr = rvs_exp_r${k1}$(0._${k1}$, 1.0_${k1}$/lambda%re) + ti = rvs_exp_r${k1}$(0._${k1}$, 1.0_${k1}$/lambda%im) + res = cmplx(tr, ti, kind=${k1}$) + end function rvs_exp_lambda_${t1[0]}$${k1}$ + + #:endfor + + + + + #:for k1, t1 in REAL_KINDS_TYPES + impure function rvs_exp_array_lambda_${t1[0]}$${k1}$(lambda, array_size) result(res) + ${t1}$, intent(in) :: lambda + integer, intent(in) :: array_size + ${t1}$ :: res(array_size) + + res = rvs_exp_array_${t1[0]}$${k1}$(0._${k1}$, 1.0_${k1}$/lambda, array_size) + end function rvs_exp_array_lambda_${t1[0]}$${k1}$ + + #:endfor + + + + + #:for k1, t1 in CMPLX_KINDS_TYPES + impure function rvs_exp_array_lambda_${t1[0]}$${k1}$(lambda, array_size) result(res) + ${t1}$, intent(in) :: lambda + integer, intent(in) :: array_size + ${t1}$ :: res(array_size) + integer :: i + real(${k1}$) :: tr, ti + + do i = 1, array_size + tr = rvs_exp_r${k1}$(0._${k1}$, 1.0_${k1}$/lambda%re) + ti = rvs_exp_r${k1}$(0._${k1}$, 1.0_${k1}$/lambda%im) + res(i) = cmplx(tr, ti, kind=${k1}$) + end do + end function rvs_exp_array_lambda_${t1[0]}$${k1}$ + + #:endfor + + + + + #:for k1, t1 in REAL_KINDS_TYPES + elemental function pdf_exp_lambda_${t1[0]}$${k1}$(x, lambda) result(res) + ! + ! Exponential Distribution Probability Density Function + ! + ${t1}$, intent(in) :: x, lambda + real(${k1}$) :: res + + res = pdf_exp_${t1[0]}$${k1}$(x, 0._${k1}$, 1.0_${k1}$/lambda) + end function pdf_exp_lambda_${t1[0]}$${k1}$ + + #:endfor + + + + + #:for k1, t1 in CMPLX_KINDS_TYPES + elemental function pdf_exp_lambda_${t1[0]}$${k1}$(x, lambda) result(res) + ${t1}$, intent(in) :: x, lambda + real(${k1}$) :: res + + res = pdf_exp_r${k1}$(x%re, 0._${k1}$, 1.0_${k1}$/lambda%re) + res = res * pdf_exp_r${k1}$(x%im, 0._${k1}$, 1.0_${k1}$/lambda%im) + end function pdf_exp_lambda_${t1[0]}$${k1}$ + + #:endfor + + + + + #:for k1, t1 in REAL_KINDS_TYPES + elemental function cdf_exp_lambda_${t1[0]}$${k1}$(x, lambda) result(res) + ! + ! Exponential Distribution Cumulative Distribution Function + ! + ${t1}$, intent(in) :: x, lambda + real(${k1}$) :: res + + res = cdf_exp_${t1[0]}$${k1}$(x, 0._${k1}$, 1.0_${k1}$/lambda) + end function cdf_exp_lambda_${t1[0]}$${k1}$ + + #:endfor + + + + + #:for k1, t1 in CMPLX_KINDS_TYPES + elemental function cdf_exp_lambda_${t1[0]}$${k1}$(x, lambda) result(res) + ${t1}$, intent(in) :: x, lambda + real(${k1}$) :: res + + res = cdf_exp_r${k1}$(x%re, 0._${k1}$, 1.0_${k1}$/lambda%re) + res = res * cdf_exp_r${k1}$(x%im, 0._${k1}$, 1.0_${k1}$/lambda%im) + end function cdf_exp_lambda_${t1[0]}$${k1}$ + + #:endfor + end module stdlib_stats_distribution_exponential diff --git a/test/stats/test_distribution_exponential.fypp b/test/stats/test_distribution_exponential.fypp index 5df3b29fb..832df2c20 100644 --- a/test/stats/test_distribution_exponential.fypp +++ b/test/stats/test_distribution_exponential.fypp @@ -51,6 +51,25 @@ contains print *, "" print *, "Test exponential random generator with chi-squared" + + ! using interface for lambda + freq = 0 + do i = 1, num + j = 1000 * (1 - exp(- expon_rvs(1.0))) + freq(j) = freq(j) + 1 + end do + chisq = 0.0_dp + expct = num / array_size + do i = 0, array_size - 1 + chisq = chisq + (freq(i) - expct) ** 2 / expct + end do + write(*,*) "The critical values for chi-squared with 1000 d. of f. is" & + //" 1143.92" + write(*,*) "Chi-squared for exponential random generator is : ", chisq + call check((chisq < 1143.9), & + msg="exponential randomness failed chi-squared test", warn=warn) + + ! using interface for loc and scale freq = 0 do i = 1, num j = 1000 * (1 - exp(- expon_rvs(0.0, 1.0))) @@ -66,72 +85,86 @@ contains write(*,*) "Chi-squared for exponential random generator is : ", chisq call check((chisq < 1143.9), & msg="exponential randomness failed chi-squared test", warn=warn) + end subroutine test_exponential_random_generator #:for k1, t1 in RC_KINDS_TYPES subroutine test_expon_rvs_${t1[0]}$${k1}$ - ${t1}$ :: res(10), loc, scale + ${t1}$ :: res(10), lambda, loc, scale integer, parameter :: k = 5 integer :: i integer :: seed, get #:if t1[0] == "r" #! for real type ${t1}$, parameter :: ans(10) = & - [ 1.37178108290154243675829093263018876_${k1}$, & - 0.309467303067628429769087006206973456_${k1}$, & - 0.303573393522570872393728791394096334_${k1}$, & - 3.00973636550766943109636031294940040_${k1}$, & - 0.250084097046766984373533659891108982_${k1}$, & - 1.20139122141795795517538181229610927_${k1}$, & - 4.43019214257137261547825346497120336_${k1}$, & - 0.835001950484080046610557701569632627_${k1}$, & - 1.82681711031524329769126779865473509_${k1}$, & - 0.910435173630070204708886194566730410_${k1}$] + [0.609680481289574416337018192280083895_${k1}$, & + 0.137541023585612635452927558314210417_${k1}$, & + 0.134921508232253721063879462841820585_${k1}$, & + 1.33766060689229752493171569464417802_${k1}$, & + 0.111148487576340881943792737729381770_${k1}$, & + 0.533951653963536868966836361020492979_${k1}$, & + 1.96897428558727671799033487332053483_${k1}$, & + 0.371111977992924465160247867364281152_${k1}$, & + 0.811918715695663687862785688290993341_${k1}$, & + 0.404637854946697868759504975362991277_${k1}$] #:else #! for complex type ${t1}$, parameter :: ans(10) = & - [(0.640164505354053137153869101894088070_${k1}$, & - 0.268204995991944639133208738712710328_${k1}$), & - (0.141667583643866407117073435983911608_${k1}$, & - 2.60843818343998017361684560455614716_${k1}$), & - (0.116705911955157926040982374615850854_${k1}$, & - 1.04120572522889689448533090398996145_${k1}$), & - (2.06742299986664055388985161698656149_${k1}$, & - 0.723668357086202707062483341360348315_${k1}$), & - (0.852514651480446872255924972705542983_${k1}$, & - 0.789043817146060844081034701957833041_${k1}$), & - (1.09098446025458177821576555288629603_${k1}$, & - 1.48569478096399141264782883808948111_${k1}$), & - (4.29633077423243321391055360436439499_${k1}$, & - 0.338216396454583145825267820328008412_${k1}$), & - (0.340462754402863337910289942556119029_${k1}$, & - 0.172319442815022222381671213042864120_${k1}$), & - (6.932352666201882229746189523211795805E-0002_${k1}$, & - 6.742518436285274002761624956292507704E-0002_${k1}$), & - (1.03231628501970258415809666985296648_${k1}$, & - 0.421413014732743429480166241773986277_${k1}$)] + [(1.30645817419194517786503898345732266_${k1}$, & + 0.158701181060322271676454874977935106_${k1}$), & + (0.289117517640543687994027420375329869_${k1}$, & + 1.54345454641418945184428733997405138_${k1}$), & + (0.238175330520730461308127295134389521_${k1}$, & + 0.616098062265619464192503493485184250_${k1}$), & + (4.21923061197273582426500329997257485_${k1}$, & + 0.428206128453374382877209077728016710_${k1}$), & + (1.73982581934785075970596933205212874_${k1}$, & + 0.466889832630805233184044202341912994_${k1}$), & + (2.22649889847873832288931745486999202_${k1}$, & + 0.879109337848515628785697537331053851_${k1}$), & + (8.76802198822945553859296653951917464_${k1}$, & + 0.200128045239398311139211728004738688_${k1}$), & + (0.694821947760945587572020290930855262_${k1}$, & + 0.101964167346166995492113143812345625_${k1}$), & + (0.141476585024528208770330398432893829_${k1}$, & + 3.989655879458742013468417133900891716E-0002_${k1}$), & + (2.10676792861163792685325850990401309_${k1}$, & + 0.249356813451327473065187125310051027_${k1}$)] #:endif print *, "Test exponential_distribution_rvs_${t1[0]}$${k1}$" seed = 593742186 - call random_seed(seed, get) + ! set args #:if t1[0] == "r" #! for real type - loc = 0.0_${k1}$ - scale = 1.5_${k1}$ + lambda = 1.5_${k1}$ + loc = 0._${k1}$ + scale = 1.0_${k1}$/lambda #:else #! for complex type - loc = (0.0_${k1}$, 0.0_${k1}$) - scale = (0.7_${k1}$, 1.3_${k1}$) + lambda = (0.7_${k1}$, 1.3_${k1}$) + loc = (0._${k1}$, 0._${k1}$) + scale = cmplx(1.0_${k1}$/lambda%re, 1.0_${k1}$/lambda%im, kind=${k1}$) #:endif + ! tests using interface for lambda + call random_seed(seed, get) + do i = 1, k + res(i) = expon_rvs(lambda) ! 1 dummy + end do + res(6:10) = expon_rvs(lambda, k) ! 2 dummies + call check(all(abs(res - ans) < ${k1}$tol), & + msg="exponential_distribution_rvs_${t1[0]}$${k1}$ failed", warn=warn) + + ! tests using interface for loc and scale + call random_seed(seed, get) do i = 1, k - res(i) = expon_rvs(loc, scale) ! 1 dummy + res(i) = expon_rvs(loc, scale) end do - res(6:10) = expon_rvs(loc, scale, k) ! 2 dummies + res(6:10) = expon_rvs(loc, scale, k) call check(all(abs(res - ans) < ${k1}$tol), & msg="exponential_distribution_rvs_${t1[0]}$${k1}$ failed", warn=warn) end subroutine test_expon_rvs_${t1[0]}$${k1}$ @@ -144,66 +177,81 @@ contains #:for k1, t1 in RC_KINDS_TYPES subroutine test_expon_pdf_${t1[0]}$${k1}$ - ${t1}$ :: x1, x2(3,4), loc, scale + ${t1}$ :: x1, x2(3,4), lambda, loc, scale integer :: seed, get real(${k1}$) :: res(3,5) #:if t1[0] == "r" #! for real type real(${k1}$), parameter :: ans(15) = & - [0.161196572913168763707695024965242901_${k1}$, & - 0.161196572913168763707695024965242901_${k1}$, & - 0.161196572913168763707695024965242901_${k1}$, & - 0.643898188440827210152396844674246781_${k1}$, & - 0.483874493503456231093693975081491447_${k1}$, & - 9.03370704401507616746076779036425611E-0002_${k1}$, & - 0.324446322474929270783681833006314743_${k1}$, & - 0.105508812337994867337369061259176307_${k1}$, & - 0.134103192260524265823113451239842656_${k1}$, & - 0.600352329443164063355621930673196730_${k1}$, & - 0.185145886685884149909301423404484569_${k1}$, & - 0.640174124005934997256906214340904849_${k1}$, & - 8.71310352316870081179003238592937984E-0002_${k1}$, & - 0.120610589741161627237827724075984870_${k1}$, & - 0.444928832931631578265476475258011855_${k1}$] + [0.362692289054629718342313806171796533_${k1}$, & + 0.362692289054629718342313806171796533_${k1}$, & + 0.362692289054629718342313806171796533_${k1}$, & + 1.44877092399186122284289290051705535_${k1}$, & + 1.08871761038277651996081144393335589_${k1}$, & + 0.203258408490339213767867275283195750_${k1}$, & + 0.730004225568590859263284124264208147_${k1}$, & + 0.237394827760488451509080387833146683_${k1}$, & + 0.301732182586179598102005265289645959_${k1}$, & + 1.35079274124711914255014934401469271_${k1}$, & + 0.416578245043239337295928202660090263_${k1}$, & + 1.44039177901335374382803898226703593_${k1}$, & + 0.196044829271295768265275728683411055_${k1}$, & + 0.271373826917613661285112379170965958_${k1}$, & + 1.00108987409617105109732206933052664_${k1}$] #:else #! for complex type real(${k1}$), parameter :: ans(15) = & - [0.486535224757776955373552792251950892_${k1}$, & - 0.486535224757776955373552792251950892_${k1}$, & - 0.486535224757776955373552792251950892_${k1}$, & - 0.204899082205378114034570971362956231_${k1}$, & - 0.160462247586374432191047938164570635_${k1}$, & - 0.377386705590797968776822125514642235_${k1}$, & - 0.555588777292678074202874575318307322_${k1}$, & - 4.92605884877423501624127196918217120E-0002_${k1}$, & - 1.27717230149361435509455319481070053_${k1}$, & - 0.202330648081272095897577464423159785_${k1}$, & - 0.123454294230420925843509908198290816_${k1}$, & - 0.702514356837609534474727020137414207_${k1}$, & - 0.183552428976738966628798837688943286_${k1}$, & - 0.764400963035895183003646698336710458_${k1}$, & - 0.465938985248816058528756692077795393_${k1}$] + [0.112097715784191810518066563334849515_${k1}$, & + 0.112097715784191810518066563334849515_${k1}$, & + 0.112097715784191810518066563334849515_${k1}$, & + 4.72087485401191174735651518020251204E-0002_${k1}$, & + 3.69705018439006691768174449531170720E-0002_${k1}$, & + 8.69498969681198520061798177185735738E-0002_${k1}$, & + 0.128007654288233028296342302153338001_${k1}$, & + 1.13496395875758374774198906169957218E-0002_${k1}$, & + 0.294260498264128747413785056084385424_${k1}$, & + 4.66169813179250908948018478030960097E-0002_${k1}$, & + 2.84438693906889813143446828488861951E-0002_${k1}$, & + 0.161859307815385236742977105439660254_${k1}$, & + 4.22904796362406579112752522035325397E-0002_${k1}$, & + 0.176117981883470250164040199296778089_${k1}$, & + 0.107352342201327219885025541854724060_${k1}$] #:endif print *, "Test exponential_distribution_pdf_${t1[0]}$${k1}$" seed = 123987654 - call random_seed(seed, get) + + ! set args #:if t1[0] == "r" #! for real type - loc = 0.0_${k1}$ - scale = 1.5_${k1}$ + lambda = 1.5_${k1}$ + loc = 0._${k1}$ + scale = 1.0_${k1}$/lambda #:else #! for complex type - loc = (0.0_${k1}$, 0.0_${k1}$) - scale = (0.3_${k1}$, 1.6_${k1}$) + lambda = (0.3_${k1}$, 1.6_${k1}$) + loc = (0._${k1}$, 0._${k1}$) + scale = cmplx(1.0_${k1}$/lambda%re, 1.0_${k1}$/lambda%im, kind=${k1}$) #:endif + ! tests using interface for lambda + call random_seed(seed, get) + x1 = expon_rvs(lambda) + x2 = reshape(expon_rvs(lambda, 12), [3,4]) + res(:,1) = expon_pdf(x1, lambda) + res(:, 2:5) = expon_pdf(x2, lambda) + call check(all(abs(res - reshape(ans, [3,5])) < ${k1}$tol), & + msg="exponential_distribution_pdf_${t1[0]}$${k1}$ failed", warn=warn) + + ! tests using interface for loc and scale + call random_seed(seed, get) x1 = expon_rvs(loc, scale) x2 = reshape(expon_rvs(loc, scale, 12), [3,4]) res(:,1) = expon_pdf(x1, loc, scale) res(:, 2:5) = expon_pdf(x2, loc, scale) call check(all(abs(res - reshape(ans, [3,5])) < ${k1}$tol), & msg="exponential_distribution_pdf_${t1[0]}$${k1}$ failed", warn=warn) + end subroutine test_expon_pdf_${t1[0]}$${k1}$ #:endfor @@ -214,7 +262,7 @@ contains #:for k1, t1 in RC_KINDS_TYPES subroutine test_expon_cdf_${t1[0]}$${k1}$ - ${t1}$ :: x1, x2(3,4), loc, scale + ${t1}$ :: x1, x2(3,4), lambda, loc, scale integer :: seed, get real(${k1}$) :: res(3,5) #:if t1[0] == "r" @@ -256,23 +304,37 @@ contains print *, "Test exponential_distribution_cdf_${t1[0]}$${k1}$" seed = 621957438 - call random_seed(seed, get) + ! set args #:if t1[0] == "r" #! for real type - loc = 0.0_${k1}$ - scale = 2.0_${k1}$ + lambda = 2.0_${k1}$ + loc = 0._${k1}$ + scale = 1.0_${k1}$/lambda #:else - loc = (0.0_${k1}$, 0.0_${k1}$) - scale = (1.3_${k1}$, 2.1_${k1}$) + lambda = (1.3_${k1}$, 2.1_${k1}$) + loc = (0._${k1}$, 0._${k1}$) + scale = cmplx(1.0_${k1}$/lambda%re, 1.0_${k1}$/lambda%im, kind=${k1}$) #:endif + ! tests using interface for lambda + call random_seed(seed, get) + x1 = expon_rvs(lambda) + x2 = reshape(expon_rvs(lambda, 12), [3,4]) + res(:,1) = expon_cdf(x1, lambda) + res(:, 2:5) = expon_cdf(x2, lambda) + call check(all(abs(res - reshape(ans,[3,5])) < ${k1}$tol), & + msg="exponential_distribution_cdf_${t1[0]}$${k1}$ failed", warn=warn) + + ! tests using interface for loc and scale + call random_seed(seed, get) x1 = expon_rvs(loc, scale) x2 = reshape(expon_rvs(loc, scale, 12), [3,4]) res(:,1) = expon_cdf(x1, loc, scale) res(:, 2:5) = expon_cdf(x2, loc, scale) call check(all(abs(res - reshape(ans,[3,5])) < ${k1}$tol), & msg="exponential_distribution_cdf_${t1[0]}$${k1}$ failed", warn=warn) + end subroutine test_expon_cdf_${t1[0]}$${k1}$ #:endfor From ba239c6a671f9f6da5a48f57119b587155c3add9 Mon Sep 17 00:00:00 2001 From: Sebastian Date: Tue, 13 May 2025 20:40:03 +0100 Subject: [PATCH 3/7] Trigger mergeability check From 9fe477b17530ba3e3c0f7d66a791fd886e840dff Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 16 May 2025 18:27:57 +0200 Subject: [PATCH 4/7] Update doc/specs/stdlib_stats_distribution_exponential.md --- doc/specs/stdlib_stats_distribution_exponential.md | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/doc/specs/stdlib_stats_distribution_exponential.md b/doc/specs/stdlib_stats_distribution_exponential.md index 27f5919a4..38177fde6 100644 --- a/doc/specs/stdlib_stats_distribution_exponential.md +++ b/doc/specs/stdlib_stats_distribution_exponential.md @@ -116,10 +116,7 @@ All arguments must have the same type. ### Return value -If `lambda` is passed, 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`. - - -If `loc` and `scale` are passed, 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`. +The result is a scalar or an array, with a shape conformable to the arguments, and the same type as the input arguments. If non-positive `lambda` or `scale`, the result is `NaN`. ### Example From 3c5114bbc07bd18d8c613f361a84bc5489b9a0b4 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 16 May 2025 18:28:10 +0200 Subject: [PATCH 5/7] Update doc/specs/stdlib_stats_distribution_exponential.md --- doc/specs/stdlib_stats_distribution_exponential.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_stats_distribution_exponential.md b/doc/specs/stdlib_stats_distribution_exponential.md index 38177fde6..3390ed7b9 100644 --- a/doc/specs/stdlib_stats_distribution_exponential.md +++ b/doc/specs/stdlib_stats_distribution_exponential.md @@ -141,7 +141,7 @@ For a complex variable \(z=(x + y i)\) with independent real \(x\) and imaginar $$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}$$ -Instead of of the inverse scale parameter `lambda`, it is possible to pass `loc` and `scale`, where \(scale = \frac{1}{\lambda}\) and `loc` specifies the value by which the distribution is shifted. +Alternative to the inverse scale parameter `lambda`, it is possible to pass `loc` and `scale`, where \(scale = \frac{1}{\lambda}\) and `loc` specifies the value by which the distribution is shifted. ### Syntax From 3ae8b58b89d2dad816d2e70bcc6e7445137201cd Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 16 May 2025 18:28:40 +0200 Subject: [PATCH 6/7] Update doc/specs/stdlib_stats_distribution_exponential.md --- doc/specs/stdlib_stats_distribution_exponential.md | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/doc/specs/stdlib_stats_distribution_exponential.md b/doc/specs/stdlib_stats_distribution_exponential.md index 3390ed7b9..f5f2e0268 100644 --- a/doc/specs/stdlib_stats_distribution_exponential.md +++ b/doc/specs/stdlib_stats_distribution_exponential.md @@ -171,10 +171,7 @@ All arguments must have the same type. ### Return value -If `lamba` is passed, 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`. - - -If `loc` and `scale` are passed, 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`. +The result is a scalar or an array, with a shape conformable to the arguments, and the same type as the input arguments. With non-positive `lambda` or `scale`, the result is `NaN`. ### Example From 9b99b53cf6747547329056f4555a3480f5d37390 Mon Sep 17 00:00:00 2001 From: Sebastian Date: Sat, 17 May 2025 21:02:48 +0100 Subject: [PATCH 7/7] new API prioritised --- doc/specs/stdlib_stats_distribution_exponential.md | 12 ++++++------ .../example_exponential_cdf.f90 | 8 ++++---- .../example_exponential_pdf.f90 | 8 ++++---- .../example_exponential_rvs.f90 | 10 +++++----- 4 files changed, 19 insertions(+), 19 deletions(-) diff --git a/doc/specs/stdlib_stats_distribution_exponential.md b/doc/specs/stdlib_stats_distribution_exponential.md index f5f2e0268..3c6577697 100644 --- a/doc/specs/stdlib_stats_distribution_exponential.md +++ b/doc/specs/stdlib_stats_distribution_exponential.md @@ -34,11 +34,11 @@ The algorithm used for generating exponential random variates is fundamentally l ### Syntax -`result = ` [[stdlib_stats_distribution_exponential(module):rvs_exp(interface)]] `([lambda] [[, array_size]])` +`result = ` [[stdlib_stats_distribution_exponential(module):rvs_exp(interface)]] `([loc, scale] [[, array_size]])` or -`result = ` [[stdlib_stats_distribution_exponential(module):rvs_exp(interface)]] `([loc, scale] [[, array_size]])` +`result = ` [[stdlib_stats_distribution_exponential(module):rvs_exp(interface)]] `([lambda] [[, array_size]])` ### Class @@ -90,11 +90,11 @@ Instead of of the inverse scale parameter `lambda`, it is possible to pass `loc` ### Syntax -`result = ` [[stdlib_stats_distribution_exponential(module):pdf_exp(interface)]] `(x, lambda)` +`result = ` [[stdlib_stats_distribution_exponential(module):pdf_exp(interface)]] `(x, loc, scale)` or -`result = ` [[stdlib_stats_distribution_exponential(module):pdf_exp(interface)]] `(x, loc, scale)` +`result = ` [[stdlib_stats_distribution_exponential(module):pdf_exp(interface)]] `(x, lambda)` ### Class @@ -145,11 +145,11 @@ Alternative to the inverse scale parameter `lambda`, it is possible to pass `loc ### Syntax -`result = ` [[stdlib_stats_distribution_exponential(module):cdf_exp(interface)]] `(x, lambda)` +`result = ` [[stdlib_stats_distribution_exponential(module):cdf_exp(interface)]] `(x, loc, scale)` or -`result = ` [[stdlib_stats_distribution_exponential(module):cdf_exp(interface)]] `(x, loc, scale)` +`result = ` [[stdlib_stats_distribution_exponential(module):cdf_exp(interface)]] `(x, lambda)` ### Class diff --git a/example/stats_distribution_exponential/example_exponential_cdf.f90 b/example/stats_distribution_exponential/example_exponential_cdf.f90 index e99edaed9..5df68d88d 100644 --- a/example/stats_distribution_exponential/example_exponential_cdf.f90 +++ b/example/stats_distribution_exponential/example_exponential_cdf.f90 @@ -12,14 +12,14 @@ program example_exponential_cdf seed_put = 1234567 call random_seed(seed_put, seed_get) - ! standard exponential cumulative distribution at x=1.0 with lambda=1.0 - print *, exp_cdf(1.0, 1.0) - ! 0.632120550 - ! standard exponential cumulative distribution at x=1.0 with loc=0.0, scale=1.0 print *, exp_cdf(1.0, 0.0, 1.0) ! 0.632120550 + ! standard exponential cumulative distribution at x=1.0 with lambda=1.0 + print *, exp_cdf(1.0, 1.0) + ! 0.632120550 + ! cumulative distribution at x=2.0 with lambda=2 print *, exp_cdf(2.0, 2.0) ! 0.981684387 diff --git a/example/stats_distribution_exponential/example_exponential_pdf.f90 b/example/stats_distribution_exponential/example_exponential_pdf.f90 index c8f9aa379..29760e4e7 100644 --- a/example/stats_distribution_exponential/example_exponential_pdf.f90 +++ b/example/stats_distribution_exponential/example_exponential_pdf.f90 @@ -12,14 +12,14 @@ program example_exponential_pdf seed_put = 1234567 call random_seed(seed_put, seed_get) - ! probability density at x=1.0 with lambda=1.0 - print *, exp_pdf(1.0, 1.0) - ! 0.367879450 - ! probability density at x=1.0 with loc=0 and scale=1.0 print *, exp_pdf(1.0, 0.0, 1.0) ! 0.367879450 + ! probability density at x=1.0 with lambda=1.0 + 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 diff --git a/example/stats_distribution_exponential/example_exponential_rvs.f90 b/example/stats_distribution_exponential/example_exponential_rvs.f90 index 65fa84a65..2f6fb04dd 100644 --- a/example/stats_distribution_exponential/example_exponential_rvs.f90 +++ b/example/stats_distribution_exponential/example_exponential_rvs.f90 @@ -9,18 +9,18 @@ program example_exponential_rvs seed_put = 1234567 call random_seed(seed_put, seed_get) - ! single standard exponential random variate + ! single standard exponential random variate print *, rexp() ! 0.358690143 - ! exponential random variate with lambda=2 - print *, rexp(2.0) - ! 0.204114929 - ! exponential random variate with loc=0 and scale=0.5 (lambda=2) print *, rexp(0.0, 0.5) ! 0.122672431 + ! exponential random variate with lambda=2 + print *, rexp(2.0) + ! 0.204114929 + ! exponential random variate with loc=0.6 and scale=0.2 (lambda=5) print *, rexp(0.6, 0.2) ! 0.681645989