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.

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

Handle Beta degen from Inf params #167

Closed
wants to merge 1 commit into from

Conversation

quildtide
Copy link

Main changes:

  • Handle beta distribution degenerate cases that emerge when alpha, beta, or both are Inf.
    • All behavior is consistent with equivalent functions for the Dirac distribution in Distributions.jl

Additional changes:

  • Implement degenerate case handling for betapdf and betalogpdf
    • Undefined (NaN) at the single point of support, 0 everywhere else
    • This diverges from the Dirac behavior in Distributions.jl, since Dirac is treated as a discrete distribution and pdf(::Dirac) is PMF instead of PDF. PMF is 1, PDF is a value which integrates to 1 at a single point.
  • Quantile functions now always return a constant on degenerate cases
    • previously, quantile functions sometimes returned values that are not in the support of the distribution when the distribution is degenerate.
    • The support of Beta(Inf, 1) does NOT include anything but 1. Docs for quantile both here and in Distributions.jl say that the returned value should be within the support.
    • This is also consistent with how quantile(::Dirac) works in Distributions.jl
  • Additional and adjusted tests to handle all changes in behavior
    • Tests on degenerate beta distributions are performed on a wider range, but with wider increments now (behavior is not going to change between betacdf(0, 0.5, 0.45) and betacdf(0, 0.5, 0.46))

Primary motivation:

I'm fixing handling of degenerate Beta distributions in Distributions.jl right now, and median(Beta(Inf, 1)) hangs because quantile(Beta(Inf, 1), 0.5) hangs.

Between this pull request and one about to be opened in Distributions.jl, degenerate Beta distribution behavior will align with the Dirac distribution behavior, EXCEPT returning pdf = NaN instead of pmf = 1.

@quildtide
Copy link
Author

JuliaStats/Distributions.jl#1881 is the Distributions.jl-side counterpart.

@andreasnoack
Copy link
Member

andreasnoack commented Aug 8, 2024

I think the current behavior (NaN) is reasonable. You don't know where the point with mass will be located. It completely depends on how you reach the limit. If you parameterize with mode and concentration and let the concentration to go infinity then both α and β go to infinity but the location of the mass point will depend on the mode parameter. Usually when the limit is indeterminate like that then you'd have to return NaN.

@quildtide
Copy link
Author

So there's a total of 5 edge cases where the Beta distribution turns into a Dirac. StatsFuns.jl currently handles 2 of them, while Distributions.jl allows the other 3 (but does not handle them correctly).

Letting c be a finite positive value, there are:

  • Beta(0, c) = Dirac(0), handled by StatsFuns.jl already
  • Beta(c, 0) = Dirac(1), handled by StatsFuns.jl already
  • Beta(Inf, c) = Dirac(1), handled by this pull request
  • Beta(c, Inf) = Dirac(0), handled by this pull request
  • Beta(Inf, Inf) = Dirac(?), handled by this pull request

Even if we cannot come to an agreement on Beta(Inf, Inf), I think the other two are worth supporting.

Admittedly, I couldn't find any source for Beta(Inf, Inf) = Dirac(0.5) aside from Wikipedia. I wound up checking the R implementation, and I found that they handle these edge cases consistently:

> rbeta(2, Inf, Inf)
[1] 0.5 0.5
> rbeta(2, Inf, 1)
[1] 1 1
> rbeta(2, 1, Inf)
[1] 0 0
> rbeta(2, 0, 1)
[1] 0 0
> rbeta(2, 1, 0)
[1] 1 1

I don't entirely agree with their implementation, since pbeta(.5, Inf, Inf) returns 1 when I think the value is more of a NaN, and dbeta(0, Inf, 0) returns 0 when I feel that 0 is not in the support.

R actually supports another edge case, Beta(0, 0), or Haldane's Prior, which converges to Bernoulli(0.5).

R's documentation actually calls out:

pl.beta(0, 0)   ## point masses at  {0, 1}

pl.beta(0, 2)   ## point mass at 0 ; the same as
pl.beta(1, Inf)

pl.beta(Inf, 2) ## point mass at 1 ; the same as
pl.beta(3, 0)

pl.beta(Inf, Inf)# point mass at 1/2

The Beta(Inf, Inf) = Dirac(0.5) idea admittedly seems to be dependent on the constraint that α and β approach infinity at the same rate. I suppose that loosening that constraint (e.g. α = 2β) would allow you to create arbitrary Dirac(α/(α+β )) from Beta(Inf, Inf).

@quildtide
Copy link
Author

The way pdf(Normal(mu, 0), mu) is handled here and in Distributions.jl is to return Inf, not NaN. I can see valid arguments for both.

@quildtide
Copy link
Author

Closing the companion pull request in Distributions.jl with intent of an alternative pull request that just blocks creation of Beta distributions with Inf parameters.

@quildtide quildtide closed this Aug 9, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants