Skip to content
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

Extension of PR#679 to stdlib_stats_distribution_exponential #717

Merged
merged 7 commits into from
Jun 15, 2023
Merged

Extension of PR#679 to stdlib_stats_distribution_exponential #717

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
c9682d0
replace calls to error_stop by NaN, convert procedures to pure when p…
HugoMVale Jun 8, 2023
129555e
Merge branch 'fortran-lang:master' into stats-distribution-exponential
HugoMVale Jun 8, 2023
f6b2b32
improve readibility, fix equation rendering, update specs
HugoMVale Jun 8, 2023
3e458d4
Merge branch 'stats-distribution-exponential' of github.com:HugoMVale…
HugoMVale Jun 8, 2023
47dfe53
clean up + add examples with NaN and concurrent call
HugoMVale Jun 8, 2023
94ed13f
pdf and cdf should be zero if x<0
HugoMVale Jun 8, 2023
afc1844
clean up + add examples with NaN and concurrent call
HugoMVale Jun 8, 2023
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
47 changes: 26 additions & 21 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 @@ -14,34 +14,37 @@ Experimental

### Description

An exponential distribution is the distribution of time between events in a Poisson point process. The inverse scale parameter `lambda` specifies the average time between events, also called the rate of events.
An exponential distribution is the distribution of time between events in a Poisson point process. The inverse scale parameter `lambda` specifies the average time between events ($\lambda$), also called the rate of events.

Without argument the function returns a random sample from the standard exponential distribution `E(1)` with `lambda = 1`.
Without argument, the function returns a random sample from the standard exponential distribution $E(\lambda=1)$.

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

With two arguments the function returns a rank one array of exponentially distributed random variates.
With two arguments, the function returns a rank-1 array of exponentially distributed random variates.

