diff --git a/src/Distributions.jl b/src/Distributions.jl
index b18dfa8086..4c7deac44e 100644
--- a/src/Distributions.jl
+++ b/src/Distributions.jl
@@ -170,6 +170,7 @@ export
     TriangularDist,
     Triweight,
     Truncated,
+    Tweedie,
     Uniform,
     UnivariateGMM,
     VonMises,
@@ -363,7 +364,7 @@ Supported distributions:
     NoncentralF, NoncentralHypergeometric, NoncentralT, Normal, NormalCanon,
     NormalInverseGaussian, Pareto, PGeneralizedGaussian, Poisson, PoissonBinomial,
     QQPair, Rayleigh, Rician, Skellam, Soliton, StudentizedRange, SymTriangularDist, TDist, TriangularDist,
-    Triweight, Truncated, Uniform, UnivariateGMM,
+    Triweight, Truncated, Tweedie, Uniform, UnivariateGMM,
     VonMises, VonMisesFisher, WalleniusNoncentralHypergeometric, Weibull,
     Wishart, ZeroMeanIsoNormal, ZeroMeanIsoNormalCanon,
     ZeroMeanDiagNormal, ZeroMeanDiagNormalCanon, ZeroMeanFullNormal,
diff --git a/src/univariate/continuous/tweedie.jl b/src/univariate/continuous/tweedie.jl
new file mode 100644
index 0000000000..6e559f290c
--- /dev/null
+++ b/src/univariate/continuous/tweedie.jl
@@ -0,0 +1,248 @@
+"""
+    Tweedie(μ,σ,p)
+
+The *Tweedie distribution* with mean `μ ≥ 0`, dispersion `σ ≥ 0` and power `1 ≥ p ≥ 2`.
+When ``p = 1`` and ``\\sigma = 1`` it is equivalent to a quasi-Poisson distribution,
+and when ``p = 2`` to the Gamma distribution. When ``1 > p > 2``, it is a compound
+Poisson-Gamma distribution, with probability density function:
+
+```math
+f(x; \\mu, \\sigma, p) = \\frac{1}{x} W(x, \\sigma^2, p) exp \\left(
+    \\frac{1}{\\sigma^2}} [ x \\frac{\\mu^(1-p)}{1-p} - \\frac{\\mu^(2-p)}{2-p} ]
+    \\right), \\quad x > 0
+```
+where ``W`` is [Wright's generalized Bessel function](https://en.wikipedia.org/wiki/Bessel%E2%80%93Maitland_function).
+
+Note that if ``1 > p > 2`` then the distribution is continuous with a point mass concentrated at zero.
+If ``p = 1`` then the distribution is discrete.
+
+Computation of [`pdf`](@ref) and [`logpdf`](@ref) is carried out using `Float64`.
+Accuracy is generally higher than 1e-11, though for some parameter values it can
+be as low as 1e-8.
+
+```julia
+Tweedie(μ, σ, p) # Tweedie distribution with location μ, scale σ and power p
+
+params(d)        # Get the parameters, i.e. (μ, σ, p)
+location(d)      # Get the location parameter, i.e. μ
+scale(d)         # Get the scale parameter, i.e. σ
+shape(d)         # Get the shape parameter, i.e. p
+
+mean(d)          # Get the mean, i.e. μ
+var(d)           # Get the variance, i.e. σ^2 * μ^p
+```
+
+External links
+
+- [Tweedie distribution on Wikipedia](https://en.wikipedia.org/wiki/Tweedie_distribution)
+- [Compound Poisson distribution on Wikipedia](https://en.wikipedia.org/wiki/Compound_Poisson_distribution)
+
+References
+
+- Dunn P. K., Smyth G. K. (2005). "Series evaluation of Tweedie exponential dispersion model densities"
+  *Statistics and Computing* 15: 267–280.
+"""
+struct Tweedie{T <: Real} <: ContinuousUnivariateDistribution
+    μ::T
+    σ::T
+    p::T
+
+    Tweedie{T}(µ::T, σ::T, p::T) where {T<:Real} = new{T}(µ, σ, p)
+end
+
+function Tweedie(μ::T, σ::T, p::T; check_args::Bool=true) where {T <: Real}
+    @check_args(
+        Tweedie,
+        (μ, μ >= 0),
+        (σ, σ >= 0),
+        (p, 1 <= p <= 2)
+    )
+    return Tweedie{T}(μ, σ, p)
+end
+
+#### Outer constructors
+Tweedie(μ::Real, σ::Real, p::Real; check_args::Bool=true) =
+    Tweedie(promote(μ, σ, p)...; check_args=check_args)
+Tweedie(μ::Integer, σ::Integer, p::Integer; check_args::Bool=true) =
+    Tweedie(float(μ), float(σ), float(p); check_args=check_args)
+
+#### Conversions
+convert(::Type{Tweedie{T}}, μ::S, σ::S, p::S) where {T <: Real, S <: Real} = Tweedie(T(μ), T(σ), T(p))
+Base.convert(::Type{Tweedie{T}}, d::Tweedie) where {T<:Real} = Tweedie{T}(T(d.μ), T(d.σ), T(d.p))
+Base.convert(::Type{Tweedie{T}}, d::Tweedie{T}) where {T<:Real} = d
+
+@distr_support Tweedie 0 Inf
+
+#### Parameters
+
+params(d::Tweedie) = (d.μ, d.σ, d.p)
+partype(::Tweedie{T}) where {T} = T
+
+location(d::Tweedie) = d.μ
+scale(d::Tweedie) = d.σ
+shape(d::Tweedie) = d.p
+
+Base.eltype(::Type{Tweedie{T}}) where {T} = float(T)
+
+#### Statistics
+
+mean(d::Tweedie) = float(d.μ)
+var(d::Tweedie) = d.σ^2 * d.μ^d.p
+std(d::Tweedie) = d.σ * d.μ^(d.p/2)
+
+# Clark, David R. and Charles A. Thayer. 2004.
+# “A Primer on the Exponential Family of Distributions.” CAS Discussion Paper Program, 117-148
+# https://www.casact.org/sites/default/files/database/dpp_dpp04_04dpp117.pdf
+skewness(d::Tweedie) = d.p * d.σ / sqrt(d.μ ^ (2 - d.p))
+kurtosis(d::Tweedie) = d.p * (2 * d.p - 1) * d.σ^2 / d.μ ^ (2 - d.p)
+
+function logpdf(d::Tweedie{T}, x::Real)::promote_type(eltype(d), typeof(x)) where {T <: Real}
+    isnan(x) && return NaN
+    x >= 0 || return -Inf
+    # See: Dunn, Smyth (2005). "Series evaluation of Tweedie exponential dispersion model densities"
+    # Statistics and Computing 15: 267–280.
+    # pdf(y, μ, p, ϕ) = f(y, θ, ϕ) = c(y, ϕ) * exp(1/ϕ (y θ - κ(θ)))
+    # κ = cumulant function
+    # θ = function of expectation μ and power p
+    # α = (2-p)/(1-p)
+    # ϕ = σ^2
+    # y = x
+    # for 1<p<2:
+    # c(y, ϕ) = 1/y * wrightbessel(a, b, z)
+    # a = -α
+    # b = 0
+    # z = (p-1)^α/(2-p) / y^α / ϕ^(1-α)
+    μ, p, ϕ = d.μ, d.p, d.σ^2
+    if p == 1
+        return logpdf(Poisson(μ / ϕ), x / ϕ)
+    elseif p == 2
+        return logpdf(Gamma(1 / ϕ, μ * ϕ), x)
+    else
+        θ = μ ^ (1 - p) / (1 - p)
+        κ = μ ^ (2 - p) / (2 - p)
+        α = (2 - p) / (1 - p)
+
+        res = (x * θ - κ) / ϕ
+        if x > 0
+            z = ((p - 1) * ϕ / x) ^ α / ((2 - p) * ϕ)
+            # Use log to reduce risks of overflow when p is close to 1
+            wb = logwrightbessel(Float64(-α), 0.0, Float64(z))
+            # Overflow in `logwrightbessel` doesn't generally indicate that the PDF
+            # value would be larger than `typemax(Float64)`
+            wb == Inf && return NaN
+            res += wb - log(x)
+        end
+        return res
+    end
+end
+
+function cdf(d::Tweedie, x::Real)::promote_type(eltype(d), typeof(x))
+    isnan(x) && return NaN
+    x == Inf && return 1
+    x >= 0 || return 0
+    μ = d.μ
+    p = d.p
+    ϕ = d.σ^2
+    if p == 1
+        return cdf(Poisson(μ / ϕ), x / ϕ)
+    elseif p == 2
+        return cdf(Gamma(1 / ϕ, μ * ϕ), x)
+    else
+        # the mass at zero has to be handled separately as `quadgk` never evaluates at bounds
+        return pdf(d, 0) + quadgk(xi -> pdf(d, xi), 0, x, rtol=1e-12)[1]
+    end
+end
+
+function rand(rng::AbstractRNG, d::Tweedie)
+    μ, p, ϕ = d.μ, d.p, d.σ^2
+    # note that sources often use β = 1/θ for Gamma distribution
+    # e.g. https://en.wikipedia.org/wiki/Compound_Poisson_distribution
+    if p == 1
+        return ϕ * rand(rng, Poisson(μ / ϕ))
+    elseif p == 2
+        return rand(rng, Gamma(1 / ϕ, μ * ϕ))
+    else
+        λ = μ^(2 - p) / ((2 - p) * ϕ)
+        α = (2 - p) / (1 - p)
+        θ = ((p - 1) * ϕ) / μ^(1 - p)
+        N = rand(rng, Poisson(λ))
+        return N == 0 ? zero(θ) : rand(rng, Gamma(- N * α, θ))
+    end
+end
+
+# Implementation inspired by `qtweedie` in R package tweedie
+# licensed under MIT with authorization from Peter Dunn
+function quantile(d::Tweedie{T}, q::Real)::eltype(d) where {T <: Real}
+    μ, ϕ, p = d.μ, d.σ^2, d.p
+
+    if q == 0
+        return zero(T)
+    elseif q == 1
+        return convert(T, Inf)
+    elseif q < 0 || q > 1
+        throw(DomainError(q, "q must be between 0 and 1"))
+    end
+
+    if p == 1
+        return ϕ * quantile(Poisson(μ / ϕ), q)
+    elseif p == 2
+        return quantile(Gamma(1 / ϕ, μ * ϕ), q)
+    else
+        # Handle point mass at zero
+        p_zero = pdf(d, 0)
+        if q <= p_zero
+            return zero(T)
+        end
+
+        # Starting values via interpolation between Poisson and Gamma quantiles
+        qp = ϕ * quantile(Poisson(μ / ϕ), q)
+        qg = quantile(Gamma(1 / ϕ, μ * ϕ), q)
+        startx = (qg - qp) * p + (2 * qp - qg)
+
+        qstart = cdf(d, startx)
+        rx = lx = startx
+        if qstart == q
+            return startx
+        elseif qstart > q
+            while true
+                lx = lx / 2
+                cdf(d, lx) < q && break
+            end
+        elseif qstart < q
+            while true
+                rx = 1.5 * (rx + 2)
+                cdf(d, rx) > q && break
+            end
+        end
+
+        # Cannot use `quantile_newton` as pdf is sometimes multimodal
+        return quantile_bisect(d, q, lx, rx)
+    end
+end
+
+function cquantile(d::Tweedie, q::Real)
+    0 <= q <= 1 || throw(DomainError(q, "q must be between 0 and 1"))
+    cq = 1.0 - q
+    # Allow for 1 eps tolerance as due to the mass at zero
+    # if `1 - q` is rounded up when storing in floating point,
+    # `cquantile(d, ccdf(d, 0))` can be very different from zero,
+    # which doesn't make mathematical sense
+    if d.p < 2 && cq <= nextfloat(pdf(d, 0))
+        return zero(eltype(d))
+    else
+        return quantile(d, cq)
+    end
+end
+
+function invlogccdf(d::Tweedie, lp::Real)
+    p = -expm1(lp)
+    # Allow for 1 eps tolerance as due to the mass at zero
+    # if `1 - q` is rounded up when storing in floating point,
+    # `invlogccdf(d, logccdf(d, 0))` can be very different from zero,
+    # which doesn't make mathematical sense
+    if d.p < 2 && p <= nextfloat(pdf(d, 0))
+        return zero(eltype(d))
+    else
+        return quantile(d, p)
+    end
+end
\ No newline at end of file
diff --git a/src/univariates.jl b/src/univariates.jl
index 36fc5a214d..86337c9d49 100644
--- a/src/univariates.jl
+++ b/src/univariates.jl
@@ -733,6 +733,7 @@ const continuous_distributions = [
     "tdist",
     "triangular",
     "triweight",
+    "tweedie",
     "uniform",
     "loguniform", # depends on Uniform
     "vonmises",
diff --git a/test/ref/continuous/tweedie.R b/test/ref/continuous/tweedie.R
new file mode 100644
index 0000000000..20a8bb29a1
--- /dev/null
+++ b/test/ref/continuous/tweedie.R
@@ -0,0 +1,32 @@
+# Tweedie Distribution
+# Using R's tweedie package for reference implementation
+
+Tweedie <- R6Class("Tweedie",
+    inherit = ContinuousDistribution,
+    public = list(
+        names = c("mu", "sigma", "p"),
+        mu = NA,
+        sigma = NA,
+        p = NA,
+        initialize = function(mu, sigma, p) {
+            self$mu <- mu
+            self$sigma <- sigma
+            self$p <- p
+        },
+        supp = function() { c(0, Inf) },
+        pdf = function(x, log=FALSE) {
+            val <- tweedie::dtweedie(x, mu=self$mu, phi=self$sigma^2, power=self$p)
+            if (log) {
+                return(log(val))
+            } else {
+                return(val)
+            }
+        },
+        cdf = function(x) {
+            tweedie::ptweedie(x, mu=self$mu, phi=self$sigma^2, power=self$p)
+        },
+        quan = function(v) {
+            tweedie::qtweedie(v, mu=self$mu, phi=self$sigma^2, power=self$p)
+        }
+    )
+)
diff --git a/test/ref/continuous_test.lst b/test/ref/continuous_test.lst
index a0a5418314..d0d206e292 100644
--- a/test/ref/continuous_test.lst
+++ b/test/ref/continuous_test.lst
@@ -209,6 +209,12 @@ TruncatedNormal(27, 3, 0, Inf)
 TruncatedNormal(-5, 1, -Inf, -10)
 TruncatedNormal(1.8, 1.2, -Inf, 0)
 
+Tweedie(2.0, 1.5, 1)
+Tweedie(2.0, 1.5, 1.01)
+Tweedie(2.0, 1.5, 1.5)
+Tweedie(2.0, 1.5, 1.99)
+Tweedie(2.0, 1.5, 2.0)
+
 Uniform()
 Uniform(0.0, 2.0)
 Uniform(3.0, 17.0)
diff --git a/test/ref/continuous_test.ref.json b/test/ref/continuous_test.ref.json
index 7f8fb4b182..f262c437a3 100644
--- a/test/ref/continuous_test.ref.json
+++ b/test/ref/continuous_test.ref.json
@@ -5457,6 +5457,136 @@
     { "q": 0.90, "x": -0.0644552402816942 }
   ]
 },
