Skip to content

modified stats exponential distribution procedures to use loc and scale #991

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 7 commits into
base: master
Choose a base branch
from
Open

modified stats exponential distribution procedures to use loc and scale #991

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
43f88eb
modified parameterisation for stats exp dist. procedures
sebastian-mutz May 8, 2025
e7984cb
reintroduced old interface through wrappers and updated specs
sebastian-mutz May 13, 2025
ba239c6
Trigger mergeability check
sebastian-mutz May 13, 2025
9fe477b
Update doc/specs/stdlib_stats_distribution_exponential.md
perazz May 16, 2025
3c5114b
Update doc/specs/stdlib_stats_distribution_exponential.md
perazz May 16, 2025
3ae8b58
Update doc/specs/stdlib_stats_distribution_exponential.md
perazz May 16, 2025
9b99b53
new API prioritised
sebastian-mutz May 17, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
54 changes: 47 additions & 7 deletions doc/specs/stdlib_stats_distribution_exponential.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,20 +15,29 @@ 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]

### Syntax

`result = ` [[stdlib_stats_distribution_exponential(module):rvs_exp(interface)]] `([loc, scale] [[, array_size]])`

or

`result = ` [[stdlib_stats_distribution_exponential(module):rvs_exp(interface)]] `([lambda] [[, array_size]])`

### Class
Expand All @@ -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
Expand All @@ -69,8 +86,14 @@ 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, loc, scale)`

or

`result = ` [[stdlib_stats_distribution_exponential(module):pdf_exp(interface)]] `(x, lambda)`

### Class
Expand All @@ -84,11 +107,17 @@ 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`.
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

Expand All @@ -112,8 +141,14 @@ 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}$$

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

`result = ` [[stdlib_stats_distribution_exponential(module):cdf_exp(interface)]] `(x, loc, scale)`

or

`result = ` [[stdlib_stats_distribution_exponential(module):cdf_exp(interface)]] `(x, lambda)`

### Class
Expand All @@ -127,11 +162,16 @@ 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`.
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

Expand Down
60 changes: 40 additions & 20 deletions example/stats_distribution_exponential/example_exponential_cdf.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,60 +4,80 @@ 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
! 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

! 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.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.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
Expand Down
59 changes: 37 additions & 22 deletions example/stats_distribution_exponential/example_exponential_pdf.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,59 +4,74 @@ 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
! 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)
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

! 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 lambda=-1.0 (out of range)
print *, exp_pdf(2.0, -1.0)
! 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
48 changes: 28 additions & 20 deletions example/stats_distribution_exponential/example_exponential_rvs.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,30 +3,38 @@ 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

! 1.84008647E-02 3.59742008E-02 0.136567295 0.262772143 3.62352766E-02
! 0.547133625 0.213591918 4.10784185E-02 0.583882213 0.671128035

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

! (1.41435969,4.081114382E-02)
! single standard exponential random variate
print *, rexp()
! 0.358690143

! 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

! 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)
! (0.426896989,2.56968451)

end program example_exponential_rvs
Loading
Loading