From efa297fed235be13f8ec1d4de174a91d8685d980 Mon Sep 17 00:00:00 2001 From: Tijs Alleman Date: Thu, 26 Dec 2024 21:59:29 +0100 Subject: [PATCH 1/3] add beta distribution, remove weibull distribution, harmonize docstrings --- .../optimization/objective_functions.py | 124 ++++++++---------- 1 file changed, 58 insertions(+), 66 deletions(-) diff --git a/src/pySODM/optimization/objective_functions.py b/src/pySODM/optimization/objective_functions.py index 2c35536..18eee64 100644 --- a/src/pySODM/optimization/objective_functions.py +++ b/src/pySODM/optimization/objective_functions.py @@ -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 @@ -133,130 +133,122 @@ 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 triangle 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=20, density=True) # `samples['my_par']` is a list containing samples from a distribution + + density_my_par_norm = density_my_par/np.sum(density_my_par) # normalise to probability + + 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 ## From f0097b34dd47cc529de7ebb6be6cc9f6aee63eb7 Mon Sep 17 00:00:00 2001 From: Tijs Alleman Date: Thu, 26 Dec 2024 22:16:30 +0100 Subject: [PATCH 2/3] update the docs --- docs/optimization.md | 70 ++++++++++++++++++++++---------------------- 1 file changed, 35 insertions(+), 35 deletions(-) diff --git a/docs/optimization.md b/docs/optimization.md index a99592a..14c2f31 100644 --- a/docs/optimization.md +++ b/docs/optimization.md @@ -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 From 5293eba42280e890887ddec28aa493d005cde9e7 Mon Sep 17 00:00:00 2001 From: Tijs Alleman Date: Thu, 26 Dec 2024 22:17:12 +0100 Subject: [PATCH 3/3] some minor stuffs --- src/pySODM/optimization/objective_functions.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/pySODM/optimization/objective_functions.py b/src/pySODM/optimization/objective_functions.py index 18eee64..7a6d7c4 100644 --- a/src/pySODM/optimization/objective_functions.py +++ b/src/pySODM/optimization/objective_functions.py @@ -163,7 +163,7 @@ def log_prior_triangle(x, args): Returns ------- - Log probability of sample x in light of a triangle prior distribution. + Log probability of sample x in light of a triangular prior distribution. """ low, high, mode = args return triang.logpdf(x, loc=low, scale=high, c=mode) @@ -235,9 +235,7 @@ def log_prior_custom(x, args): Example use: ------------ - density_my_par, bins_my_par = np.histogram(samples['my_par'], bins=20, density=True) # `samples['my_par']` is a list containing samples from a distribution - - density_my_par_norm = density_my_par/np.sum(density_my_par) # normalise to probability + 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