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

Constructors for Gaussian priors with constraints #161

Merged
merged 1 commit into from
Jun 30, 2022

Conversation

tsj5
Copy link
Collaborator

@tsj5 tsj5 commented May 17, 2022

It's desirable to streamline the UI for specifying parameter priors, in particular those for parameters that live on a constrained domain and hence need a nontrivial Constraint().

This PR does this by adding a new constructor, bounded_gaussian constrained_gaussian, for ParameterizedDistribution. The unconstrained distribution is a Gaussian, with parameters (μ_u, σ_u) selected to reproduce user-required statistics in the constrained space [currently only (mean, standard deviation) are supported]. This therefore includes the lognormal and logit-normal distributions as special cases.

In the cases where the constraint interval is infinite or semi-infinite, closed-form expressions are used for (μ_u, σ_u). When the interval is finite, no closed-form solution exists, and (μ_u, σ_u) are optimized numerically via the simplex method. This is still less reliable than it should be for large |μ_u| or σ_u: for these cases the distribution becomes concentrated near one or both interval endpoints, causing instabilities for both numerical integration and optimization (since parameter changes have minimal effect on the objective function).

This PR fixes a sign error in bounded_above() which caused its Jacobian to be negative.

Tasks remaining as of this writing:

  • constructor for distributions given as Cartesian products of the 1D case currently implemented.
  • interface for using constrained_gaussian to specify priors from parameter file conventions.

@tsj5
Copy link
Collaborator Author

tsj5 commented May 17, 2022

Does anyone know how to fix the failing JuliaFormatter/format test? I've run julia .dev/climaformat.jl locally, and running it again doesn't change any files. Is the version of JuliaFormatter that's run in CI using different settings than climaformat.jl?

@tsj5 tsj5 force-pushed the tsj/prior-init-2 branch 2 times, most recently from 91e0b5d to dffb677 Compare May 17, 2022 21:10
Copy link
Collaborator

@odunbar odunbar left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi Tom - this looks much more along the lines of what I was looking for.
In the todo's, could you also add something to the docs to show how it works?

src/ParameterDistributions.jl Outdated Show resolved Hide resolved
src/ParameterDistributions.jl Outdated Show resolved Hide resolved
if ~isapprox(μ_c, m_c, atol = 1e-3, rtol = 1e-2)
@warn "Unable to set constrained mean for `$(name)`: target = $(μ_c), got $(m_c)"
end
if ~isapprox(σ_c, s_c, atol = 1e-3, rtol = 1e-2)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If these are taken in the constrained space, perhaps they should be normalized? (or is that what rtol does?)

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry, I'm lost on the proper context for "normalized" -- normalization, svd. etc. are only done in CES, right?

m_c and s_c are the mean/stddev of the distribution at its final, hopefully optimized parameter values, computed as 1d numerical integrals in the physical, constrained space. The highlighted lines are a sanity check to verify that the optimization set those parameters to values that satisfy the user's requested μ_c, σ_c for the prior, to very loose tolerances. rtol is just the tolerance on the relative difference between the two values.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry i was just asking about tolerances. Namely why do we need atol here? (do both tolerances needed to be satisfied?) Because the constrained space is not normalized so if your domain is on [0,1e-5] or [0, 1e5] in the first case atol is satisfied always, and in the second, one could imagine it would never be satisfied

# optimize in log σ to avoid constraint σ>0
init_logσ_u = log(init_σ_u)

# Optimize parameters; could speed up with gradient info but don't bother to.
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If we don't need gradients then they can be removed from the constraints?

@tsj5
Copy link
Collaborator Author

tsj5 commented May 31, 2022

Resolution of the climaformat issue mentioned above, written up here for discoverability: @odunbar mentioned that the version of JuliaFormatter had been updated to 0.22; this was in .dev/Project.toml but for some reason .dev/Manifest.toml was still using the older version (I would have thought that that would've been caught and updated when that environment was activate'ed). Deleting the Manifest and re-instantiating solved the issue.

@tsj5
Copy link
Collaborator Author

tsj5 commented Jun 2, 2022

bors try

bors bot added a commit that referenced this pull request Jun 2, 2022
@bors
Copy link
Contributor

bors bot commented Jun 2, 2022

try

Build succeeded:

Copy link
Collaborator

@odunbar odunbar left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As this is user facing - could you also provide an example snippet and some information in the docs priors page?

Otherwise, LGTM!

@tsj5
Copy link
Collaborator Author

tsj5 commented Jun 9, 2022

Finally added docs and examples for constrained_gaussian. See GitHub's markdown preview here -- I don't think there's any way to trigger the docs build for PRs, right?

Copy link
Collaborator

