Add Wrapped distribution wrapper #1724
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,317 @@ | ||
| using Distributions | ||
|
|
||
| """ | ||
| wrapped(d0::UnivariateDistribution, l::Real, u::Real; k::Int=1, kwargs...) | ||
|
|
||
| Construct a ``k``-wrapped distribution with support ``[l, u)``. | ||
|
|
||
| The resulting density wraps around the interval. With the period ``\\tau = u - l``, the | ||
| density can be expressed as two infinite series in terms of either the density ``f_0`` of | ||
| the unwrapped distribution `d0` or in terms of its characteristic function [`cf`](@ref) | ||
| ``\\phi_0``: | ||
| ````math | ||
| \\begin{aligned} | ||
| f(x) &= \\frac{1}{\\tau}\\sum_{j=-\\infty}^\\infty f_0\\left(x + \\frac{j}{k}\\tau\\right)\\ | ||
| &= \\frac{1}{\\tau}\\sum_{j=-\\infty}^\\infty | ||
| \\phi_0\\left(jk\\frac{2\\pi}{\\tau}\\right) | ||
| \\exp\\left(-ijk\\frac{2\\pi}{\\tau}x\\right) | ||
| \\end{aligned} | ||
| ```` | ||
| If `k>1`, then the resulting distribution has ``k``-fold symmetry in the interval. | ||
|
|
||
| # Keywords | ||
| - `characteristic::Bool`: whether to use the characteristic function to compute the density. | ||
| If `false`, the density is computed using the logpdf of the unwrapped distribution. | ||
| If not specified, the `wrapped` method called may choose based on the parameters of | ||
| `l`, `u`, `tol`, and `d0`. | ||
| - `tol`: target absolute error of `pdf`. May be used to determine whether to use the | ||
| characteristic function in computations or not. | ||
|
|
||
| # Examples | ||
|
|
||
| """ | ||
| function wrapped(d::UnivariateDistribution, l::Real, u::Real; kwargs...) | ||
| return Wrapped(d, promote(l, u)...; kwargs...) | ||
| end | ||
| function wrapped(d::UnivariateDistribution, l::T, u::T; kwargs...) where {T<:Real} | ||
| return Wrapped(d, l, u; kwargs...) | ||
| end | ||
|
|
||
| """ | ||
| Wrapped | ||
|
|
||
| Generic wrapper for a wrapped distribution. | ||
| """ | ||
| struct Wrapped{D<:UnivariateDistribution,S<:ValueSupport,T<:Real} <: | ||
| UnivariateDistribution{S} | ||
| unwrapped::D # the original distribution (untruncated) | ||
| lower::T # lower bound | ||
| upper::T # upper bound | ||
| characteristic::Bool # whether to use characteristic function | ||
| k::Int # k-fold symmetry | ||
| function Wrapped( | ||
| d::UnivariateDistribution, | ||
| l::T, | ||
| u::T; | ||
| tol::Real=sqrt(eps(float(T))), | ||
| characteristic::Bool=_wrapped_use_characteristic(d, l, u, tol), | ||
| k::Int=1, | ||
| ) where {T<:Real} | ||
| return new{typeof(d),value_support(typeof(d)),T}( | ||
| d, l, u, characteristic, k | ||
| ) | ||
| end | ||
| end | ||
|
|
||
| _wrapped_use_characteristic(d::UnivariateDistribution, l, u, tol) = false | ||
|
|
||
| minimum(d::Wrapped) = d.lower | ||
| maximum(d::Wrapped) = d.upper - eps(float(d.upper)) | ||
|
There was a problem hiding this comment. For There was a problem hiding this comment. For a wrapped distribution, See for example the wrapped Poisson example in #1716 (comment). Because 0 and 15 are the same point, the PMF sums to greater than 1 because A compromise could be returning There was a problem hiding this comment.
Distributions defines the support of a There was a problem hiding this comment.
Oh really? But this doesn't seem to be consistently enforced, e.g. julia> d = Binomial(10, 0.5) * 1//2;
julia> support(d)
0//1:1//2:5//1
julia> d = DiscreteNonParametric(randn(5), normalize(rand(5), 1));
julia> support(d)
5-element Vector{Float64}:
-0.9108476704164078
-0.5194957770100691
-0.024410217492521903
0.15370458561285605
1.338024968131234There was a problem hiding this comment. Hah, amazing. 🤦 I was going based on the generic definition for |
||
|
|
||
| function logpdf( | ||
| d::Wrapped, x::Real; characteristic::Bool=d.characteristic, kwargs... | ||
| ) | ||
| if characteristic | ||
| log(pdf_wrapped_cf(d, x; kwargs...)) | ||
| else | ||
| logpdf_wrapped_logpdf(d, x; kwargs...) | ||
| end | ||
| end | ||
|
|
||
| logpdf_wrapped_cf(d::Wrapped, x::Real; kwargs...) = log(pdf_wrapped_cf(d, x; kwargs...)) | ||
|
|
||
| function logpdf_wrapped_logpdf( | ||
| d::Wrapped, x::Real; tol=sqrt(eps(float(typeof(x)))), maxiter=1_000 | ||
| ) | ||
| d_unwrapped = d.unwrapped | ||
| period = d.upper - d.lower | ||
| k = d.k | ||
| logp = logpdf(d_unwrapped, x) | ||
| insupport(d, x) || return oftype(logp, -Inf) | ||
| Δx = period / k | ||
| x₋ = x - Δx | ||
| x₊ = x + Δx | ||
| x₋_in_support = insupport(d_unwrapped, x₋) | ||
| x₊_in_support = insupport(d_unwrapped, x₊) | ||
| for i in 1:maxiter | ||
| Δlogp = oftype(logp, -Inf) | ||
| if x₋_in_support | ||
| Δlogp = logaddexp(Δlogp, logpdf(d_unwrapped, x₋)) | ||
| x₋ -= Δx | ||
| x₋_in_support = insupport(d_unwrapped, x₋) | ||
| end | ||
| if x₊_in_support | ||
| Δlogp = logaddexp(Δlogp, logpdf(d_unwrapped, x₊)) | ||
| x₊ += Δx | ||
| x₊_in_support = insupport(d_unwrapped, x₊) | ||
| end | ||
| logp = logaddexp(logp, Δlogp) | ||
| i > k && abs(Δlogp) ≤ tol && isfinite(logp) && break | ||
| end | ||
| return logp - log(oftype(logp, k)) | ||
| end | ||
|
|
||
| function pdf( | ||
| d::Wrapped, x::Real; characteristic::Bool=d.characteristic, kwargs... | ||
| ) | ||
| if characteristic | ||
| pdf_wrapped_cf(d, x; kwargs...) | ||
| else | ||
| pdf_wrapped_pdf(d, x; kwargs...) | ||
| end | ||
| end | ||
|
|
||
| function pdf_wrapped_pdf( | ||
| d::Wrapped, x::Real; tol=sqrt(eps(float(typeof(x)))), maxiter=1_000 | ||
| ) | ||
| d_unwrapped = d.unwrapped | ||
| period = d.upper - d.lower | ||
| k = d.k | ||
| p = pdf(d_unwrapped, x) | ||
| insupport(d, x) || return zero(p) | ||
| Δx = period / k | ||
| x₋ = x - Δx | ||
| x₊ = x + Δx | ||
| x₋_in_support = insupport(d_unwrapped, x₋) | ||
| x₊_in_support = insupport(d_unwrapped, x₊) | ||
| for i in 1:maxiter | ||
| Δp = zero(p) | ||
| if x₋_in_support | ||
| Δp += pdf(d_unwrapped, x₋) | ||
| x₋ -= Δx | ||
| x₋_in_support = insupport(d_unwrapped, x₋) | ||
| end | ||
| if x₊_in_support | ||
| Δp += pdf(d_unwrapped, x₊) | ||
| x₊ += Δx | ||
| x₊_in_support = insupport(d_unwrapped, x₊) | ||
| end | ||
| p += Δp | ||
| i > k && abs(Δp) ≤ tol && p > 0 && break | ||
| end | ||
| return p / k | ||
| end | ||
|
|
||
| function cdf( | ||
| d::Wrapped, x::Real; characteristic::Bool=d.characteristic, kwargs... | ||
| ) | ||
| if characteristic && !(d.unwrapped isa DiscreteDistribution) | ||
| cdf_wrapped_cf(d, x; kwargs...) | ||
| else | ||
| cdf_wrapped_cdf(d, x; kwargs...) | ||
| end | ||
| end | ||
|
|
||
| function logcdf( | ||
| d::Wrapped, x::Real; characteristic::Bool=d.characteristic, kwargs... | ||
| ) | ||
| if characteristic && !(d.unwrapped isa DiscreteDistribution) | ||
| log(cdf_wrapped_cf(d, x; kwargs...)) | ||
| else | ||
| logcdf_wrapped_logcdf(d, x; kwargs...) | ||
| end | ||
| end | ||
|
|
||
| function cdf_wrapped_cdf( | ||
| d::Wrapped, x::Real; tol=sqrt(eps(float(typeof(x)))), maxiter=1_000 | ||
| ) | ||
| d_unwrapped = d.unwrapped | ||
| l = d.lower | ||
| period = d.upper - l | ||
| k = d.k | ||
| p = cdf(d_unwrapped, x) - cdf(d_unwrapped, l) | ||
| x < l && return zero(p) | ||
| Δx = period / k | ||
| x₋ = x - Δx | ||
| x₊ = x + Δx | ||
| x₋_in_support = insupport(d_unwrapped, x₋) | ||
| l₋ = l - Δx | ||
| l₊ = l + Δx | ||
| l₊_in_support = insupport(d_unwrapped, l₊) | ||
| for i in 1:maxiter | ||
| Δp = zero(p) | ||
| if x₋_in_support | ||
| Δp += cdf(d_unwrapped, x₋) - cdf(d_unwrapped, l₋) | ||
| if d_unwrapped isa DiscreteDistribution | ||
| Δp += pdf(d_unwrapped, l₋) | ||
| end | ||
| x₋ -= Δx | ||
| l₋ -= Δx | ||
| x₋_in_support = insupport(d_unwrapped, x₋) | ||
| end | ||
| if l₊_in_support | ||
| Δp += cdf(d_unwrapped, x₊) - cdf(d_unwrapped, l₊) | ||
| if d_unwrapped isa DiscreteDistribution | ||
| Δp += pdf(d_unwrapped, l₊) | ||
| end | ||
| x₊ += Δx | ||
| l₊ += Δx | ||
| l₊_in_support = insupport(d_unwrapped, l₊) | ||
| end | ||
| p += Δp | ||
| i > k && abs(Δp) ≤ tol && p > 0 && break | ||
| end | ||
| return p / k | ||
| end | ||
|
|
||
| function logcdf_wrapped_logcdf( | ||
| d::Wrapped, x::Real; tol=sqrt(eps(float(typeof(x)))), maxiter=1_000 | ||
| ) | ||
| d_unwrapped = d.unwrapped | ||
| l = d.lower | ||
| period = d.upper - l | ||
| k = d.k | ||
| x < l && return oftype(logcdf(d_unwrapped, l), -Inf) | ||
| logp = logdiffcdf(d_unwrapped, x, l) | ||
| x < l && return oftype(logp, -Inf) | ||
| Δx = period / k | ||
| x₋ = x - Δx | ||
| x₊ = x + Δx | ||
| x₋_in_support = insupport(d_unwrapped, x₋) | ||
| l₋ = l - Δx | ||
| l₊ = l + Δx | ||
| l₊_in_support = insupport(d_unwrapped, l₊) | ||
| for i in 1:maxiter | ||
| Δlogp = oftype(logp, -Inf) | ||
| if x₋_in_support | ||
| Δlogp = logaddexp(Δlogp, logdiffcdf(d_unwrapped, x₋, l₋)) | ||
| if d_unwrapped isa DiscreteDistribution | ||
| Δlogp = logaddexp(Δlogp, logpdf(d_unwrapped, l₋)) | ||
| end | ||
| x₋ -= Δx | ||
| l₋ -= Δx | ||
| x₋_in_support = insupport(d_unwrapped, x₋) | ||
| end | ||
| if l₊_in_support | ||
| Δlogp = logaddexp(Δlogp, logdiffcdf(d_unwrapped, x₊, l₊)) | ||
| if d_unwrapped isa DiscreteDistribution | ||
| Δlogp = logaddexp(Δlogp, logpdf(d_unwrapped, l₊)) | ||
| end | ||
| x₊ += Δx | ||
| l₊ += Δx | ||
| l₊_in_support = insupport(d_unwrapped, l₊) | ||
| end | ||
| logp = logaddexp(logp, Δlogp) | ||
| i > k && abs(Δlogp) ≤ tol && isfinite(logp) && break | ||
| end | ||
| return min(logp - log(oftype(logp, k)), 0) | ||
| end | ||
|
|
||
| function pdf_wrapped_cf(d::Wrapped, x::Real; tol=sqrt(eps(float(typeof(x)))), maxiter=1_000) | ||
| d_unwrapped = d.unwrapped | ||
| period = d.upper - d.lower | ||
| k = d.k | ||
| scale = twoπ * (k / period) | ||
| iscale = zero(scale) | ||
| z₀ = cis(-x * scale) | ||
| z = one(z₀) | ||
| p = real(cf(d_unwrapped, iscale) * z) | ||
| insupport(d, x) || return zero(p) | ||
| for i in 1:maxiter | ||
| z *= z₀ | ||
| iscale += scale | ||
| Δp = 2 * real(cf(d_unwrapped, iscale) * z) | ||
| p += Δp | ||
| i > k && abs(Δp) ≤ tol && p > 0 && break | ||
| end | ||
| # ensure density is never negative | ||
| return max(0, p) / period | ||
| end | ||
|
|
||
| function cdf_wrapped_cf(d::Wrapped, x::Real; tol=sqrt(eps(float(typeof(x)))), maxiter=1_000) | ||
| d_unwrapped = d.unwrapped | ||
| period = d.upper - d.lower | ||
| k = d.k | ||
| scale = twoπ * (k / period) | ||
| iscale = zero(scale) | ||
| z₀ = cis(-x * scale) | ||
| z₀ₗ = cis(-d.lower * scale) | ||
| z = zₗ = one(z₀) | ||
| p = real(cf(d_unwrapped, iscale)) * (x - d.lower) | ||
| for i in 1:maxiter | ||
| z *= z₀ | ||
| zₗ *= z₀ₗ | ||
| iscale += scale | ||
| Δp = 2 * imag(cf(d_unwrapped, iscale) * (zₗ - z)) / iscale | ||
| p += Δp | ||
| i > k && abs(Δp) ≤ tol && p > 0 && break | ||
| end | ||
| # ensure probability is never negative | ||
| return clamp(p / period, 0, 1) | ||
| end | ||
|
|
||
| function rand(rng::AbstractRNG, d::Wrapped) | ||
| l = d.lower | ||
| k = d.k | ||
| period = d.upper - l | ||
| # choose which of the `k` folds the point will be from | ||
| i = k == 1 ? k : rand(rng, 1:k) | ||
| z = rand(rng, d.unwrapped) | ||
| return mod(z - l + ((k - i)//k) * period, period) + l | ||
| end | ||
|
|
||
| ### specialized wrapped distributions | ||
|
|
||
| include(joinpath("wrapped", "normal.jl")) | ||
| include(joinpath("wrapped", "cauchy.jl")) | ||
| include(joinpath("wrapped", "exponential.jl")) | ||
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,47 @@ | ||
| function pdf(d::Wrapped{<:Cauchy}, x::Real; kwargs...) | ||
| d0 = d.unwrapped | ||
| l = d.lower | ||
| u = d.upper | ||
| period = u - l | ||
| scale = d.k * 2π / period | ||
| s = d0.σ * scale | ||
| z = (x - d0.μ) * scale | ||
| p = sinh(s) / (cosh(s) - cos(z)) / period | ||
| l ≤ x < u || return zero(p) | ||
| return p | ||
| end | ||
|
|
||
| function logpdf(d::Wrapped{<:Cauchy}, x::Real; kwargs...) | ||
| d0 = d.unwrapped | ||
| l = d.lower | ||
| u = d.upper | ||
| period = u - l | ||
| scale = d.k * 2π / period | ||
| s = d0.σ * scale | ||
| z = (x - d0.μ) * scale | ||
| logp = log1mexp(-2s) - logaddexp(2log1mexp(-s), logtwo - s + log1p(-cos(z))) | ||
| l ≤ x < u || return oftype(logp, -Inf) | ||
| return logp - log(oftype(logp, period)) | ||
| end | ||
|
|
||
| function cdf(d::Wrapped{<:Cauchy}, x::Real; kwargs...) | ||
| d0 = d.unwrapped | ||
| l = d.lower | ||
| u = d.upper | ||
| k = d.k | ||
| period = u - l | ||
| scale = k * π / period | ||
| s = d0.σ * scale | ||
| z = (x - d0.μ) * scale | ||
| zₗ = (l - d0.μ) * scale | ||
| tanhs = tanh(s) | ||
| p = (atan(tan(z), tanhs) - atan(tan(zₗ), tanhs)) / π / k | ||
| x < l && return zero(p) | ||
| x ≥ u && return one(p) | ||
| # handle discontinuities in the antiderivative at tan((2i+1)/2 π) for integers i | ||
| num_discontinuities = Int(div(z, π, RoundNearest)) - Int(div(zₗ, π, RoundNearest)) | ||
| return p + num_discontinuities // k | ||
| end | ||
|
|
||
| logcdf(d::Wrapped{<:Cauchy}, x::Real; kwargs...) = log(cdf(d, x; kwargs...)) | ||
|
|
Oops, something went wrong.
Oops, something went wrong.
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This seems like an odd thing to store in a field and expose to the user. Naively, it seems like an implementation detail that an average user would not need to necessarily care about or modify. I would expect that whether the characteristic function is used would be up to methods defined for particular
Ds inWrapped{D}.There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I agree that this design seems non-ideal. But this field isn't so much for the user as it is for a developer who has some some analysis of
Wrapped{D}for some specificD. In short, the two series representations of pdf have different convergence properties depending onDand its parameters. The one that sums the shifted pdf, etc ofDconverges quickly when the tails of the distribution quickly decay. The one using the characteristic function decays quickly whenDis not very peaked. These conditions are often opposites.Ideally, this would be automatically chosen for every distribution based on its parameters, but this requires someone to sit down and work out bounds for the number of terms needed for each series representation. Having a field allows this to be done at construction time instead of every time
pdfis called. See e.g.Wrapped{Normal}, where we use known bounds to switch between them.For cases where bounds are not known but the user knows one converges faster than the other for their specific
Dand its parameters, the field can be useful, but this is not the primary use.There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Okay, that makes sense, thanks. Would it be reasonable then, do you think, to not document the field and remove it altogether once someone has gone through to work out the numbers of terms?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'm not sure I understand the question. There is no one number of terms. For a given wrapped distribution
distwith values ofparams, one or both of the series will have fewer terms. But if you change any of theparamsvalues or change thedist, then the other series might now have fewer terms. The difference can be quite extreme. For wrapped normal, for example, ifσ=1e+3, the series usingcfneeds only 1 term, while the series using pdf needs >2,000 terms. But forσ=1e-3, the series usingcfneeds >6,000 terms, while the one usingpdfneeds 1 term. So we need to choose whether to use thepdfseries orcfseries. Options:wrapped(::SomeDistribution, ...)if they have a heuristic for choosing the strategypdf,cdf, etc. Likely wastefulcfis not implemented, and does twice as much works as might be necessary.There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Sorry, I was referring to what you had said in the preceding message: