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.

Sign up for GitHub

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

Overhaul log prior distributions #118

Merged
merged 3 commits into from
Dec 26, 2024
+91 −101
Merged

Overhaul log prior distributions #118

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
efa297f
add beta distribution, remove weibull distribution, harmonize docstrings
twallema Dec 26, 2024
f0097b3
update the docs
twallema Dec 26, 2024
5293eba
some minor stuffs
twallema Dec 26, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
70 changes: 35 additions & 35 deletions docs/optimization.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -86,75 +86,75 @@

> **Parameters:**

> * **x** (float) - Parameter value whos probability we want to test.
> * **bounds** (tuple) - Contains the upper and lower bounds of the parameter value.
> * **x** (float) - Parameter value.
> * **bounds** (tuple) - Lower and upper bound of the uniform probability distribution.

> **Returns:**
> * **lp** (float) Log probability of sample x in light of a uniform prior distribution.
> * **lp** (float) Log probability of x in light of a uniform prior distribution.

***function* log_prior_custom(x, args)**
***function* log_prior_triangle(x, pars)**

> Computes the probability of a sample in light of a list containing samples.
> Triangular log prior distribution.

> **Parameters:**
> * **x** (float) - Parameter value whos probability we want to test.
> * **args** (tuple) - Contains the density of each bin in the first position and the bounds of the bins in the second position. Contains a weight given to the custom prior in the third position of the tuple.
> * **x** (float) - Parameter value.
> * **pars** (tuple) - Tuple containg the lower bound, upper bound and mode of the triangle distribution.

> **Returns:**
> * **lp** (float) Log probability of sample x in light of a custom prior distribution.

> **Example use:**

>```python
>density_my_par, bins_my_par = np.histogram([sample_0, sample_1, ..., sample_n], bins=20, density=True)
>density_my_par_norm = density_my_par/np.sum(density_my_par)
>prior_fcn = prior_custom
>prior_fcn_args = (density_my_par_norm, bins_my_par, weight)
>```
> * **lp** (float) Log probability of sample x in light of a triangular prior distribution.

***function* log_prior_normal(x, norm_pars)**
***function* log_prior_normal(x, pars)**

> Normal log prior distribution.

> **Parameters:**
> * **x** (float) - Parameter value whos probability we want to test.
> * **norm_pars** (tuple) - Tuple containg average and standard deviation of normal distribution.
> * **x** (float) - Parameter value.
> * **pars** (tuple) - Tuple containg the average and standard deviation of a normal distribution.

> **Returns:**
> * **lp** (float) Log probability of sample x in light of a normal prior distribution.

***function* log_prior_triangle(x, triangle_pars)**
***function* log_prior_gamma(x, pars)**

> Triangular log prior distribution.
> Gamma log prior distribution.

> **Parameters:**
> * **x** (float) - Parameter value whos probability we want to test.
> * **triangle_pars** (tuple) - Tuple containg lower bound, upper bound and mode of the triangle distribution.
> * **x** (float) - Parameter value.
> * **pars** (tuple) - Tuple containg the parameters `a`, `loc` and `scale` of `scipy.stats.gamma.logpdf`.

> **Returns:**
> * **lp** (float) Log probability of sample x in light of a triangular prior distribution.
> * **lp** (float) Log probability of sample x in light of a gamma prior distribution.

***function* log_prior_gamma(x, gamma_pars)**
***function* log_prior_beta(x, pars)**

> Gamma log prior distribution.
> Beta log prior distribution.

> **Parameters:**
> * **x** (float) - Parameter value whos probability we want to test.
> * **gamma_pars** (tuple) - Tuple containg gamma parameters alpha and beta.
> * **x** (float) - Parameter value.
> * **pars** (tuple) - Tuple containg the parameters `a`, `b`, `loc` and `scale` of `scipy.stats.beta.logpdf`.

> **Returns:**
> * **lp** (float) Log probability of sample x in light of a gamma prior distribution.
> * **lp** (float) Log probability of sample x in light of a beta prior distribution.

***function* log_prior_weibull(x, weibull_params)**

> Weibull log prior distribution.
***function* log_prior_custom(x, args)**

> A custom log prior distribution: compute the probability of a sample in light of a list containing samples from a distribution

> **Parameters:**
> * **x** (float) - Parameter value whos probability we want to test.
> * **weibull_params** (tuple) - Contains the weibull parameters k and lambda.
> * **x** (float) - Parameter value.
> * **args** (tuple) - Must contain the density of each bin in the first position and the bounds of the bins in the second position.

> **Returns:**
> * **lp** (float) Log probability of sample x in light of a weibull prior distribution.
> * **lp** (float) Log probability of x in light of a custom distribution of data.

> **Example use:**

>```python
>density_my_par, bins_my_par = np.histogram([sample_0, sample_1, ..., sample_n], bins=50, density=True) # convert to a list of samples to a binned PDF
>prior_fcn = prior_custom
>prior_fcn_args = (density_my_par_norm, bins_my_par, weight)
>```

## nelder_mead.py

Expand Down
122 changes: 56 additions & 66 deletions src/pySODM/optimization/objective_functions.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,8 +4,8 @@
import numpy as np
import xarray as xr
from datetime import datetime
from scipy.stats import norm, triang, gamma
from scipy.special import gammaln
from scipy.stats import norm, triang, gamma, beta
from pySODM.models.utils import list_to_dict
from pySODM.models.validation import validate_initial_states

Expand Down Expand Up @@ -133,130 +133,120 @@ def ll_negative_binomial(ymodel, ydata, alpha):
#####################################

def log_prior_uniform(x, bounds):
""" Uniform log prior distribution
""" A uniform log prior distribution