@odunbar odunbar left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry it's taken a while for review - this looks great. Please merge after addressing my comments.

  • My only real point is that now we have a "recommended" constructor, we should present this first, likewise with the simple example for it, so could you move everything for this up to the top? Then most users will only have to read the first part of this page. After we can talk about general constructors, and all the details for the constraints etc.

docs/src/parameter_distributions.md Outdated Show resolved Hide resolved
docs/src/parameter_distributions.md Outdated Show resolved Hide resolved
docs/src/parameter_distributions.md Show resolved Hide resolved
docs/src/parameter_distributions.md Outdated Show resolved Hide resolved
docs/src/parameter_distributions.md Outdated Show resolved Hide resolved
docs/src/parameter_distributions.md Outdated Show resolved Hide resolved
docs/src/parameter_distributions.md Outdated Show resolved Hide resolved
@tsj5
Copy link
Collaborator Author

tsj5 commented Jun 28, 2022

Sorry it's taken a while for review - this looks great. Please merge after addressing my comments.

  • My only real point is that now we have a "recommended" constructor, we should present this first, likewise with the simple example for it, so could you move everything for this up to the top? Then most users will only have to read the first part of this page. After we can talk about general constructors, and all the details for the constraints etc.

Good idea -- I definitely agree that it'll be best for users if we put all the info pertaining to the "recommended" constructor up front.

@tsj5
Copy link
Collaborator Author

tsj5 commented Jun 28, 2022

Current CI failures are described by issue #175 and fixed in PR #176.

@odunbar
Copy link
Collaborator

odunbar commented Jun 28, 2022

@tsj5 Much prefer the structure now! LGTM - merge away!

PS - you have been extraordinarily unlucky with the CI failures on your PRs :(

@tsj5
Copy link
Collaborator Author

tsj5 commented Jun 28, 2022

bors try

bors bot added a commit that referenced this pull request Jun 28, 2022
@bors
Copy link
Contributor

bors bot commented Jun 28, 2022

try

Build succeeded:

Add Optim, QuadGK as deps

Add xform jacobian to Constraint struct

Properly handle +/-Inf in Constraints

Fix sign of bounded_above()

Initial commit of bounded_gaussian

Rename bounded_gaussian -> constrained_gaussian

Extra tests for codecov

Add optional kwargs passed to Optim.optimize

Reorder sections of parameter_distributions.md doc

Add docs & examples for constrained_gaussian()
@tsj5
Copy link
Collaborator Author

tsj5 commented Jun 29, 2022

bors try

bors bot added a commit that referenced this pull request Jun 29, 2022
@tsj5 tsj5 changed the title [WIP] Constructors for priors with constraints Constructors for Gaussian priors with constraints Jun 29, 2022
@bors
Copy link
Contributor

bors bot commented Jun 29, 2022

try

Build succeeded:

@tsj5
Copy link
Collaborator Author

tsj5 commented Jun 29, 2022

bors r+

bors bot added a commit that referenced this pull request Jun 29, 2022
161: Constructors for Gaussian priors with constraints r=tsj5 a=tsj5

It's desirable to streamline the UI for specifying parameter priors, in particular those for parameters that live on a constrained domain and hence need a nontrivial Constraint().

This PR does this by adding a new constructor, ~~bounded_gaussian~~ `constrained_gaussian`, for ParameterizedDistribution. The unconstrained distribution is a Gaussian, with parameters (μ_u, σ_u) selected to reproduce user-required statistics in the constrained space [currently only (mean, standard deviation) are supported]. This therefore includes the [lognormal](https://en.wikipedia.org/wiki/Log-normal_distribution) and [logit-normal](https://en.wikipedia.org/wiki/Logit-normal_distribution) distributions as special cases.

In the cases where the constraint interval is infinite or semi-infinite, closed-form expressions are used for (μ_u, σ_u). When the interval is finite, no closed-form solution exists, and (μ_u, σ_u) are optimized numerically via the simplex method. This is still less reliable than it should be for large |μ_u| or σ_u: for these cases the distribution becomes concentrated near one or both interval endpoints, causing instabilities for both numerical integration and optimization (since parameter changes have minimal effect on the objective function).

This PR fixes a sign error in bounded_above() which caused its Jacobian to be negative.

Tasks remaining as of this writing:
- constructor for distributions given as Cartesian products of the 1D case currently implemented.
- interface for using `constrained_gaussian` to specify priors from parameter file conventions.


Co-authored-by: Thomas Jackson <[email protected]>
@bors
Copy link
Contributor

bors bot commented Jun 29, 2022

Build failed:

@tsj5
Copy link
Collaborator Author

tsj5 commented Jun 30, 2022

bors r+

@bors
Copy link
Contributor

bors bot commented Jun 30, 2022

Build succeeded:

@bors bors bot merged commit 6585e00 into CliMA:main Jun 30, 2022
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