Note: the algorithm used for generating exponetial random variates is fundamentally limited to double precision. Ref.: Marsaglia, G. & Tsang, W.W. (2000) `The ziggurat method for generating random variables', J. Statist. Software, v5(8).
@note
The algorithm used for generating exponential random variates is fundamentally limited to double precision.[^1]

### Syntax

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

### Class

Function
Elemental function

### Arguments

`lambda`: optional argument has `intent(in)` and is a scalar of type `real` or `complex`. The value of `lambda` has to be non-negative.
`lambda`: optional argument has `intent(in)` and is a scalar of type `real` or `complex`.
If `lambda` is `real`, its value must be positive. If `lambda` is `complex`, both the real and imaginary components must be positive.

`array_size`: optional argument has `intent(in)` and is a scalar of type `integer` with default kind.

### Return value

The result is a scalar or rank one array with a size of `array_size`, and has the same type of `lambda`.
The result is a scalar or rank-1 array with a size of `array_size`, and the same type as `lambda`.
If `lambda` is non-positive, the result is `NaN`.

### Example

Expand All @@ -57,15 +60,13 @@ Experimental

### Description

The probability density function (pdf) of the single real variable exponential distribution:
The probability density function (pdf) of the single real variable exponential distribution is:

$$f(x)=\begin{cases} \lambda e^{-\lambda x} &x\geqslant 0 \\\\ 0 &x< 0\end{}$$
$$f(x)=\begin{cases} \lambda e^{-\lambda x} &x\geqslant 0 \\\\ 0 &x< 0\end{cases}$$

For a complex variable (x + y i) with independent real x and imaginary y parts, the joint probability density function
is the product of the corresponding marginal pdf of real and imaginary pdf (for more details, see
"Probability and Random Processes with Applications to Signal Processing and Communications", 2nd ed., Scott L. Miller and Donald Childers, 2012, p.197):
For a complex variable $z=(x + y i)$ with independent real $x$ and imaginary $y$ parts, the joint probability density function is the product of the corresponding real and imaginary marginal pdfs:[^2]

$$f(x+\mathit{i}y)=f(x)f(y)=\begin{cases} \lambda_{x} \lambda_{y} e^{-(\lambda_{x} x + \lambda_{y} y)} &x\geqslant 0, y\geqslant 0 \\\\ 0 &otherwise\end{}$$
$$f(x+\mathit{i}y)=f(x)f(y)=\begin{cases} \lambda_{x} \lambda_{y} e^{-(\lambda_{x} x + \lambda_{y} y)} &x\geqslant 0, y\geqslant 0 \\\\ 0 &\text{otherwise}\end{cases}$$

### Syntax

Expand All @@ -80,12 +81,13 @@ Elemental function
`x`: has `intent(in)` and is a scalar of type `real` or `complex`.

`lambda`: has `intent(in)` and is a scalar of type `real` or `complex`.
If `lambda` is `real`, its value must be positive. If `lambda` is `complex`, both the real and imaginary components must be positive.

All arguments must have the same type.

### Return value

The result is a scalar or an array, with a shape conformable to arguments, and has the same type of input arguments.
The result is a scalar or an array, with a shape conformable to the arguments, and the same type as the input arguments. If `lambda` is non-positive, the result is `NaN`.

### Example

Expand All @@ -103,13 +105,11 @@ Experimental

Cumulative distribution function (cdf) of the single real variable exponential distribution:

$$F(x)=\begin{cases}1 - e^{-\lambda x} &x\geqslant 0 \\\\ 0 &x< 0\end{}$$
$$F(x)=\begin{cases}1 - e^{-\lambda x} &x\geqslant 0 \\\\ 0 &x< 0\end{cases}$$

For a complex variable (x + y i) with independent real x and imaginary y parts, the joint cumulative distribution
function is the product of corresponding marginal cdf of real and imaginary cdf (for more details, see
"Probability and Random Processes with Applications to Signal Processing and Communications", 2nd ed., Scott L. Miller and Donald Childers, 2012, p.197):
For a complex variable $z=(x + y i)$ with independent real $x$ and imaginary $y$ parts, the joint cumulative distribution function is the product of corresponding real and imaginary marginal cdfs:[^2]

$$F(x+\mathit{i}y)=F(x)F(y)=\begin{cases} (1 - e^{-\lambda_{x} x})(1 - e^{-\lambda_{y} y}) &x\geqslant 0, \;\; y\geqslant 0 \\\\ 0 &otherwise \end{}$$
$$F(x+\mathit{i}y)=F(x)F(y)=\begin{cases} (1 - e^{-\lambda_{x} x})(1 - e^{-\lambda_{y} y}) &x\geqslant 0, \;\; y\geqslant 0 \\\\ 0 & \text{otherwise} \end{cases}$$

### Syntax

Expand All @@ -124,15 +124,20 @@ Elemental function
`x`: has `intent(in)` and is a scalar of type `real` or `complex`.

`lambda`: has `intent(in)` and is a scalar of type `real` or `complex`.
If `lambda` is `real`, its value must be positive. If `lambda` is `complex`, both the real and imaginary components must be positive.

All arguments must have the same type.

### Return value

The result is a scalar or an array, with a shape conformable to arguments, and has the same type of input arguments.
The result is a scalar or an array, with a shape conformable to the arguments, and the same type as the input arguments. If `lambda` is non-positive, the result is `NaN`.

### Example

```fortran
{!example/stats_distribution_exponential/example_exponential_cdf.f90!}
```

[^1] Marsaglia, George, and Wai Wan Tsang. "The ziggurat method for generating random variables." _Journal of statistical software_ 5 (2000): 1-7.

[^2] Miller, Scott, and Donald Childers. _Probability and random processes: With applications to signal processing and communications_. Academic Press, 2012 (p. 197).
59 changes: 41 additions & 18 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,38 +4,61 @@ program example_exponential_cdf
rexp => rvs_exp

implicit none
real :: x(2, 3, 4), a(2, 3, 4)
real, dimension(2, 3, 4) :: x, lambda
real :: xsum
complex :: scale
integer :: seed_put, seed_get
integer :: seed_put, seed_get, i

seed_put = 1234567
call random_seed(seed_put, seed_get)

print *, exp_cdf(1.0, 1.0) ! a standard exponential cumulative at 1.0
! standard exponential cumulative distribution at x=1.0
print *, exp_cdf(1.0, 1.0)
! 0.632120550

! 0.632120550
! cumulative distribution at x=2.0 with lambda=2
print *, exp_cdf(2.0, 2.0)
! 0.981684387

print *, exp_cdf(2.0, 2.0) ! a cumulative at 2.0 with lambda=2

! 0.981684387
! cumulative distribution at x=2.0 with lambda=-1.0 (out of range)
print *, exp_cdf(2.0, -1.0)
! NaN

! standard exponential random variates array
x = reshape(rexp(0.5, 24), [2, 3, 4])
! standard exponential random variates array
a(:, :, :) = 0.5
print *, exp_cdf(x, a) ! a rank 3 array of standard exponential cumulative

! 8.57694745E-02 9.70223546E-02 1.52170658E-02 2.95336246E-02
! 0.107568979 0.196659625 2.97447443E-02 0.366151094 0.163051903
! 3.36527228E-02 0.385267735 0.428375721 0.103964329 0.147119939
! 0.192206264 0.330693483 0.179247737 2.92580128E-02 0.332765043
! 0.472417951 0.500440359 8.56802464E-02 8.72612000E-03 3.55126858E-02
! a rank-3 exponential cumulative distribution
lambda(:, :, :) = 0.5
print *, exp_cdf(x, lambda)
! 0.301409245 0.335173965 5.94930053E-02 0.113003314
! 0.365694344 0.583515942 0.113774836 0.838585377
! 0.509324908 0.127967060 0.857194781 0.893231630
! 0.355383813 0.470882893 0.574203610 0.799321830
! 0.546216846 0.111995399 0.801794767 0.922525287
! 0.937719882 0.301136374 3.44503522E-02 0.134661376

! cumulative distribution array where lambda<=0.0 for certain elements
print *, exp_cdf([1.0, 1.0, 1.0], [1.0, 0.0, -1.0])
! 0.632120550 NaN NaN

! `cdf_exp` is pure and, thus, can be called concurrently
xsum = 0.0
do concurrent (i=1:size(x,3))
xsum = xsum + sum(exp_cdf(x(:,:,i), lambda(:,:,i)))
end do
print *, xsum
! 11.0886612

! complex exponential cumulative distribution at (0.5,0.5) with real part of
! lambda=0.5 and imaginary part of lambda=1.0
scale = (0.5, 1.0)
print *, exp_cdf((0.5, 0.5), scale)
!complex exponential cumulative distribution at (0.5,0.5) with real part of
!lambda=0.5 and imaginary part of lambda=1.0
! 8.70351046E-02

! 8.70351046E-02
! As above, but with lambda%im < 0
scale = (1.0, -2.0)
print *, exp_cdf((1.5, 1.0), scale)
! NaN

end program example_exponential_cdf

71 changes: 47 additions & 24 deletions example/stats_distribution_exponential/example_exponential_pdf.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,39 +1,62 @@
program example_exponential_pdf
use stdlib_random, only: random_seed
use stdlib_stats_distribution_exponential, only: exp_pdf => pdf_exp, &
rexp => rvs_exp
rexp => rvs_exp

implicit none
real :: x(2, 3, 4), a(2, 3, 4)
real, dimension(2, 3, 4) :: x, lambda
real :: xsum
complex :: scale
integer :: seed_put, seed_get
integer :: seed_put, seed_get, i

seed_put = 1234567
call random_seed(seed_put, seed_get)

print *, exp_pdf(1.0, 1.0) !a probability density at 1.0 in standard expon

! 0.367879450

print *, exp_pdf(2.0, 2.0) !a probability density at 2.0 with lambda=2.0

! 3.66312787E-02

x = reshape(rexp(0.5, 24), [2, 3, 4]) ! standard expon random variates array
a(:, :, :) = 0.5
print *, exp_pdf(x, a) ! a rank 3 standard expon probability density

! 0.457115263 0.451488823 0.492391467 0.485233188 0.446215510
! 0.401670188 0.485127628 0.316924453 0.418474048 0.483173639
! 0.307366133 0.285812140 0.448017836 0.426440030 0.403896868
! 0.334653258 0.410376132 0.485370994 0.333617479 0.263791025
! 0.249779820 0.457159877 0.495636940 0.482243657

! probability density at x=1.0 in standard exponential
print *, exp_pdf(1.0, 1.0)
! 0.367879450

! probability density at x=2.0 with lambda=2.0
print *, exp_pdf(2.0, 2.0)
! 3.66312787E-02

! probability density at x=2.0 with lambda=-1.0 (out of range)
print *, exp_pdf(2.0, -1.0)
! NaN

! standard exponential random variates array
x = reshape(rexp(0.5, 24), [2, 3, 4])

! a rank-3 exponential probability density
lambda(:, :, :) = 0.5
print *, exp_pdf(x, lambda)
! 0.349295378 0.332413018 0.470253497 0.443498343 0.317152828
! 0.208242029 0.443112582 8.07073265E-02 0.245337561 0.436016470
! 7.14025944E-02 5.33841923E-02 0.322308093 0.264558554 0.212898195
! 0.100339092 0.226891592 0.444002301 9.91026312E-02 3.87373678E-02
! 3.11400592E-02 0.349431813 0.482774824 0.432669312

! probability density array where lambda<=0.0 for certain elements
print *, exp_pdf([1.0, 1.0, 1.0], [1.0, 0.0, -1.0])
! 0.367879450 NaN NaN

! `pdf_exp` is pure and, thus, can be called concurrently
xsum = 0.0
do concurrent (i=1:size(x,3))
xsum = xsum + sum(exp_pdf(x(:,:,i), lambda(:,:,i)))
end do
print *, xsum
! 6.45566940

! complex exponential probability density function at (1.5,1.0) with real part
! of lambda=1.0 and imaginary part of lambda=2.0
scale = (1.0, 2.)
print *, exp_pdf((1.5, 1.0), scale)
! a complex expon probability density function at (1.5,1.0) with real part
!of lambda=1.0 and imaginary part of lambda=2.0
! 6.03947677E-02

! 6.03947677E-02
! As above, but with lambda%re < 0
scale = (-1.0, 2.)
print *, exp_pdf((1.5, 1.0), scale)
! NaN

end program example_exponential_pdf
Loading