Adding stable distributions

Project.toml

@@ -7,6 +7,7 @@ version = "0.25.86"
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 DensityInterface = "b429d917-457f-4dbc-8f4c-0cc954292b1d"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
+Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 PDMats = "90014a1f-27ba-587c-ab20-58faa44d9150"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"

docs/src/univariate.md

@@ -431,6 +431,13 @@ SkewNormal
 plotdensity((-4, 4), SkewNormal, (0, 1, -1)) # hide
 ```
 
+```@docs
+Stable
+```
+```@example plotdensity
+plotdensity((-4, 5), Stable, (1.5, 1, 1,0)) # hide
+```
+
 ```@docs
 StudentizedRange
 SymTriangularDist

src/Distributions.jl

@@ -4,6 +4,8 @@ using StatsBase, PDMats, StatsFuns, Statistics
 using StatsFuns: logtwo, invsqrt2, invsqrt2π
 
 import QuadGK: quadgk
+import Interpolations: interpolate, Gridded, Linear
+
 import Base: size, length, convert, show, getindex, rand, vec, inv
 import Base: sum, maximum, minimum, extrema, +, -, *, ==
 import Base.Math: @horner
@@ -153,6 +155,7 @@ export
     Skellam,
     SkewNormal,
     Soliton,
+    Stable,
     StudentizedRange,
     SymTriangularDist,
     TDist,
@@ -352,7 +355,7 @@ Supported distributions:
     MvNormalKnownCov, MvTDist, NegativeBinomial, NoncentralBeta, NoncentralChisq,
     NoncentralF, NoncentralHypergeometric, NoncentralT, Normal, NormalCanon,
     NormalInverseGaussian, Pareto, PGeneralizedGaussian, Poisson, PoissonBinomial,
-    QQPair, Rayleigh, Rician, Skellam, Soliton, StudentizedRange, SymTriangularDist, TDist, TriangularDist,
+    QQPair, Rayleigh, Rician, Skellam, Soliton, Stable, StudentizedRange, SymTriangularDist, TDist, TriangularDist,
     Triweight, Truncated, TruncatedNormal, Uniform, UnivariateGMM,
     VonMises, VonMisesFisher, WalleniusNoncentralHypergeometric, Weibull,
     Wishart, ZeroMeanIsoNormal, ZeroMeanIsoNormalCanon,

src/conversion.jl

@@ -7,3 +7,6 @@ convert(::Type{Binomial}, d::Bernoulli) = Binomial(1, d.p)
 
 convert(::Type{Gamma}, d::Exponential) = Gamma(1.0, d.θ)
 convert(::Type{Gamma}, d::Erlang) = Gamma(d.α, d.θ)
+convert(::Type{Stable}, d::Normal) = Stable(2, 0, d.σ/√2, d.μ)
+convert(::Type{Stable}, d::Cauchy) = Stable(1, 0, d.σ, d.μ)
+convert(::Type{Stable}, d::Levy) = Stable(1/2, 1, d.σ, d.μ)

src/convolution.jl

@@ -4,9 +4,7 @@
 Convolve two distributions and return the distribution corresponding to the sum of
 independent random variables drawn from the underlying distributions.
 
-Currently, the function is only defined in cases where the convolution has a closed form.
-More precisely, the function is defined if the distributions of `d1` and `d2` are the same
-and one of
+Currently, the function is only defined in cases where the convolution has a closed form and is one of
 * [`Bernoulli`](@ref)
 * [`Binomial`](@ref)
 * [`NegativeBinomial`](@ref)
@@ -18,6 +16,7 @@ and one of
 * [`Exponential`](@ref)
 * [`Gamma`](@ref)
 * [`MvNormal`](@ref)
+* [`Stable`](@ref)
 
 External links: [List of convolutions of probability distributions on Wikipedia](https://en.wikipedia.org/wiki/List_of_convolutions_of_probability_distributions)
 """
@@ -62,6 +61,19 @@ function convolve(d1::Gamma, d2::Gamma)
     return Gamma(d1.α + d2.α, d1.θ)
 end
 
+function convolve(d1::Stable, d2::Stable)
+    d1.α ≈ d2.α || throw(ArgumentError("$(d1.α) !≈ $(d2.α): α parameters must be approximately equal"))
+    α, β₁, σ₁, μ₁ = params(d1)
+    _, β₂, σ₂, μ₂ = params(d2)
+
+    return Stable(
+        α,                                   # α
+        (β₁*σ₁^α + β₂*σ₂^α) / (σ₁^α + σ₂^α), # β
+        (σ₁^α + σ₂^α)^(1/α),                 # σ
+        μ₁ + μ₂                              # μ
+    )
+end
+
 # continuous multivariate
 function convolve(d1::MvNormal, d2::MvNormal)
     _check_convolution_shape(d1, d2)