Parameters
----------
x: float
Parameter value whos probability we want to test.
Parameter value.
bounds: tuple
Tuple containg the upper and lower bounds of the parameter value.
Lower and upper bound of the uniform probability distribution.

Returns
-------
Log probability of sample x in light of a uniform prior distribution.

Log probability of x in light of a uniform prior distribution.
"""
prob = 1/(bounds[1]-bounds[0])
condition = bounds[0] <= x <= bounds[1]
if condition == True:
# should be set to zero to accomodate bounds = np.inf --> prob is inf!!!
return 0
if bounds[0] <= x <= bounds[1]:
return 0 # technically prob = 1/(upper - lower)
else:
return -np.inf

def log_prior_custom(x, args):
""" Custom log prior distribution: computes the probability of a sample in light of a list containing samples from a previous MCMC run
def log_prior_triangle(x, args):
""" A triangular log prior distribution

Parameters
----------
x: float
Parameter value whos probability we want to test.
Parameter value.
args: tuple
Tuple containg the density of each bin in the first position and the bounds of the bins in the second position.
Contains a weight given to the custom prior in the third position of the tuple.
Tuple containg the lower bound, upper bound and mode of the triangle distribution.

Returns
-------
Log probability of sample x in light of a list with previously sampled parameter values.

Example use:
------------
# Posterior of 'my_par' in samples['my_par']
density_my_par, bins_my_par = np.histogram(samples['my_par'], bins=20, density=True)
density_my_par_norm = density_my_par/np.sum(density_my_par)
prior_fcn = prior_custom
prior_fcn_args = (density_my_par_norm, bins_my_par, weight)
# Prior_fcn and prior_fcn_args must then be passed on to the function log_probability
Log probability of sample x in light of a triangular prior distribution.
"""
density, bins, weight = args
if x <= bins.min() or x >= bins.max():
return -np.inf
else:
idx = np.digitize(x, bins)
return weight*np.log(density[idx-1])
low, high, mode = args
return triang.logpdf(x, loc=low, scale=high, c=mode)

def log_prior_normal(x,norm_params):
""" Normal log prior distribution
def log_prior_normal(x, args):
""" A normal log prior distribution

Parameters
----------
x: float
Parameter value whos probability we want to test.
norm_params: tuple
Tuple containg mu and stdev.
Parameter value.
args: tuple
Tuple containg the average and standard deviation of a normal distribution.

Returns
-------
Log probability of sample x in light of a normal prior distribution.

"""
mu,stdev=norm_params
return np.sum(norm.logpdf(x, loc = mu, scale = stdev))
mu, stdev = args
return np.sum(norm.logpdf(x, loc=mu, scale=stdev))

def log_prior_triangle(x,triangle_params):
""" Triangle log prior distribution
def log_prior_gamma(x, args):
""" A gamma distributed log prior distribution

Parameters
----------
x: float
Parameter value whos probability we want to test.
triangle_params: tuple
Tuple containg lower bound, upper bound and mode of the triangle distribution.
Parameter value.
args: tuple
Tuple containg the parameters `a`, `loc` and `scale` of `scipy.stats.gamma.logpdf`.

Returns
-------
Log probability of sample x in light of a triangle prior distribution.

Log probability of sample x in light of a gamma prior distribution.
"""
low,high,mode = triangle_params
return triang.logpdf(x, loc=low, scale=high, c=mode)
a, loc, scale = args
return gamma.logpdf(x, a=a, loc=loc, scale=scale)

def log_prior_gamma(x,gamma_params):
""" Gamma log prior distribution
def log_prior_beta(x, args):
""" A beta distributed log prior distribution

Parameters
----------
x: float
Parameter value whos probability we want to test.
gamma_params: tuple
Tuple containg gamma parameters alpha, beta and loc (shift along x-axis).
Parameter value.
args: tuple
Tuple containg the parameters `a`, `b`, `loc` and `scale` of `scipy.stats.beta.logpdf`.

Returns
-------
Log probability of sample x in light of a gamma prior distribution.

Log probability of sample x in light of a beta prior distribution.
"""
a,b,loc = gamma_params
return gamma.logpdf(x, a=a, scale=1/b, loc=loc)
a, b, loc, scale = args
return beta.logpdf(x, a, b, loc=loc, scale=scale)

def log_prior_weibull(x,weibull_params):
""" Weibull log prior distribution
def log_prior_custom(x, args):
""" A custom log prior distribution: compute the probability of a sample in light of a list containing samples from a distribution

Parameters
----------
x: float
Parameter value whos probability we want to test.
weibull_params: tuple
Tuple containg weibull parameters k, lambda and loc (shift along x-axis).
Parameter value.
args: tuple
Must contain the density of each bin in the first position and the bounds of the bins in the second position.

Returns
-------
Log probability of sample x in light of a weibull prior distribution.
Log probability of x in light of a custom distribution of data.

Example use:
------------
density_my_par, bins_my_par = np.histogram(samples['my_par'], bins=50, density=True) # a list of samples is converted to a binned PDF

prior_fcn = log_prior_custom # this function

prior_fcn_args = (density_my_par_norm, bins_my_par) # `Prior_fcn` and `prior_fcn_args` must then be passed on to the pySODM `log_probability` class
"""
k,lam,loc = weibull_params
return gamma.logpdf(x, k, shape=lam, loc=loc)
density, bins = args
if x <= bins.min() or x >= bins.max():
return -np.inf
else:
idx = np.digitize(x, bins)
return np.log(density[idx-1])

#############################################
## Computing the log posterior probability ##
Expand Down
Loading