Skip to content

Commit 3d0243e

Browse files
perazzperazz
authoredMay 22, 2025··
[stats] exponential distribution procedures: use loc and scale (#991)
2 parents 09bccec + 9b99b53 commit 3d0243e

File tree

6 files changed

+442
-130
lines changed

6 files changed

+442
-130
lines changed
 

‎doc/specs/stdlib_stats_distribution_exponential.md

Lines changed: 47 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -15,20 +15,29 @@ Experimental
1515
### Description
1616

1717
An exponential distribution is the distribution of time between events in a Poisson point process.
18-
The inverse scale parameter `lambda` specifies the average time between events (\(\lambda\)), also called the rate of events.
18+
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.
1919

20-
Without argument, the function returns a random sample from the standard exponential distribution \(E(\lambda=1)\).
20+
Without argument, the function returns a random sample from the unshifted standard exponential distribution \(E(\lambda=1)\) or \(E(loc=0, scale=1)\).
2121

22-
With a single argument, the function returns a random sample from the exponential distribution \(E(\lambda=\text{lambda})\).
22+
With a single argument of type `real`, the function returns a random sample from the exponential distribution \(E(\lambda=\text{lambda})\).
2323
For complex arguments, the real and imaginary parts are sampled independently of each other.
2424

25-
With two arguments, the function returns a rank-1 array of exponentially distributed random variates.
25+
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})\).
26+
27+
With two arguments of type `real`, the function returns a random sample from the exponential distribution \(E(loc=loc, scale=scale)\).
28+
For complex arguments, the real and imaginary parts are sampled independently of each other.
29+
30+
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)\).
2631

2732
@note
2833
The algorithm used for generating exponential random variates is fundamentally limited to double precision.[^1]
2934

3035
### Syntax
3136

37+
`result = ` [[stdlib_stats_distribution_exponential(module):rvs_exp(interface)]] `([loc, scale] [[, array_size]])`
38+
39+
or
40+
3241
`result = ` [[stdlib_stats_distribution_exponential(module):rvs_exp(interface)]] `([lambda] [[, array_size]])`
3342

3443
### Class
@@ -40,13 +49,21 @@ Elemental function
4049
`lambda`: optional argument has `intent(in)` and is a scalar of type `real` or `complex`.
4150
If `lambda` is `real`, its value must be positive. If `lambda` is `complex`, both the real and imaginary components must be positive.
4251

52+
`loc`: optional argument has `intent(in)` and is a scalar of type `real` or `complex`.
53+
54+
`scale`: optional argument has `intent(in)` and is a positive scalar of type `real` or `complex`.
55+
If `scale` is `real`, its value must be positive. If `scale` is `complex`, both the real and imaginary components must be positive.
56+
4357
`array_size`: optional argument has `intent(in)` and is a scalar of type `integer` with default kind.
4458

4559
### Return value
4660

47-
The result is a scalar or rank-1 array with a size of `array_size`, and the same type as `lambda`.
61+
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`.
4862
If `lambda` is non-positive, the result is `NaN`.
4963

64+
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`.
65+
If `scale` is non-positive, the result is `NaN`.
66+
5067
### Example
5168