+{
+  "expr": "Tweedie(2.0, 1.5, 1)",
+  "dtype": "Tweedie",
+  "minimum": 0,
+  "maximum": "inf",
+  "properties": {
+  },
+  "points": [
+    { "x": 0, "pdf": 0.411112290507187, "logpdf": -0.888888888888889, "cdf": 0.411112290507187 },
+    { "x": 0, "pdf": 0.411112290507187, "logpdf": -0.888888888888889, "cdf": 0.411112290507187 },
+    { "x": 0, "pdf": 0.411112290507187, "logpdf": -0.888888888888889, "cdf": 0.411112290507187 },
+    { "x": 0, "pdf": 0.411112290507187, "logpdf": -0.888888888888889, "cdf": 0.411112290507187 },
+    { "x": 2.25, "pdf": 0.3654331471175, "logpdf": -1.00667192454527, "cdf": 0.776545437624687 },
+    { "x": 2.25, "pdf": 0.3654331471175, "logpdf": -1.00667192454527, "cdf": 0.776545437624687 },
+    { "x": 2.25, "pdf": 0.3654331471175, "logpdf": -1.00667192454527, "cdf": 0.776545437624687 },
+    { "x": 4.5, "pdf": 0.162414732052222, "logpdf": -1.8176021407616, "cdf": 0.93896016967691 },
+    { "x": 4.5, "pdf": 0.162414732052222, "logpdf": -1.8176021407616, "cdf": 0.93896016967691 }
+  ],
+  "quans": [
+    { "q": 0.10, "x": -0 },
+    { "q": 0.25, "x": 0 },
+    { "q": 0.50, "x": 2.25 },
+    { "q": 0.75, "x": 2.25 },
+    { "q": 0.90, "x": 4.5 }
+  ]
+},
+{
+  "expr": "Tweedie(2.0, 1.5, 1.01)",
+  "dtype": "Tweedie",
+  "minimum": 0,
+  "maximum": "inf",
+  "properties": {
+  },
+  "points": [
+    { "x": 0, "pdf": 0.409972358341713, "logpdf": -0.891665540235274, "cdf": 0.409972358341713 },
+    { "x": 0, "pdf": 0.409972358341713, "logpdf": -0.891665540235274, "cdf": 0.409972358341713 },
+    { "x": 0, "pdf": 0.409972358341713, "logpdf": -0.891665540235274, "cdf": 0.409972358341713 },
+    { "x": 0, "pdf": 0.409972358341713, "logpdf": -0.891665540235274, "cdf": 0.409972358341713 },
+    { "x": 2.08458174204291, "pdf": 0.536709862522414, "logpdf": -0.622297623738618, "cdf": 0.499999999999887 },
+    { "x": 2.24665165372521, "pdf": 0.64524670586641, "logpdf": -0.438122545601544, "cdf": 0.599999999999865 },
+    { "x": 2.42462213545382, "pdf": 0.439428418580859, "logpdf": -0.822280445244871, "cdf": 0.7 },
+    { "x": 4.15658364531784, "pdf": 0.12549162469383, "logpdf": -2.07551625814455, "cdf": 0.8 },
+    { "x": 4.71112894150918, "pdf": 0.152496114134358, "logpdf": -1.8806161643455, "cdf": 0.9 }
+  ],
+  "quans": [
+    { "q": 0.10, "x": 0 },
+    { "q": 0.25, "x": 0 },
+    { "q": 0.50, "x": 2.08458174204291 },
+    { "q": 0.75, "x": 2.58441784358517 },
+    { "q": 0.90, "x": 4.71112894150918 }
+  ]
+},
+{
+  "expr": "Tweedie(2.0, 1.5, 1.5)",
+  "dtype": "Tweedie",
+  "minimum": 0,
+  "maximum": "inf",
+  "properties": {
+  },
+  "points": [
+    { "x": 0, "pdf": 0.284483870247908, "logpdf": -1.25707872210942, "cdf": 0.284483870247908 },
+    { "x": 0, "pdf": 0.284483870247908, "logpdf": -1.25707872210942, "cdf": 0.284483870247908 },
+    { "x": 0.0695940988351711, "pdf": 0.221126951073659, "logpdf": -1.50901830314413, "cdf": 0.3 },
+    { "x": 0.548859776983353, "pdf": 0.196301027378934, "logpdf": -1.62810594399496, "cdf": 0.400000000000001 },
+    { "x": 1.09635850238778, "pdf": 0.169308648560591, "logpdf": -1.77603190691513, "cdf": 0.5 },
+    { "x": 1.74368937890802, "pdf": 0.140219032553443, "logpdf": -1.96454956071761, "cdf": 0.6 },
+    { "x": 2.54916644046212, "pdf": 0.10903389348092, "logpdf": -2.21609649572149, "cdf": 0.7 },
+    { "x": 3.6418146049751, "pdf": 0.0756550324262397, "logpdf": -2.581571318472, "cdf": 0.8 },
+    { "x": 5.42635624081533, "pdf": 0.0397731565622146, "logpdf": -3.22456305247684, "cdf": 0.9 }
+  ],
+  "quans": [
+    { "q": 0.10, "x": 0 },
+    { "q": 0.25, "x": 0 },
+    { "q": 0.50, "x": 1.09635850238778 },
+    { "q": 0.75, "x": 3.04590854984609 },
+    { "q": 0.90, "x": 5.42635624081533 }
+  ]
+},
+{
+  "expr": "Tweedie(2.0, 1.5, 1.99)",
+  "dtype": "Tweedie",
+  "minimum": 0,
+  "maximum": "inf",
+  "properties": {
+  },
+  "points": [
+    { "x": 0.0180116841854315, "pdf": 2.36802160266867, "logpdf": 0.862054839728742, "cdf": 0.100000000000001 },
+    { "x": 0.0909089850716691, "pdf": 0.942832617388424, "logpdf": -0.0588665122220102, "cdf": 0.2 },
+    { "x": 0.235225094656205, "pdf": 0.539396455646777, "logpdf": -0.617304439216196, "cdf": 0.3 },
+    { "x": 0.468504550181304, "pdf": 0.350419323103255, "logpdf": -1.04862477559775, "cdf": 0.399999999999998 },
+    { "x": 0.81685153909239, "pdf": 0.23900152765095, "logpdf": -1.43128533522608, "cdf": 0.5 },
+    { "x": 1.32470072354968, "pdf": 0.163772302491998, "logpdf": -1.80927821531652, "cdf": 0.6 },
+    { "x": 2.07733165294402, "pdf": 0.108243564930168, "logpdf": -2.22337136023913, "cdf": 0.7 },
+    { "x": 3.2704062085269, "pdf": 0.064685097882654, "logpdf": -2.73822443039313, "cdf": 0.8 },
+    { "x": 5.54191528143532, "pdf": 0.0291464642496897, "logpdf": -3.53542166857885, "cdf": 0.9 }
+  ],
+  "quans": [
+    { "q": 0.10, "x": 0.0180116841854315 },
+    { "q": 0.25, "x": 0.153164414123896 },
+    { "q": 0.50, "x": 0.81685153909239 },
+    { "q": 0.75, "x": 2.59770635997387 },
+    { "q": 0.90, "x": 5.54191528143532 }
+  ]
+},
+{
+  "expr": "Tweedie(2.0, 1.5, 2.0)",
+  "dtype": "Tweedie",
+  "minimum": 0,
+  "maximum": "inf",
+  "properties": {
+  },
+  "points": [
+    { "x": 0.0193166098973791, "pdf": 2.29401151609469, "logpdf": 0.830302038788168, "cdf": 0.1 },
+    { "x": 0.0929286024211608, "pdf": 0.942933429728374, "logpdf": -0.0587595929780245, "cdf": 0.2 },
+    { "x": 0.236509737467784, "pdf": 0.543544813747258, "logpdf": -0.609643121823222, "cdf": 0.3 },
+    { "x": 0.467817526522702, "pdf": 0.353460636956834, "logpdf": -1.03998315216011, "cdf": 0.4 },
+    { "x": 0.813387643444831, "pdf": 0.240732079165742, "logpdf": -1.42407066859999, "cdf": 0.5 },
+    { "x": 1.31820953609136, "pdf": 0.164558740147814, "logpdf": -1.80448768955023, "cdf": 0.6 },
+    { "x": 2.06837862454134, "pdf": 0.108450167200076, "logpdf": -2.22146449998544, "cdf": 0.7 },
+    { "x": 3.26120962635577, "pdf": 0.0646025808808938, "logpdf": -2.73950091727501, "cdf": 0.8 },
+    { "x": 5.54018226362739, "pdf": 0.0290033364805552, "logpdf": -3.5403444045599, "cdf": 0.9 }
+  ],
+  "quans": [
+    { "q": 0.10, "x": 0.0193166098973791 },
+    { "q": 0.25, "x": 0.15499278218233 },
+    { "q": 0.50, "x": 0.813387643444831 },
+    { "q": 0.75, "x": 2.58817864849412 },
+    { "q": 0.90, "x": 5.54018226362739 }
+  ]
+},
 {
   "expr": "Uniform()",
   "dtype": "Uniform",
diff --git a/test/ref/rdistributions.R b/test/ref/rdistributions.R
index 8df19b70a7..4a8d1dacd0 100644
--- a/test/ref/rdistributions.R
+++ b/test/ref/rdistributions.R
@@ -77,6 +77,7 @@ source("continuous/studentizedrange.R")
 source("continuous/tdist.R")
 source("continuous/triangulardist.R")
 source("continuous/truncatednormal.R")
+source("continuous/tweedie.R")
 source("continuous/uniform.R")
 source("continuous/vonmises.R")
 source("continuous/weibull.R")
diff --git a/test/ref/renv.lock b/test/ref/renv.lock
index 1e520649f8..8aadd1755a 100644
--- a/test/ref/renv.lock
+++ b/test/ref/renv.lock
@@ -1,6 +1,6 @@
 {
   "R": {
-    "Version": "4.5.1",
+    "Version": "4.5.2",
     "Repositories": [
       {
         "Name": "CRAN",
@@ -1335,6 +1335,12 @@
       "Maintainer": "Georgi N. Boshnakov <georgi.boshnakov@manchester.ac.uk>",
       "Repository": "CRAN"
     },
+    "tweedie": {
+      "Package": "tweedie",
+      "Version": "3.0.17",
+      "Source": "Repository",
+      "Repository": "CRAN"
+    },
     "vctrs": {
       "Package": "vctrs",
       "Version": "0.6.5",
diff --git a/test/runtests.jl b/test/runtests.jl
index 69d3327479..28e25bcaf4 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -27,6 +27,7 @@ const tests = [
     "univariate/continuous/cauchy",
     "univariate/continuous/uniform",
     "univariate/continuous/lognormal",
+    "univariate/continuous/tweedie",
     "multivariate/mvnormal",
     "multivariate/mvlogitnormal",
     "multivariate/mvlognormal",
diff --git a/test/testutils.jl b/test/testutils.jl
index 8a62c52d3b..b7d785aaf4 100644
--- a/test/testutils.jl
+++ b/test/testutils.jl
@@ -540,7 +540,7 @@ function test_evaluation(d::ContinuousUnivariateDistribution, vs::AbstractVector
         end
     end
 
-    if !isa(d, StudentizedRange)
+    if !isa(d, StudentizedRange) && !(isa(d, Tweedie) && d.p == 1)
         # check: pdf should be the derivative of cdf
         for i = 2:(nv-1)
             if p[i] > 1.0e-6
diff --git a/test/univariate/continuous/tweedie.jl b/test/univariate/continuous/tweedie.jl
new file mode 100644
index 0000000000..414428f0c1
--- /dev/null
+++ b/test/univariate/continuous/tweedie.jl
@@ -0,0 +1,89 @@
+@testset "Tweedie quantile bounds" begin
+    d = Tweedie(2.0, 1.5, 1.0)
+    @test quantile(d, 0.0) == 0.0
+    @test quantile(d, 1.0) == Inf
+end
+
+@testset "Tweedie pdf" begin
+    d = Tweedie(2.0, 1.5, 1.0)
+    @test logpdf(d, -0.3) == -Inf
+    @test pdf(d, -0.3) == 0
+    @test logpdf(d, -0.0) == logpdf(d, 0.0)
+    @test pdf(d, -0.0) == pdf(d, 0.0)
+end
+
+@testset "Tweedie overflows" begin
+    d = Tweedie(1.1, 1.0, 1.01)
+    # This case overflows when using `wrightbessel` rather than `logwrightbessel`
+    @test pdf(d, 10) ≈ 3.1785795511027634e-7
+    # This one would return Inf if we don't manually return NaN
+    @test isnan(pdf(d, 15))
+end
+
+
+@testset "Tweedie elementary statistics" begin
+    d = Tweedie(2.0, 1.5, 1.0)
+    @test mean(d) == 2.0
+    @test var(d) == 4.5
+    @test std(d) ≈ sqrt(4.5)
+    @test skewness(d) ≈ 1.0606601717798212
+    @test kurtosis(d) == 1.125
+end
+
+@testset "Tweedie custom type" begin
+    d1 = Tweedie(2.0, 1.5, 1.0)
+    d2 = Tweedie(2.0f0, 1.5f0, 1.0f0)::Tweedie{Float32}
+    @test @inferred(pdf(d2, 1))::Float32 ≈ pdf(d1, 1)
+    @test @inferred(cdf(d2, 1))::Float32 ≈ cdf(d1, 1)
+    @test @inferred(quantile(d2, 0.1))::Float32 ≈ quantile(d1, 0.1)
+    @test @inferred(mean(d2))::Float32 ≈ mean(d1)
+    @test @inferred(median(d2))::Float32 ≈ median(d1)
+    @test @inferred(skewness(d2))::Float32 ≈ skewness(d1)
+    @test @inferred(kurtosis(d2))::Float32 ≈ kurtosis(d1)
+    # return type depends on distribution and argument type via promotion
+    @test @inferred(pdf(d2, 1.0))::Float64 ≈ pdf(d1, 1.0)
+    @test @inferred(cdf(d2, 1.0))::Float64 ≈ cdf(d1, 1.0)
+
+
+    d1 = Tweedie(2.0, 2.0, 2.0)
+    d2 = Tweedie{Int}(2, 2, 2)::Tweedie{Int}
+    @test @inferred(pdf(d2, 1))::Float64 ≈ pdf(d1, 1)
+    @test @inferred(cdf(d2, 1))::Float64 ≈ cdf(d1, 1)
+    @test @inferred(quantile(d2, 0.1))::Float64 ≈ quantile(d1, 0.1)
+    @test @inferred(mean(d2))::Float64 ≈ mean(d1)
+    @test @inferred(median(d2))::Float64 ≈ median(d1)
+    @test @inferred(skewness(d2))::Float64 ≈ skewness(d1)
+    @test @inferred(kurtosis(d2))::Float64 ≈ kurtosis(d1)
+    # return type depends on distribution and argument type via promotion
+    @test @inferred(pdf(d2, big(1)))::BigFloat ≈ pdf(d1, 1.0)
+    @test @inferred(cdf(d2, big(1)))::BigFloat ≈ cdf(d1, 1.0)
+
+    for d2 in (
+        Tweedie(2.0, 1.5f0, 1.0f0),
+        Tweedie(2.0f0, 1.5, 1.0f0),
+        Tweedie(2.0f0, 1.5f0, 1.0),
+        Tweedie(2.0, 1.5, 1.0f0),
+        Tweedie(2.0f0, 1.5, 1.0),
+
+        Tweedie(2, 2.0, 1.0),
+        Tweedie(2.0, 2, 1.0),
+        Tweedie(2.0, 2.0, 1),
+        Tweedie(2, 2, 1.0),
+        Tweedie(2.0, 2, 1),
+        )
+        @test d2 isa Tweedie{Float64}
+    end
+end
+
+@testset "Tweedie (c)quantile and invlog(c)cdf" begin
+    # check that rounding from exact computation to floating point
+    # does not break round-tripping due to mass at zero
+    d = Tweedie(1.1, 1, 1.01)
+    @test quantile(d, cdf(d, 0)) == 0
+    @test cquantile(d, ccdf(d, 0)) == 0
+    @test invlogcdf(d, logcdf(d, 0)) == 0
+    @test invlogccdf(d, logccdf(d, 0)) == 0
+
+    @test_throws DomainError cquantile(d, -0.1)
+    @test_throws DomainError cquantile(d, 1.1)
+end
diff --git a/test/univariates.jl b/test/univariates.jl
index a806deb9c0..e3c5726037 100644
--- a/test/univariates.jl
+++ b/test/univariates.jl
@@ -99,12 +99,18 @@ function verify_and_test(D::Type, d::UnivariateDistribution, dct::AbstractDict,
         # pdf method is not implemented for StudentizedRange
         if !isa(d, StudentizedRange)
             @test Base.Fix1(pdf, d).(x) ≈ p atol=1e-16 rtol=1e-8
-            @test Base.Fix1(logpdf, d).(x) ≈ lp atol=isa(d, NoncentralHypergeometric) ? 1e-4 : 1e-12
+            atol = isa(d, NoncentralHypergeometric) ? 1e-4 :
+                   isa(d, Tweedie) ? 1e-9 :
+                   1e-12
+            @test Base.Fix1(logpdf, d).(x) ≈ lp atol=atol
         end
 
         # cdf method is not implemented for NormalInverseGaussian
         if !isa(d, NormalInverseGaussian)
-            @test isapprox(cdf(d, x), cf; atol=isa(d, NoncentralHypergeometric) ? 1e-8 : 1e-12)
+            atol = isa(d, NoncentralHypergeometric) ? 1e-8 :
+                   isa(d, Tweedie) ? 1e-11 :
+                   1e-12
+            @test isapprox(cdf(d, x), cf; atol=atol)
         end
     end