Folded Distributions #1631

Potatoasad · 2022-10-26T04:08:02Z

Folded Distributions

Here's a Generic implementation for Folded continuous univariate distributions. I worked on this a while ago and only recently cleaned it up. This has been useful to me in my research and hopefully is helpful to others as well.

Definition

A folded distribution $F_{c}(D)$ of a distribution $D$ at crease $c \in \mathbb{R}$ is the reflection of the distribution below $c$ onto the the distribution above $c$.

The pdf of such a distribution is given by:

$$ p(x | F_{c}(D)) = p(x | D) + p(x' | D) $$

where $x' = 2c - x$

Interface

The interface is very very similar to the truncated interface.
You can take any continuous univariate distribution and make it into a folded distribution like so:

julia> normal_absolute = folded(Normal(2.0,1.0), 0.0)
# Folded(Normal{Float64}(μ=2.0, σ=1.0); crease = 0.0 | fold left onto right)
julia> plot(x -> pdf(normal_absolute, x))

This gives the result:

The folded function creates a folded distribution from the original distribution d about the value crease.
This function defaults to folding the left side onto the right, but by using keep_right=false one can reflect the right side onto the left.

This can be implemented for any univariate continuous distribution by using the following method:

folded(d::ContinuousUnivariateDistribution, crease::Real)

If one wants to reflect points above the crease $c$ onto points below $c$ (such that the resultant distribution lives on points below $c$) one can do that using the keep_right=false tag:

folded(d::ContinuousUnivariateDistribution, crease::Real, keep_right=false)

A very useful and oft-occuring example of this is the action of the absolute value function ($|\cdot|$) on random variables.
For example if we have a random variable obeying the normal distribution:
$$\hat{X} \sim \mathcal{N}(\alpha,\beta)$$

Taking the absolute value of this random variable always follows a folded normal distribution with a crease at $0$

$$|\hat{X}| \sim F_{c=0}(\mathcal{N}(\alpha,\beta))$$
In julia that looks like:

folded_normal = folded(Normal(μ,σ), 0.0)

and then we can ask for many things like the pdf, or sample the distribution. The samples of course will then be positive:

all(rand(folded_normal, 1000) .≥ 0) # True

In general, one can take any univariate continous distribution and fold it using this function.

Example

We can create the absolute value of the laplace distribution by:

folded_laplace = folded(Laplace(μ,b), 0.0)

as a quick check we can see if the 1000 samples generated from this are all greater than 0:

all(rand(folded(Lapace(μ,b),0),1000) .≥ 0) # -> true

For more one can check the docs in the PR.

With this addition one should be able to define stuff like the distribution of the absolute value of any univariate random variable.

Distributions this could lead to

In issue #124 ,there was a desire to implement the following distributions:

foldcauchy - A folded Cauchy continuous random variable
foldnorm - A folded normal continuous random variable
halfcauchy - A Half-Cauchy continuous random variable
halflogistic - A half-logistic continuous random variable
halfnorm - A half-normal continuous random variable
and all of these can be implemented very easily using the folded function.

For example, an implementation of the HalfNormal distribution, would be as simple as

halfnorm(σ) = folded(Normal(zero(σ), σ), zero(σ))

The others can be implemented similarly very easily. The Folded generic subsumes all these cases.

Conclusion & Why would this be useful?

I've added tests and docs, and wanted to see if this is an appropriate inclusion to the package.
Let me know if this is would be welcomed by this package, and if there's anything I should add to the tests or the docs, or change the code style.

Thank you for all the work you all do on this amazing package!

Added all the functions necessary to make Folded Distributions work for univariate continuous distributions. Need to figure out how to make this work for discrete distributions

…se is finalized

azev77 · 2022-10-28T15:17:24Z

This is a brilliant incredibly useful idea!!!

Potatoasad · 2022-10-29T23:20:11Z

This is a brilliant incredibly useful idea!!!

Thank you so much! Hoping this can get merged soon 🤞

Potatoasad · 2022-12-30T20:41:32Z

Hi! Just wanted to get an update on this if this is good enough to merge - I think it's ready for review from my end.
I understand the maintainers have a limited time so no worries at all, just wanted to bump this.

sethaxen · 2022-12-30T21:07:01Z

Hey @Potatoasad, I had just hacked together a Folded distribution for my own research and just now saw your PR. I'd be happy to review it and will try to do so this weekend.

ParadaCarleton · 2022-12-30T23:02:11Z

src/fold.jl

+```
+as a quick check we can see if the 1000 samples generated from this are all greater than 0:
+```julia
+all(rand(folded(Lapace(μ,b),0),1000) .≥ 0) # -> true


Suggested change

all(rand(folded(Lapace(μ,b),0),1000) .≥ 0) # -> true

all(≥(0), rand(folded_laplace, 1000)) # -> true

I think this is easier to read; it's also more efficient, since it avoids allocating an intermediate array.

ParadaCarleton · 2022-12-30T23:03:31Z

src/fold.jl

@@ -0,0 +1,320 @@
+"""
+    folded(d::ContinuousUnivariateDistribution, crease::Real)


Suggested change

folded(d::ContinuousUnivariateDistribution, crease::Real)

Not necessary to document the no-keyword method, since the keyword method already covers it (by showing the default value).

ParadaCarleton · 2022-12-30T23:43:40Z

src/fold.jl

+
+
+##### Defining the function and the initialization functions #####
+"""


Should follow the Blue style guide for docstrings--see here.

No. JuliaStats does not use Blue style. One should just be consistent with the existing style until the org decides on a specific style. It does not matter anyway as all these docstrings should be removed (see my review).

ParadaCarleton · 2022-12-30T23:47:23Z

src/fold.jl

+"""
+struct Folded{D<:ContinuousUnivariateDistribution, K<:Truncated, S<:Continuous, T <: Real} <: ContinuousUnivariateDistribution
+    original::D      # the original distribution (unfolded)
+    included::K     # The part of the old distribution left unchanged


Not sure about the names here. I feel like I'd have a hard time remembering included/excluded 😅

ParadaCarleton · 2022-12-30T23:49:06Z

src/fold.jl

+    folded(d::ContinuousUnivariateDistribution, crease::Real)
+    folded(d::ContinuousUnivariateDistribution, crease::Real, keep_right=false)
+
+Creates a _folded distribution_ from the original distribution `d` about the value `crease`.


Suggested change

Creates a _folded distribution_ from the original distribution `d` about the value `crease`.

Fold a distribution `d` at the value `crease`.

ParadaCarleton · 2022-12-31T02:13:49Z

src/fold.jl

+# >> MethodError: no method matching iterate(::Truncated{Beta{...}...})
+```
+"""
+mean(d::Folded) = mean(d.included)*d.included.tp - mean(d.excluded)*d.excluded.tp + 2*d.crease


Can you add comments or expand this out? It's not immediately obvious to me why this is correct.

ParadaCarleton · 2022-12-31T02:24:11Z

test/fold.jl

+end 
+
+# Testing if the insupport function works
+@testset "Folded insupport, minimum and maximum functions" begin


Can you add more tests? We should be covering at least mean, var, median, and at least two different quantiles.

ParadaCarleton · 2022-12-31T02:25:10Z

test/fold.jl

+	@test all(((x->Distributions.maximum(folded(Normal(2,1),x,keep_right=false))).(-10:0.1:10)) .≈ collect(-10:0.1:10))
+end
+
+@testset "Folded Rand functions" begin


Same as for the above--we need tests that check if rand is working correctly, e.g. testing if a big sample from a folded normal distribution has the right distribution.

ParadaCarleton · 2022-12-31T02:32:19Z

test/fold.jl

+@testset "Folded pdf, logpdf, cdf and logcdf functions" begin
+	@test pdf(truncated(Normal(3,4),-1,1),-0.1) + pdf(truncated(Normal(3,4),-1,1),0.1) ≈ pdf(folded(truncated(Normal(3,4),-1,1),0.0),0.1)
+	@test pdf(folded(truncated(Normal(3,4),-1,1),0.0),-0.1) == 0
+end


Suggested change

@testset "Folded pdf, logpdf, cdf and logcdf functions" begin

@test pdf(truncated(Normal(3,4),-1,1),-0.1) + pdf(truncated(Normal(3,4),-1,1),0.1) ≈ pdf(folded(truncated(Normal(3,4),-1,1),0.0),0.1)

@test pdf(folded(truncated(Normal(3,4),-1,1),0.0),-0.1) == 0

end

@testset "Folded pdf, logpdf, cdf and logcdf functions" begin

trunc_dist = truncated(Normal(3, 4), -1, 1)

fold_dist = folded(trunc_dist, 0.0)

@test pdf(truncated(trunc_dist, -0.1) + pdf(trunc_dist, 0.1) ≈ pdf(fold_dist, 0.1)

@test pdf(fold_dist, -0.1) ≈ 0

end

ParadaCarleton · 2022-12-31T02:33:44Z

test/fold.jl

+	@test all(rand(folded(Normal(1.0,1.0),0.0),100) .> 0.0)
+end
+
+@testset "Folded pdf, logpdf, cdf and logcdf functions" begin


Can you add tests that deal with cases other than folding at 0?

devmotion

Thank you for the PR! I added some initial comments.

src/Distributions.jl

devmotion · 2022-12-30T21:40:10Z

src/fold.jl

+
+# This is the user facing initialization function that will be used in most cases
+"""
+   folded(d::ContinuousUnivariateDistribution, crease::T; keep_right=true) 


Suggested change

folded(d::ContinuousUnivariateDistribution, crease::T; keep_right=true)

folded(d::ContinuousUnivariateDistribution, crease::Real; keep_right::Bool=true)

src/fold.jl

devmotion · 2022-12-30T22:39:52Z

src/fold.jl

@@ -0,0 +1,320 @@
+"""
+    folded(d::ContinuousUnivariateDistribution, crease::Real)
+    folded(d::ContinuousUnivariateDistribution, crease::Real, keep_right=false)


Suggested change

folded(d::ContinuousUnivariateDistribution, crease::Real, keep_right=false)

folded(d::ContinuousUnivariateDistribution, crease::Real; keep_right::Bool=true)

Maybe a (better?) name for the keyword argument would be

Suggested change

folded(d::ContinuousUnivariateDistribution, crease::Real, keep_right=false)

folded(d::ContinuousUnivariateDistribution, crease::Real; above::Bool=true)

or

Suggested change

folded(d::ContinuousUnivariateDistribution, crease::Real, keep_right=false)

folded(d::ContinuousUnivariateDistribution, crease::Real; support_above::Bool=true)

devmotion · 2022-12-30T22:44:00Z

src/fold.jl

+Holds the original distribution, and two truncated copies of the 
+distribution: one above the crease and one below. 
+"""
+struct Folded{D<:ContinuousUnivariateDistribution, K<:Truncated, S<:Continuous, T <: Real} <: ContinuousUnivariateDistribution


Why do you want to restrict K to subtypes of Truncated? S is not needed?

Suggested change

struct Folded{D<:ContinuousUnivariateDistribution, K<:Truncated, S<:Continuous, T <: Real} <: ContinuousUnivariateDistribution

struct Folded{D<:ContinuousUnivariateDistribution, K, T <: Real} <: ContinuousUnivariateDistribution

devmotion · 2022-12-30T23:06:08Z

src/fold.jl

+"""
+function cdf(d::Folded, x::Real) 
+    if !insupport(d,x)
+        (x ≤ minimum(d)) && return zero(x) # Is below the support


Again type unstable. Again best to always compute the result (possibly with adjusted x) below. Maybe it's already sufficient to just perform the evaluation for clamp(x, extrema(d)...).

devmotion · 2022-12-31T07:30:23Z

src/fold.jl

+        (x ≤ minimum(d)) && return zero(x) # Is below the support
+        (x ≥ maximum(d)) && return one(x) # Is above the support
+    end
+    d.included.tp*cdf(d.included, x) + d.excluded.tp*ccdf(d.excluded, unfold_value(x,d))


d.included.tp does not exist for general truncated distributions. Could we use something like

Suggested change

d.included.tp*cdf(d.included, x) + d.excluded.tp*ccdf(d.excluded, unfold_value(x,d))

cdf(d::Folded, x::Real) = exp(logcdf(d, x)) # not sure, maybe the default definition already?

function logcdf(d::Folded, x::Real)

_x = clamp(x, extrema(d.included)...)

y = unfolded_value(d, _x)

a, b = d.keep_right ? (y, _x) : (_x, y)

return logdiffcdf(d.original, b, a)

end

My understanding, by looking at the source code was that the Truncated object has a field tp in its definition, and having this would prevent recalculation each time.

Is there a way to infer from the types whether a truncated distribution has the tp field? So if it doesn't exist for a particular type, we could dispatch to the code you provided?

Distributions.jl/src/truncate.jl

Line 93 in fa8c30d

struct Truncated{D<:UnivariateDistribution, S<:ValueSupport, T <: Real} <: UnivariateDistribution{S}

devmotion · 2022-12-31T07:31:47Z

src/fold.jl

+    d.included.tp*cdf(d.included, x) + d.excluded.tp*ccdf(d.excluded, unfold_value(x,d))
+end
+
+


ccdf and logccdf are not defined, ie onky fallback definitions will be used?

devmotion · 2022-12-31T07:32:11Z

src/fold.jl

+######################################################
+
+"""
+    mean(d::Folded)


Please remove.

devmotion · 2022-12-31T07:36:04Z

src/fold.jl

+# >> MethodError: no method matching iterate(::Truncated{Beta{...}...})
+```
+"""
+mean(d::Folded) = mean(d.included)*d.included.tp - mean(d.excluded)*d.excluded.tp + 2*d.crease


Again tp does not exist in general.

Potatoasad · 2023-01-01T03:46:09Z

Thank you so much for the detailed review! I'll pour over these and make the changes as needed thank you so much for the feedback and explanations.

azev77 · 2023-05-10T19:45:00Z

Thank you so much for the detailed review! I'll pour over these and make the changes as needed thank you so much for the feedback and explanations.

@Potatoasad have you had any luck with this?

Potatoasad · 2023-05-11T16:00:50Z

Thank you so much for the detailed review! I'll pour over these and make the changes as needed thank you so much for the feedback and explanations.

@Potatoasad have you had any luck with this?

Hey, sorry. Yes I did make some slight progress, but got sidetracked into a paper I was working on. I think I can finish this up in the next two weeks. Apologies for the very long delay on this.

ParadaCarleton · 2023-09-03T16:28:40Z

Thank you so much for the detailed review! I'll pour over these and make the changes as needed thank you so much for the feedback and explanations.

@Potatoasad have you had any luck with this?

Hey, sorry. Yes I did make some slight progress, but got sidetracked into a paper I was working on. I think I can finish this up in the next two weeks. Apologies for the very long delay on this.

Hi, just wanted to know if you're still interested in working on this :)

Co-authored-by: David Widmann <[email protected]>

codecov-commenter · 2023-09-04T22:28:47Z

Codecov Report

Patch coverage is 59.25% of modified lines.

❗ Current head 009f83d differs from pull request most recent head f79504f. Consider uploading reports for the commit f79504f to get more accurate results

Files Changed	Coverage
src/Distributions.jl	`ø`
src/fold.jl	`59.25%`

📢 Thoughts on this report? Let us know!.

Co-authored-by: Carlos Parada <[email protected]>

Potatoasad added 7 commits March 31, 2022 14:47

Added the Folded Distribution and relevant definitions

a3fe4e1

Added all the functions necessary to make Folded Distributions work for univariate continuous distributions. Need to figure out how to make this work for discrete distributions

Specialized everything only for ContinuousUnivariate until general ca…

6cb3eac

…se is finalized

Merge branch 'JuliaStats:master' into folded-generic-continuous-only

d3433f9

use exp. of loglikelihood for numerical stability

5d1a914

Added Docs for folded distributions

3e8575f

rand function for folded distributions now passes the AbstractRNG object

8a2363d

Merge branch 'master' into folded-generic-continuous-only

009f83d

ParadaCarleton requested changes Dec 31, 2022

View reviewed changes

devmotion reviewed Dec 31, 2022

View reviewed changes

Potatoasad and others added 2 commits September 4, 2023 17:27

Update src/fold.jl

6afa096

Co-authored-by: David Widmann <[email protected]>

Update src/Distributions.jl

c54b8f3

Co-authored-by: David Widmann <[email protected]>

Potatoasad and others added 2 commits September 4, 2023 17:29

Update src/fold.jl

b0ec45a

Co-authored-by: Carlos Parada <[email protected]>

Update src/fold.jl

f79504f

Co-authored-by: Carlos Parada <[email protected]>

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Folded Distributions #1631

Folded Distributions #1631

Potatoasad commented Oct 26, 2022 •

edited

Loading

azev77 commented Oct 28, 2022

Potatoasad commented Oct 29, 2022

Potatoasad commented Dec 30, 2022

sethaxen commented Dec 30, 2022

ParadaCarleton Dec 30, 2022

ParadaCarleton Dec 30, 2022

ParadaCarleton Dec 30, 2022

devmotion Jan 1, 2023

ParadaCarleton Dec 30, 2022

ParadaCarleton Dec 30, 2022

ParadaCarleton Dec 31, 2022

ParadaCarleton Dec 31, 2022

ParadaCarleton Dec 31, 2022

ParadaCarleton Dec 31, 2022

ParadaCarleton Dec 31, 2022

devmotion left a comment

devmotion Dec 30, 2022

devmotion Dec 30, 2022

devmotion Dec 30, 2022

devmotion Dec 30, 2022

devmotion Dec 31, 2022

Potatoasad May 11, 2023

devmotion Dec 31, 2022

devmotion Dec 31, 2022

devmotion Dec 31, 2022

Potatoasad commented Jan 1, 2023

azev77 commented May 10, 2023

Potatoasad commented May 11, 2023

ParadaCarleton commented Sep 3, 2023

codecov-commenter commented Sep 4, 2023 •

edited

Loading

	all(rand(folded(Lapace(μ,b),0),1000) .≥ 0) # -> true
	all(≥(0), rand(folded_laplace, 1000)) # -> true

		@@ -0,0 +1,320 @@
		"""
		folded(d::ContinuousUnivariateDistribution, crease::Real)



		##### Defining the function and the initialization functions #####
		"""

	Creates a _folded distribution_ from the original distribution `d` about the value `crease`.
	Fold a distribution `d` at the value `crease`.

	struct Folded{D<:ContinuousUnivariateDistribution, K<:Truncated, S<:Continuous, T <: Real} <: ContinuousUnivariateDistribution
	struct Folded{D<:ContinuousUnivariateDistribution, K, T <: Real} <: ContinuousUnivariateDistribution

-    d.included.tp*cdf(d.included, x) + d.excluded.tp*ccdf(d.excluded, unfold_value(x,d))
+cdf(d::Folded, x::Real) = exp(logcdf(d, x)) # not sure, maybe the default definition already?
+function logcdf(d::Folded, x::Real)
+    _x = clamp(x, extrema(d.included)...)
+    y = unfolded_value(d, _x)
+    a, b = d.keep_right ? (y, _x) : (_x, y)
+    return logdiffcdf(d.original, b, a)
+end

		d.included.tpcdf(d.included, x) + d.excluded.tpccdf(d.excluded, unfold_value(x,d))
		end

Folded Distributions #1631

Are you sure you want to change the base?

Folded Distributions #1631

Conversation

Potatoasad commented Oct 26, 2022 • edited Loading

Folded Distributions

Definition

Interface

Example

Distributions this could lead to

Conclusion & Why would this be useful?

azev77 commented Oct 28, 2022

Potatoasad commented Oct 29, 2022

Potatoasad commented Dec 30, 2022

sethaxen commented Dec 30, 2022

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

devmotion left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Potatoasad commented Jan 1, 2023

azev77 commented May 10, 2023

Potatoasad commented May 11, 2023

ParadaCarleton commented Sep 3, 2023

codecov-commenter commented Sep 4, 2023 • edited Loading

Codecov Report

Potatoasad commented Oct 26, 2022 •

edited

Loading

codecov-commenter commented Sep 4, 2023 •

edited

Loading