5269
```fortran
@@ -69,8 +86,14 @@ For a complex variable \(z=(x + y i)\) with independent real \(x\) and imaginary
6986

7087
$$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}$$
7188

89+
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.
90+
7291
### Syntax
7392

93+
`result = ` [[stdlib_stats_distribution_exponential(module):pdf_exp(interface)]] `(x, loc, scale)`
94+
95+
or
96+
7497
`result = ` [[stdlib_stats_distribution_exponential(module):pdf_exp(interface)]] `(x, lambda)`
7598

7699
### Class
@@ -84,11 +107,17 @@ Elemental function
84107
`lambda`: has `intent(in)` and is a scalar of type `real` or `complex`.
85108
If `lambda` is `real`, its value must be positive. If `lambda` is `complex`, both the real and imaginary components must be positive.
86109

110+
`loc`: has `intent(in)` and is a scalar of type `real` or `complex`.
111+
112+
`scale`: has `intent(in)` and is a positive scalar of type `real` or `complex`.
113+
If `scale` is `real`, its value must be positive. If `scale` is `complex`, both the real and imaginary components must be positive.
114+
87115
All arguments must have the same type.
88116

89117
### Return value
90118

91-
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`.
119+
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`.
120+
92121

93122
### Example
94123

@@ -112,8 +141,14 @@ For a complex variable \(z=(x + y i)\) with independent real \(x\) and imaginar
112141

113142
$$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}$$
114143

144+
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.
145+
115146
### Syntax
116147

148+
`result = ` [[stdlib_stats_distribution_exponential(module):cdf_exp(interface)]] `(x, loc, scale)`
149+
150+
or
151+
117152
`result = ` [[stdlib_stats_distribution_exponential(module):cdf_exp(interface)]] `(x, lambda)`
118153

119154
### Class
@@ -127,11 +162,16 @@ Elemental function
127162
`lambda`: has `intent(in)` and is a scalar of type `real` or `complex`.
128163
If `lambda` is `real`, its value must be positive. If `lambda` is `complex`, both the real and imaginary components must be positive.
129164

165+
`loc`: has `intent(in)` and is a scalar of type `real` or `complex`.
166+
167+
`scale`: has `intent(in)` and is a positive scalar of type `real` or `complex`.
168+
If `scale` is `real`, its value must be positive. If `scale` is `complex`, both the real and imaginary components must be positive.
169+
130170
All arguments must have the same type.
131171

132172
### Return value
133173

134-
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`.
174+
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`.
135175

136176
### Example
137177

‎example/stats_distribution_exponential/example_exponential_cdf.f90

Lines changed: 40 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -4,60 +4,80 @@ program example_exponential_cdf
44
rexp => rvs_exp
55

66
implicit none
7-
real, dimension(2, 3, 4) :: x, lambda
7+
real, dimension(2, 3, 4) :: x, loc, scale
88
real :: xsum
9-
complex :: scale
9+
complex :: cloc, cscale
1010
integer :: seed_put, seed_get, i
1111

1212
seed_put = 1234567
1313
call random_seed(seed_put, seed_get)
1414

15-
! standard exponential cumulative distribution at x=1.0
15+
! standard exponential cumulative distribution at x=1.0 with loc=0.0, scale=1.0
16+
print *, exp_cdf(1.0, 0.0, 1.0)
17+
! 0.632120550
18+
19+
! standard exponential cumulative distribution at x=1.0 with lambda=1.0
1620
print *, exp_cdf(1.0, 1.0)
1721
! 0.632120550
1822

1923
! cumulative distribution at x=2.0 with lambda=2
2024
print *, exp_cdf(2.0, 2.0)
2125
! 0.981684387
2226

23-
! cumulative distribution at x=2.0 with lambda=-1.0 (out of range)
24-
print *, exp_cdf(2.0, -1.0)
27+
! cumulative distribution at x=2.0 with loc=0.0 and scale=0.5 (equivalent of lambda=2)
28+
print *, exp_cdf(2.0, 0.0, 0.5)
29+
! 0.981684387
30+
31+
! cumulative distribution at x=2.5 with loc=0.5 and scale=0.5 (equivalent of lambda=2)
32+
print *, exp_cdf(2.5, 0.5, 0.5)
33+
! 0.981684387
34+
35+
! cumulative distribution at x=2.0 with loc=0.0 and scale=-1.0 (out of range)
36+
print *, exp_cdf(2.0, 0.0, -1.0)
2537
! NaN
2638

39+
! cumulative distribution at x=0.5 with loc=1.0 and scale=1.0, putting x below the minimum
40+
print *, exp_cdf(0.5, 1.0, 1.0)
41+
! 0.00000000
42+
2743
! standard exponential random variates array
28-
x = reshape(rexp(0.5, 24), [2, 3, 4])
44+
x = reshape(rexp(0.0, 2.0, 24), [2, 3, 4])
2945

3046
! a rank-3 exponential cumulative distribution
31-
lambda(:, :, :) = 0.5
32-
print *, exp_cdf(x, lambda)
47+
loc(:, :, :) = 0.0
48+
scale(:, :, :) = 2.0
49+
print *, exp_cdf(x, loc, scale)
3350
! 0.301409245 0.335173965 5.94930053E-02 0.113003314
3451
! 0.365694344 0.583515942 0.113774836 0.838585377
3552
! 0.509324908 0.127967060 0.857194781 0.893231630
3653
! 0.355383813 0.470882893 0.574203610 0.799321830
3754
! 0.546216846 0.111995399 0.801794767 0.922525287
38-
! 0.937719882 0.301136374 3.44503522E-02 0.134661376
55+
! 0.937719882 0.301136374 3.44503522E-02 0.134661376
56+
3957

40-
! cumulative distribution array where lambda<=0.0 for certain elements
41-
print *, exp_cdf([1.0, 1.0, 1.0], [1.0, 0.0, -1.0])
58+
! cumulative distribution array where scale<=0.0 for certain elements
59+
print *, exp_cdf([1.0, 1.0, 1.0], [0.0, 0.0, 0.0], [1.0, 0.0, -1.0])
4260
! 0.632120550 NaN NaN
4361

44-
! `cdf_exp` is pure and, thus, can be called concurrently
62+
! `cdf_exp` is pure and, thus, can be called concurrently
4563
xsum = 0.0
4664
do concurrent (i=1:size(x,3))
47-
xsum = xsum + sum(exp_cdf(x(:,:,i), lambda(:,:,i)))
65+
xsum = xsum + sum(exp_cdf(x(:,:,i), loc(:,:,i), scale(:,:,i)))
4866
end do
4967
print *, xsum
5068
! 11.0886612
5169

52-
! complex exponential cumulative distribution at (0.5,0.5) with real part of
53-
! lambda=0.5 and imaginary part of lambda=1.0
54-
scale = (0.5, 1.0)
55-
print *, exp_cdf((0.5, 0.5), scale)
70+
! complex exponential cumulative distribution at (0.5, 0.0, 2) with real part of
71+
! scale=2 and imaginary part of scale=1.0
72+
cloc = (0.0, 0.0)
73+
cscale = (2, 1.0)
74+
print *, exp_cdf((0.5, 0.5), cloc, cscale)
5675
! 8.70351046E-02
5776

58-
! As above, but with lambda%im < 0
59-
scale = (1.0, -2.0)
60-
print *, exp_cdf((1.5, 1.0), scale)
77+
! As above, but with scale%im < 0
78+
cloc = (0.0, 0.0)
79+
cscale = (1.0, -2.0)
80+
print *, exp_cdf((1.5, 1.0), cloc, cscale)
6181
! NaN
6282

6383
end program example_exponential_cdf

‎example/stats_distribution_exponential/example_exponential_pdf.f90

Lines changed: 37 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -4,59 +4,74 @@ program example_exponential_pdf
44
rexp => rvs_exp
55

66
implicit none
7-
real, dimension(2, 3, 4) :: x, lambda
7+
real, dimension(2, 3, 4) :: x, loc, scale
88
real :: xsum
9-
complex :: scale
9+
complex :: cloc, cscale
1010
integer :: seed_put, seed_get, i
1111

1212
seed_put = 1234567
1313
call random_seed(seed_put, seed_get)
1414

15-
! probability density at x=1.0 in standard exponential
15+
! probability density at x=1.0 with loc=0 and scale=1.0
16+
print *, exp_pdf(1.0, 0.0, 1.0)
17+
! 0.367879450
18+
19+
! probability density at x=1.0 with lambda=1.0
1620
print *, exp_pdf(1.0, 1.0)
1721
! 0.367879450
1822

1923
! probability density at x=2.0 with lambda=2.0
20-
print *, exp_pdf(2.0, 2.0)
24+
print *, exp_pdf(2.0, 2.0)
25+
! 3.66312787E-02
26+
27+
! probability density at x=2.0 with loc=0.0 and scale=0.5 (lambda=2.0)
28+
print *, exp_pdf(2.0, 0.0, 0.5)
29+
! 3.66312787E-02
30+
31+
! probability density at x=1.5 with loc=0.5 and scale=0.5 (lambda=2.0)
32+
print *, exp_pdf(2.5, 0.5, 0.5)
2133
! 3.66312787E-02
2234

23-
! probability density at x=2.0 with lambda=-1.0 (out of range)
24-
print *, exp_pdf(2.0, -1.0)
35+
! probability density at x=2.0 with loc=0.0 and scale=-1.0 (out of range)
36+
print *, exp_pdf(2.0, 0.0, -1.0)
2537
! NaN
2638

27-
! standard exponential random variates array
28-
x = reshape(rexp(0.5, 24), [2, 3, 4])
39+
! standard exponential random variates array
40+
x = reshape(rexp(0.0, 2.0, 24), [2, 3, 4])
2941

3042
! a rank-3 exponential probability density
31-
lambda(:, :, :) = 0.5
32-
print *, exp_pdf(x, lambda)
43+
loc(:, :, :) = 0.0
44+
scale(:, :, :) = 2.0
45+
print *, exp_pdf(x, loc, scale)
3346
! 0.349295378 0.332413018 0.470253497 0.443498343 0.317152828
3447
! 0.208242029 0.443112582 8.07073265E-02 0.245337561 0.436016470
3548
! 7.14025944E-02 5.33841923E-02 0.322308093 0.264558554 0.212898195
3649
! 0.100339092 0.226891592 0.444002301 9.91026312E-02 3.87373678E-02
37-
! 3.11400592E-02 0.349431813 0.482774824 0.432669312
50+
! 3.11400592E-02 0.349431813 0.482774824 0.432669312
3851

39-
! probability density array where lambda<=0.0 for certain elements
40-
print *, exp_pdf([1.0, 1.0, 1.0], [1.0, 0.0, -1.0])
52+
! probability density array where scale<=0.0 for certain elements (loc = 0.0)
53+
print *, exp_pdf([1.0, 1.0, 1.0], [0.0, 0.0, 0.0], [1.0, 0.0, -1.0])
4154
! 0.367879450 NaN NaN
4255

43-
! `pdf_exp` is pure and, thus, can be called concurrently
56+
! `pdf_exp` is pure and, thus, can be called concurrently
4457
xsum = 0.0
4558
do concurrent (i=1:size(x,3))
46-
xsum = xsum + sum(exp_pdf(x(:,:,i), lambda(:,:,i)))
59+
xsum = xsum + sum(exp_pdf(x(:,:,i), loc(:,:,i), scale(:,:,i)))
4760
end do
4861
print *, xsum
4962
! 6.45566940
5063

51-
! complex exponential probability density function at (1.5,1.0) with real part
52-
! of lambda=1.0 and imaginary part of lambda=2.0
53-
scale = (1.0, 2.)
54-
print *, exp_pdf((1.5, 1.0), scale)
64+
! complex exponential probability density function at (1.5, 0.0, 1.0) with real part
65+
! of scale=1.0 and imaginary part of scale=0.5
66+
cloc = (0.0, 0.0)
67+
cscale = (1.0, 0.5)
68+
print *, exp_pdf((1.5, 1.0), cloc, cscale)
5569
! 6.03947677E-02
5670

57-
! As above, but with lambda%re < 0
58-
scale = (-1.0, 2.)
59-
print *, exp_pdf((1.5, 1.0), scale)
71+
! As above, but with scale%re < 0
72+
cloc = (0.0, 0.0)
73+
cscale = (-1.0, 2.0)
74+
print *, exp_pdf((1.5, 1.0), cloc, cscale)
6075
! NaN
6176

6277
end program example_exponential_pdf

‎example/stats_distribution_exponential/example_exponential_rvs.f90

Lines changed: 28 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -3,30 +3,38 @@ program example_exponential_rvs
33
use stdlib_stats_distribution_exponential, only: rexp => rvs_exp
44

55
implicit none
6-
complex :: scale
6+
complex :: cloc, cscale
77
integer :: seed_put, seed_get
88

99
seed_put = 1234567
1010
call random_seed(seed_put, seed_get)
1111

12-
print *, rexp() !single standard exponential random variate
13-
14-
! 0.358690143
15-
16-
print *, rexp(2.0) !exponential random variate with lambda=2.0
17-
18-
! 0.816459715
19-
20-
print *, rexp(0.3, 10) !an array of 10 variates with lambda=0.3
21-
22-
! 1.84008647E-02 3.59742008E-02 0.136567295 0.262772143 3.62352766E-02
23-
! 0.547133625 0.213591918 4.10784185E-02 0.583882213 0.671128035
24-
25-
scale = (2.0, 0.7)
26-
print *, rexp(scale)
27-
!single complex exponential random variate with real part of lambda=2.0;
28-
!imagainary part of lambda=0.7
29-
30-
! (1.41435969,4.081114382E-02)
12+
! single standard exponential random variate
13+
print *, rexp()
14+
! 0.358690143
15+
16+
! exponential random variate with loc=0 and scale=0.5 (lambda=2)
17+
print *, rexp(0.0, 0.5)
18+
! 0.122672431
19+
20+
! exponential random variate with lambda=2
21+
print *, rexp(2.0)
22+
! 0.204114929
23+
24+
! exponential random variate with loc=0.6 and scale=0.2 (lambda=5)
25+
print *, rexp(0.6, 0.2)
26+
! 0.681645989
27+
28+
! an array of 10 variates with loc=0.0 and scale=3.0 (lambda=1/3)
29+
print *, rexp(0.0, 3.0, 10)
30+
! 1.36567295 2.62772131 0.362352759 5.47133636 2.13591909
31+
! 0.410784155 5.83882189 6.71128035 1.31730068 1.90963650
32+
33+
! single complex exponential random variate with real part of scale=0.5 (lambda=2.0);
34+
! imagainary part of scale=1.6 (lambda=0.625)
35+
cloc = (0.0, 0.0)
36+
cscale = (0.5, 1.6)
37+
print *, rexp(cloc, cscale)
38+
! (0.426896989,2.56968451)
3139

3240
end program example_exponential_rvs

‎src/stdlib_stats_distribution_exponential.fypp

Lines changed: 200 additions & 39 deletions
Original file line numberDiff line numberDiff line change
@@ -9,7 +9,7 @@ module stdlib_stats_distribution_exponential
99
implicit none
1010
private
1111

12-
integer :: ke(0:255)
12+
integer :: ke(0:255)
1313
real(dp) :: we(0:255), fe(0:255)
1414
logical :: zig_exp_initialized = .false.
1515

@@ -26,14 +26,26 @@ module stdlib_stats_distribution_exponential
2626
!! ([Specification](../page/specs/stdlib_stats_distribution_exponential.html#
2727
!! rvs_exp-exponential-distribution-random-variates))
2828
!!
29-
module procedure rvs_exp_0_rsp !0 dummy variable
29+
module procedure rvs_exp_0_rsp !0 dummy variable
30+
31+
! new interfaces using loc and scale
3032

3133
#:for k1, t1 in RC_KINDS_TYPES
32-
module procedure rvs_exp_${t1[0]}$${k1}$ !1 dummy variable
34+
module procedure rvs_exp_${t1[0]}$${k1}$ !1 dummy variable
3335
#:endfor
3436

3537
#:for k1, t1 in RC_KINDS_TYPES
36-
module procedure rvs_exp_array_${t1[0]}$${k1}$ !2 dummy variables
38+
module procedure rvs_exp_array_${t1[0]}$${k1}$ !2 dummy variables
39+
#:endfor
40+
41+
! original interfaces using lambda
42+
43+
#:for k1, t1 in RC_KINDS_TYPES
44+
module procedure rvs_exp_lambda_${t1[0]}$${k1}$ !1 dummy variable
45+
#:endfor
46+
47+
#:for k1, t1 in RC_KINDS_TYPES
48+
module procedure rvs_exp_array_lambda_${t1[0]}$${k1}$ !2 dummy variables
3749
#:endfor
3850
end interface rvs_exp
3951

@@ -46,9 +58,17 @@ module stdlib_stats_distribution_exponential
4658
!! ([Specification](../page/specs/stdlib_stats_distribution_exponential.html#
4759
!! pdf_exp-exponential-distribution-probability-density-function))
4860
!!
61+
! new interfaces using loc and scale
62+
4963
#:for k1, t1 in RC_KINDS_TYPES
5064
module procedure pdf_exp_${t1[0]}$${k1}$
5165
#:endfor
66+
67+
! original interfaces using lambda
68+
69+
#:for k1, t1 in RC_KINDS_TYPES
70+
module procedure pdf_exp_lambda_${t1[0]}$${k1}$
71+
#:endfor
5272
end interface pdf_exp
5373

5474

@@ -60,9 +80,17 @@ module stdlib_stats_distribution_exponential
6080
!! ([Specification](../page/specs/stdlib_stats_distribution_exponential.html#
6181
!! cdf_exp-exponential-distribution-cumulative-distribution-function))
6282
!!
83+
! new interfaces using loc and scale
84+
6385
#:for k1, t1 in RC_KINDS_TYPES
6486
module procedure cdf_exp_${t1[0]}$${k1}$
6587
#:endfor
88+
89+
! original interfaces using lambda
90+
91+
#:for k1, t1 in RC_KINDS_TYPES
92+
module procedure cdf_exp_lambda_${t1[0]}$${k1}$
93+
#:endfor
6694
end interface cdf_exp
6795

6896

@@ -114,7 +142,7 @@ contains
114142
#:for k1, t1 in REAL_KINDS_TYPES
115143
impure function rvs_exp_0_${t1[0]}$${k1}$( ) result(res)
116144
!
117-
! Standard exponential random variate (lambda=1)
145+
! Standard exponential random variate (lambda=1; scale=1/lambda=1)
118146
!
119147
${t1}$ :: res, x
120148
${t1}$, parameter :: r = 7.69711747013104972_${k1}$
@@ -153,18 +181,18 @@ contains
153181

154182

155183
#:for k1, t1 in REAL_KINDS_TYPES
156-
impure elemental function rvs_exp_${t1[0]}$${k1}$(lambda) result(res)
184+
impure elemental function rvs_exp_${t1[0]}$${k1}$(loc, scale) result(res)
157185
!
158186
! Exponential distributed random variate
159187
!
160-
${t1}$, intent(in) :: lambda
188+
${t1}$, intent(in) :: loc, scale
161189
${t1}$ :: res
162190

163-
if (lambda <= 0._${k1}$) then
191+
if (scale <= 0._${k1}$) then
164192
res = ieee_value(1._${k1}$, ieee_quiet_nan)
165193
else
166194
res = rvs_exp_0_${t1[0]}$${k1}$( )
167-
res = res / lambda
195+
res = res * scale + loc
168196
end if
169197
end function rvs_exp_${t1[0]}$${k1}$
170198

@@ -174,13 +202,13 @@ contains
174202

175203

176204
#:for k1, t1 in CMPLX_KINDS_TYPES
177-
impure elemental function rvs_exp_${t1[0]}$${k1}$(lambda) result(res)
178-
${t1}$, intent(in) :: lambda
205+
impure elemental function rvs_exp_${t1[0]}$${k1}$(loc, scale) result(res)
206+
${t1}$, intent(in) :: loc, scale
179207
${t1}$ :: res
180208
real(${k1}$) :: tr, ti
181209

182-
tr = rvs_exp_r${k1}$(lambda % re)
183-
ti = rvs_exp_r${k1}$(lambda % im)
210+
tr = rvs_exp_r${k1}$(loc%re, scale%re)
211+
ti = rvs_exp_r${k1}$(loc%im, scale%im)
184212
res = cmplx(tr, ti, kind=${k1}$)
185213
end function rvs_exp_${t1[0]}$${k1}$
186214

@@ -190,14 +218,14 @@ contains
190218

191219

192220
#:for k1, t1 in REAL_KINDS_TYPES
193-
impure function rvs_exp_array_${t1[0]}$${k1}$(lambda, array_size) result(res)
194-
${t1}$, intent(in) :: lambda
221+
impure function rvs_exp_array_${t1[0]}$${k1}$(loc, scale, array_size) result(res)
222+
${t1}$, intent(in) :: loc, scale
195223
integer, intent(in) :: array_size
196224
${t1}$ :: res(array_size), x, re
197225
${t1}$, parameter :: r = 7.69711747013104972_${k1}$
198226
integer :: jz, iz, i
199227

200-
if (lambda <= 0._${k1}$) then
228+
if (scale <= 0._${k1}$) then
201229
res = ieee_value(1._${k1}$, ieee_quiet_nan)
202230
return
203231
end if
@@ -228,7 +256,7 @@ contains
228256
end if
229257
end do L1
230258
endif
231-
res(i) = re / lambda
259+
res(i) = re * scale + loc
232260
end do
233261
end function rvs_exp_array_${t1[0]}$${k1}$
234262

@@ -238,16 +266,16 @@ contains
238266

239267

240268
#:for k1, t1 in CMPLX_KINDS_TYPES
241-
impure function rvs_exp_array_${t1[0]}$${k1}$(lambda, array_size) result(res)
242-
${t1}$, intent(in) :: lambda
269+
impure function rvs_exp_array_${t1[0]}$${k1}$(loc, scale, array_size) result(res)
270+
${t1}$, intent(in) :: loc, scale
243271
integer, intent(in) :: array_size
244272
${t1}$ :: res(array_size)
245273
integer :: i
246274
real(${k1}$) :: tr, ti
247275

248276
do i = 1, array_size
249-
tr = rvs_exp_r${k1}$(lambda % re)
250-
ti = rvs_exp_r${k1}$(lambda % im)
277+
tr = rvs_exp_r${k1}$(loc%re, scale%re)
278+
ti = rvs_exp_r${k1}$(loc%im, scale%im)
251279
res(i) = cmplx(tr, ti, kind=${k1}$)
252280
end do
253281
end function rvs_exp_array_${t1[0]}$${k1}$
@@ -258,19 +286,19 @@ contains
258286

259287

260288
#:for k1, t1 in REAL_KINDS_TYPES
261-
elemental function pdf_exp_${t1[0]}$${k1}$(x, lambda) result(res)
289+
elemental function pdf_exp_${t1[0]}$${k1}$(x, loc, scale) result(res)
262290
!
263291
! Exponential Distribution Probability Density Function
264292
!
265-
${t1}$, intent(in) :: x, lambda
293+
${t1}$, intent(in) :: x, loc, scale
266294
real(${k1}$) :: res
267295

268-
if (lambda <= 0._${k1}$) then
296+
if (scale <= 0._${k1}$) then
269297
res = ieee_value(1._${k1}$, ieee_quiet_nan)
270-
else if (x < 0._${k1}$) then
298+
else if (x < loc) then
271299
res = 0._${k1}$
272300
else
273-
res = exp(- x * lambda) * lambda
301+
res = exp(- (x - loc) / scale) / scale
274302
end if
275303
end function pdf_exp_${t1[0]}$${k1}$
276304

@@ -280,12 +308,12 @@ contains
280308

281309

282310
#:for k1, t1 in CMPLX_KINDS_TYPES
283-
elemental function pdf_exp_${t1[0]}$${k1}$(x, lambda) result(res)
284-
${t1}$, intent(in) :: x, lambda
311+
elemental function pdf_exp_${t1[0]}$${k1}$(x, loc, scale) result(res)
312+
${t1}$, intent(in) :: x, loc, scale
285313
real(${k1}$) :: res
286314

287-
res = pdf_exp_r${k1}$(x % re, lambda % re)
288-
res = res * pdf_exp_r${k1}$(x % im, lambda % im)
315+
res = pdf_exp_r${k1}$(x%re, loc%re, scale%re)
316+
res = res * pdf_exp_r${k1}$(x%im, loc%im, scale%im)
289317
end function pdf_exp_${t1[0]}$${k1}$
290318

291319
#:endfor
@@ -294,19 +322,19 @@ contains
294322

295323

296324
#:for k1, t1 in REAL_KINDS_TYPES
297-
elemental function cdf_exp_${t1[0]}$${k1}$(x, lambda) result(res)
325+
elemental function cdf_exp_${t1[0]}$${k1}$(x, loc, scale) result(res)
298326
!
299327
! Exponential Distribution Cumulative Distribution Function
300328
!
301-
${t1}$, intent(in) :: x, lambda
329+
${t1}$, intent(in) :: x, loc, scale
302330
real(${k1}$) :: res
303331

304-
if (lambda <= 0._${k1}$) then
332+
if (scale <= 0._${k1}$) then
305333
res = ieee_value(1._${k1}$, ieee_quiet_nan)
306-
else if (x < 0._${k1}$) then
334+
else if (x < loc) then
307335
res = 0._${k1}$
308336
else
309-
res = 1.0_${k1}$ - exp(- x * lambda)
337+
res = 1.0_${k1}$ - exp(- (x - loc) / scale)
310338
end if
311339
end function cdf_exp_${t1[0]}$${k1}$
312340

@@ -316,14 +344,147 @@ contains
316344

317345

318346
#:for k1, t1 in CMPLX_KINDS_TYPES
319-
elemental function cdf_exp_${t1[0]}$${k1}$(x, lambda) result(res)
320-
${t1}$, intent(in) :: x, lambda
347+
elemental function cdf_exp_${t1[0]}$${k1}$(x, loc, scale) result(res)
348+
${t1}$, intent(in) :: x, loc, scale
321349
real(${k1}$) :: res
322350

323-
res = cdf_exp_r${k1}$(x % re, lambda % re)
324-
res = res * cdf_exp_r${k1}$(x % im, lambda % im)
351+
res = cdf_exp_r${k1}$(x%re, loc%re, scale%re)
352+
res = res * cdf_exp_r${k1}$(x%im, loc%im, scale%im)
325353
end function cdf_exp_${t1[0]}$${k1}$
326354

327355
#:endfor
328356

357+
358+
359+
360+
!
361+
! below: wrapper functions for interfaces using lambda (for backward compatibility)
362+
!
363+
364+
365+
366+
367+
#:for k1, t1 in REAL_KINDS_TYPES
368+
impure elemental function rvs_exp_lambda_${t1[0]}$${k1}$(lambda) result(res)
369+
!
370+
! Exponential distributed random variate
371+
!
372+
${t1}$, intent(in) :: lambda
373+
${t1}$ :: res
374+
375+
res = rvs_exp_${t1[0]}$${k1}$(0._${k1}$, 1.0_${k1}$/lambda)
376+
end function rvs_exp_lambda_${t1[0]}$${k1}$
377+
378+
#:endfor
379+
380+
381+
382+
383+
#:for k1, t1 in CMPLX_KINDS_TYPES
384+
impure elemental function rvs_exp_lambda_${t1[0]}$${k1}$(lambda) result(res)
385+
${t1}$, intent(in) :: lambda
386+
${t1}$ :: res
387+
real(${k1}$) :: tr, ti
388+
389+
tr = rvs_exp_r${k1}$(0._${k1}$, 1.0_${k1}$/lambda%re)
390+
ti = rvs_exp_r${k1}$(0._${k1}$, 1.0_${k1}$/lambda%im)
391+
res = cmplx(tr, ti, kind=${k1}$)
392+
end function rvs_exp_lambda_${t1[0]}$${k1}$
393+
394+
#:endfor
395+
396+
397+
398+
399+
#:for k1, t1 in REAL_KINDS_TYPES
400+
impure function rvs_exp_array_lambda_${t1[0]}$${k1}$(lambda, array_size) result(res)
401+
${t1}$, intent(in) :: lambda
402+
integer, intent(in) :: array_size
403+
${t1}$ :: res(array_size)
404+
405+
res = rvs_exp_array_${t1[0]}$${k1}$(0._${k1}$, 1.0_${k1}$/lambda, array_size)
406+
end function rvs_exp_array_lambda_${t1[0]}$${k1}$
407+
408+
#:endfor
409+
410+
411+
412+
413+
#:for k1, t1 in CMPLX_KINDS_TYPES
414+
impure function rvs_exp_array_lambda_${t1[0]}$${k1}$(lambda, array_size) result(res)
415+
${t1}$, intent(in) :: lambda
416+
integer, intent(in) :: array_size
417+
${t1}$ :: res(array_size)
418+
integer :: i
419+
real(${k1}$) :: tr, ti
420+
421+
do i = 1, array_size
422+
tr = rvs_exp_r${k1}$(0._${k1}$, 1.0_${k1}$/lambda%re)
423+
ti = rvs_exp_r${k1}$(0._${k1}$, 1.0_${k1}$/lambda%im)
424+
res(i) = cmplx(tr, ti, kind=${k1}$)
425+
end do
426+
end function rvs_exp_array_lambda_${t1[0]}$${k1}$
427+
428+
#:endfor
429+
430+
431+
432+
433+
#:for k1, t1 in REAL_KINDS_TYPES
434+
elemental function pdf_exp_lambda_${t1[0]}$${k1}$(x, lambda) result(res)
435+
!
436+
! Exponential Distribution Probability Density Function
437+
!
438+
${t1}$, intent(in) :: x, lambda
439+
real(${k1}$) :: res
440+
441+
res = pdf_exp_${t1[0]}$${k1}$(x, 0._${k1}$, 1.0_${k1}$/lambda)
442+
end function pdf_exp_lambda_${t1[0]}$${k1}$
443+
444+
#:endfor
445+
446+
447+
448+
449+
#:for k1, t1 in CMPLX_KINDS_TYPES
450+
elemental function pdf_exp_lambda_${t1[0]}$${k1}$(x, lambda) result(res)
451+
${t1}$, intent(in) :: x, lambda
452+
real(${k1}$) :: res
453+
454+
res = pdf_exp_r${k1}$(x%re, 0._${k1}$, 1.0_${k1}$/lambda%re)
455+
res = res * pdf_exp_r${k1}$(x%im, 0._${k1}$, 1.0_${k1}$/lambda%im)
456+
end function pdf_exp_lambda_${t1[0]}$${k1}$
457+
458+
#:endfor
459+
460+
461+
462+
463+
#:for k1, t1 in REAL_KINDS_TYPES
464+
elemental function cdf_exp_lambda_${t1[0]}$${k1}$(x, lambda) result(res)
465+
!
466+
! Exponential Distribution Cumulative Distribution Function
467+
!
468+
${t1}$, intent(in) :: x, lambda
469+
real(${k1}$) :: res
470+
471+
res = cdf_exp_${t1[0]}$${k1}$(x, 0._${k1}$, 1.0_${k1}$/lambda)
472+
end function cdf_exp_lambda_${t1[0]}$${k1}$
473+
474+
#:endfor
475+
476+
477+
478+
479+
#:for k1, t1 in CMPLX_KINDS_TYPES
480+
elemental function cdf_exp_lambda_${t1[0]}$${k1}$(x, lambda) result(res)
481+
${t1}$, intent(in) :: x, lambda
482+
real(${k1}$) :: res
483+
484+
res = cdf_exp_r${k1}$(x%re, 0._${k1}$, 1.0_${k1}$/lambda%re)
485+
res = res * cdf_exp_r${k1}$(x%im, 0._${k1}$, 1.0_${k1}$/lambda%im)
486+
end function cdf_exp_lambda_${t1[0]}$${k1}$
487+
488+
#:endfor
489+
329490
end module stdlib_stats_distribution_exponential

‎test/stats/test_distribution_exponential.fypp

Lines changed: 90 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -51,6 +51,8 @@ contains
5151

5252
print *, ""
5353
print *, "Test exponential random generator with chi-squared"
54+
55+
! using interface for lambda
5456
freq = 0
5557
do i = 1, num
5658
j = 1000 * (1 - exp(- expon_rvs(1.0)))
@@ -66,13 +68,31 @@ contains
6668
write(*,*) "Chi-squared for exponential random generator is : ", chisq
6769
call check((chisq < 1143.9), &
6870
msg="exponential randomness failed chi-squared test", warn=warn)
71+
72+
! using interface for loc and scale
73+
freq = 0
74+
do i = 1, num
75+
j = 1000 * (1 - exp(- expon_rvs(0.0, 1.0)))
76+
freq(j) = freq(j) + 1
77+
end do
78+
chisq = 0.0_dp
79+
expct = num / array_size
80+
do i = 0, array_size - 1
81+
chisq = chisq + (freq(i) - expct) ** 2 / expct
82+
end do
83+
write(*,*) "The critical values for chi-squared with 1000 d. of f. is" &
84+
//" 1143.92"
85+
write(*,*) "Chi-squared for exponential random generator is : ", chisq
86+
call check((chisq < 1143.9), &
87+
msg="exponential randomness failed chi-squared test", warn=warn)
88+
6989
end subroutine test_exponential_random_generator
7090

7191

7292

7393
#:for k1, t1 in RC_KINDS_TYPES
7494
subroutine test_expon_rvs_${t1[0]}$${k1}$
75-
${t1}$ :: res(10), scale
95+
${t1}$ :: res(10), lambda, loc, scale
7696
integer, parameter :: k = 5
7797
integer :: i
7898
integer :: seed, get
@@ -116,20 +136,35 @@ contains
116136

117137
print *, "Test exponential_distribution_rvs_${t1[0]}$${k1}$"
118138
seed = 593742186
119-
call random_seed(seed, get)
120139

140+
! set args
121141
#:if t1[0] == "r"
122142
#! for real type
123-
scale = 1.5_${k1}$
143+
lambda = 1.5_${k1}$
144+
loc = 0._${k1}$
145+
scale = 1.0_${k1}$/lambda
124146
#:else
125147
#! for complex type
126-
scale = (0.7_${k1}$, 1.3_${k1}$)
148+
lambda = (0.7_${k1}$, 1.3_${k1}$)
149+
loc = (0._${k1}$, 0._${k1}$)
150+
scale = cmplx(1.0_${k1}$/lambda%re, 1.0_${k1}$/lambda%im, kind=${k1}$)
127151
#:endif
128152

153+
! tests using interface for lambda
154+
call random_seed(seed, get)
155+
do i = 1, k
156+
res(i) = expon_rvs(lambda) ! 1 dummy
157+
end do
158+
res(6:10) = expon_rvs(lambda, k) ! 2 dummies
159+
call check(all(abs(res - ans) < ${k1}$tol), &
160+
msg="exponential_distribution_rvs_${t1[0]}$${k1}$ failed", warn=warn)
161+
162+
! tests using interface for loc and scale
163+
call random_seed(seed, get)
129164
do i = 1, k
130-
res(i) = expon_rvs(scale) ! 1 dummy
165+
res(i) = expon_rvs(loc, scale)
131166
end do
132-
res(6:10) = expon_rvs(scale, k) ! 2 dummies
167+
res(6:10) = expon_rvs(loc, scale, k)
133168
call check(all(abs(res - ans) < ${k1}$tol), &
134169
msg="exponential_distribution_rvs_${t1[0]}$${k1}$ failed", warn=warn)
135170
end subroutine test_expon_rvs_${t1[0]}$${k1}$
@@ -142,7 +177,7 @@ contains
142177
#:for k1, t1 in RC_KINDS_TYPES
143178
subroutine test_expon_pdf_${t1[0]}$${k1}$
144179

145-
${t1}$ :: x1, x2(3,4), scale
180+
${t1}$ :: x1, x2(3,4), lambda, loc, scale
146181
integer :: seed, get
147182
real(${k1}$) :: res(3,5)
148183
#:if t1[0] == "r"
@@ -185,21 +220,38 @@ contains
185220

186221
print *, "Test exponential_distribution_pdf_${t1[0]}$${k1}$"
187222
seed = 123987654
188-
call random_seed(seed, get)
223+
224+
! set args
189225
#:if t1[0] == "r"
190226
#! for real type
191-
scale = 1.5_${k1}$
227+
lambda = 1.5_${k1}$
228+
loc = 0._${k1}$
229+
scale = 1.0_${k1}$/lambda
192230
#:else
193231
#! for complex type
194-
scale = (0.3_${k1}$, 1.6_${k1}$)
232+
lambda = (0.3_${k1}$, 1.6_${k1}$)
233+
loc = (0._${k1}$, 0._${k1}$)
234+
scale = cmplx(1.0_${k1}$/lambda%re, 1.0_${k1}$/lambda%im, kind=${k1}$)
195235
#:endif
196236

197-
x1 = expon_rvs(scale)
198-
x2 = reshape(expon_rvs(scale, 12), [3,4])
199-
res(:,1) = expon_pdf(x1, scale)
200-
res(:, 2:5) = expon_pdf(x2, scale)
237+
! tests using interface for lambda
238+
call random_seed(seed, get)
239+
x1 = expon_rvs(lambda)
240+
x2 = reshape(expon_rvs(lambda, 12), [3,4])
241+
res(:,1) = expon_pdf(x1, lambda)
242+
res(:, 2:5) = expon_pdf(x2, lambda)
243+
call check(all(abs(res - reshape(ans, [3,5])) < ${k1}$tol), &
244+
msg="exponential_distribution_pdf_${t1[0]}$${k1}$ failed", warn=warn)
245+
246+
! tests using interface for loc and scale
247+
call random_seed(seed, get)
248+
x1 = expon_rvs(loc, scale)
249+
x2 = reshape(expon_rvs(loc, scale, 12), [3,4])
250+
res(:,1) = expon_pdf(x1, loc, scale)
251+
res(:, 2:5) = expon_pdf(x2, loc, scale)
201252
call check(all(abs(res - reshape(ans, [3,5])) < ${k1}$tol), &
202253
msg="exponential_distribution_pdf_${t1[0]}$${k1}$ failed", warn=warn)
254+
203255
end subroutine test_expon_pdf_${t1[0]}$${k1}$
204256

205257
#:endfor
@@ -210,7 +262,7 @@ contains
210262
#:for k1, t1 in RC_KINDS_TYPES
211263
subroutine test_expon_cdf_${t1[0]}$${k1}$
212264

213-
${t1}$ :: x1, x2(3,4), scale
265+
${t1}$ :: x1, x2(3,4), lambda, loc, scale
214266
integer :: seed, get
215267
real(${k1}$) :: res(3,5)
216268
#:if t1[0] == "r"
@@ -252,21 +304,37 @@ contains
252304

253305
print *, "Test exponential_distribution_cdf_${t1[0]}$${k1}$"
254306
seed = 621957438
255-
call random_seed(seed, get)
256307

308+
! set args
257309
#:if t1[0] == "r"
258310
#! for real type
259-
scale = 2.0_${k1}$
311+
lambda = 2.0_${k1}$
312+
loc = 0._${k1}$
313+
scale = 1.0_${k1}$/lambda
260314
#:else
261-
scale = (1.3_${k1}$, 2.1_${k1}$)
315+
lambda = (1.3_${k1}$, 2.1_${k1}$)
316+
loc = (0._${k1}$, 0._${k1}$)
317+
scale = cmplx(1.0_${k1}$/lambda%re, 1.0_${k1}$/lambda%im, kind=${k1}$)
262318
#:endif
263319

264-
x1 = expon_rvs(scale)
265-
x2 = reshape(expon_rvs(scale, 12), [3,4])
266-
res(:,1) = expon_cdf(x1, scale)
267-
res(:, 2:5) = expon_cdf(x2, scale)
320+
! tests using interface for lambda
321+
call random_seed(seed, get)
322+
x1 = expon_rvs(lambda)
323+
x2 = reshape(expon_rvs(lambda, 12), [3,4])
324+
res(:,1) = expon_cdf(x1, lambda)
325+
res(:, 2:5) = expon_cdf(x2, lambda)
326+
call check(all(abs(res - reshape(ans,[3,5])) < ${k1}$tol), &
327+
msg="exponential_distribution_cdf_${t1[0]}$${k1}$ failed", warn=warn)
328+
329+
! tests using interface for loc and scale
330+
call random_seed(seed, get)
331+
x1 = expon_rvs(loc, scale)
332+
x2 = reshape(expon_rvs(loc, scale, 12), [3,4])
333+
res(:,1) = expon_cdf(x1, loc, scale)
334+
res(:, 2:5) = expon_cdf(x2, loc, scale)
268335
call check(all(abs(res - reshape(ans,[3,5])) < ${k1}$tol), &
269336
msg="exponential_distribution_cdf_${t1[0]}$${k1}$ failed", warn=warn)
337+
270338
end subroutine test_expon_cdf_${t1[0]}$${k1}$
271339

272340
#:endfor

0 commit comments

Comments
 (0)
Please sign in to comment.