From 26db9761ed147bfc76fa7706f70d092c3eb7a106 Mon Sep 17 00:00:00 2001 From: Polina Lakrisenko Date: Wed, 11 Jan 2023 10:06:38 +0100 Subject: [PATCH 1/5] add possibility to provide parameters for noise distributiond; add left-truncated-normal, left-censored-normal --- python/sdist/amici/import_utils.py | 47 +++++++++++++++++++++--------- python/sdist/amici/petab_import.py | 13 +++++++-- python/sdist/amici/sbml_import.py | 28 +++++++++++------- 3 files changed, 60 insertions(+), 28 deletions(-) diff --git a/python/sdist/amici/import_utils.py b/python/sdist/amici/import_utils.py index 65d7285a65..17ff0a01f2 100644 --- a/python/sdist/amici/import_utils.py +++ b/python/sdist/amici/import_utils.py @@ -52,7 +52,7 @@ class ObservableTransformation(str, enum.Enum): def noise_distribution_to_observable_transformation( - noise_distribution: Union[str, Callable] + noise_distribution: Union[Dict[str, Union[str, list]], Callable] ) -> ObservableTransformation: """ Parse noise distribution string and extract observable transformation @@ -63,17 +63,18 @@ def noise_distribution_to_observable_transformation( :return: observable transformation """ - if isinstance(noise_distribution, str): - if noise_distribution.startswith('log-'): + if isinstance(noise_distribution, dict): + noise_distribution_type = noise_distribution['type'] + if noise_distribution_type.startswith('log-'): return ObservableTransformation.LOG - if noise_distribution.startswith('log10-'): + if noise_distribution_type.startswith('log10-'): return ObservableTransformation.LOG10 return ObservableTransformation.LIN def noise_distribution_to_cost_function( - noise_distribution: Union[str, Callable] + noise_distribution: Union[dict, Callable] ) -> Callable[[str], str]: """ Parse noise distribution string to a cost function definition amici can @@ -104,6 +105,10 @@ def noise_distribution_to_cost_function( \\pi(m|y,\\sigma) = \\frac{1}{\\sqrt{2\\pi}\\sigma m \\log(10)}\\ exp\\left(-\\frac{(\\log_{10} m - \\log_{10} y)^2}{2\\sigma^2}\\right) + - `'truncated-normal'`: A truncated normal distribution: + + .. math:: + - `'laplace'`, `'lin-laplace'`: A laplace distribution: .. math:: @@ -181,27 +186,37 @@ def noise_distribution_to_cost_function( if isinstance(noise_distribution, Callable): return noise_distribution - if noise_distribution in ['normal', 'lin-normal']: + noise_distribution_type = noise_distribution['dist_type'] + + if noise_distribution_type in ['normal', 'lin-normal']: y_string = '0.5*log(2*pi*{sigma}**2) + 0.5*(({y} - {m}) / {sigma})**2' - elif noise_distribution == 'log-normal': + elif noise_distribution_type == 'log-normal': y_string = '0.5*log(2*pi*{sigma}**2*{m}**2) ' \ '+ 0.5*((log({y}) - log({m})) / {sigma})**2' - elif noise_distribution == 'log10-normal': + elif noise_distribution_type == 'log10-normal': y_string = '0.5*log(2*pi*{sigma}**2*{m}**2*log(10)**2) ' \ '+ 0.5*((log({y}, 10) - log({m}, 10)) / {sigma})**2' - elif noise_distribution in ['laplace', 'lin-laplace']: + elif noise_distribution_type == ['left-truncated-normal', + 'lin-left-truncated-normal']: + # truncated at a + a = noise_distribution['dist_params'][0] # TODO + y_string = f'log(1-0.5*(1+erf({a}-{{y}}/(sqrt(2)*{{sigma}}))))' \ + ' + 0.5*log(2*pi*{sigma}**2)' \ + ' + 0.5*(({y} - {m}) / {sigma})**2' + elif noise_distribution_type in ['laplace', 'lin-laplace']: y_string = 'log(2*{sigma}) + Abs({y} - {m}) / {sigma}' - elif noise_distribution == 'log-laplace': + elif noise_distribution_type == 'log-laplace': y_string = 'log(2*{sigma}*{m}) + Abs(log({y}) - log({m})) / {sigma}' - elif noise_distribution == 'log10-laplace': + elif noise_distribution_type == 'log10-laplace': y_string = 'log(2*{sigma}*{m}*log(10)) ' \ '+ Abs(log({y}, 10) - log({m}, 10)) / {sigma}' - elif noise_distribution in ['binomial', 'lin-binomial']: + elif noise_distribution_type in ['binomial', 'lin-binomial']: # Binomial noise model parameterized via success probability p y_string = '- log(Heaviside({y} - {m})) - loggamma({y}+1) ' \ '+ loggamma({m}+1) + loggamma({y}-{m}+1) ' \ '- {m} * log({sigma}) - ({y} - {m}) * log(1-{sigma})' - elif noise_distribution in ['negative-binomial', 'lin-negative-binomial']: + elif noise_distribution_type in ['negative-binomial', + 'lin-negative-binomial']: # Negative binomial noise model of the number of successes m # (data) before r=(1-sigma)/sigma * y failures occur, # with mean number of successes y (simulation), @@ -210,9 +225,13 @@ def noise_distribution_to_cost_function( y_string = f'- loggamma({{m}}+{r}) + loggamma({{m}}+1) ' \ f'+ loggamma({r}) - {r} * log(1-{{sigma}}) ' \ f'- {{m}} * log({{sigma}})' + elif noise_distribution_type == 'left-censored-normal': + # left-censored at v (detection limit) + v = noise_distribution['dist_params'][0] # TODO + y_string = f'log(0.5*(1+erf(({v}-{{y}})/(sqrt(2)*{{sigma}}))))' else: raise ValueError( - f"Cost identifier {noise_distribution} not recognized.") + f"Cost identifier {noise_distribution_type} not recognized.") def nllh_y_string(str_symbol): y, m, sigma = _get_str_symbol_identifiers(str_symbol) diff --git a/python/sdist/amici/petab_import.py b/python/sdist/amici/petab_import.py index 0fc9ed251d..850c576ad9 100644 --- a/python/sdist/amici/petab_import.py +++ b/python/sdist/amici/petab_import.py @@ -683,7 +683,7 @@ def import_model_sbml( def get_observation_model( observable_df: pd.DataFrame, -) -> Tuple[Dict[str, Dict[str, str]], Dict[str, str], +) -> Tuple[Dict[str, Dict[str, str]], Dict[str, Dict[str, Union[str, list]]], Dict[str, Union[str, float]]]: """ Get observables, sigmas, and noise distributions from PEtab observation @@ -730,7 +730,8 @@ def get_observation_model( def petab_noise_distributions_to_amici(observable_df: pd.DataFrame - ) -> Dict[str, str]: + ) -> Dict[str, + Dict[str, Union[str, list]]]: """ Map from the petab to the amici format of noise distribution identifiers. @@ -756,7 +757,13 @@ def petab_noise_distributions_to_amici(observable_df: pd.DataFrame amici_val += observable[NOISE_DISTRIBUTION] else: amici_val += 'normal' - amici_distrs[observable.name] = amici_val + amici_distrs[observable.name] = {'type': amici_val, 'parameters': []} + + DISTRIBUTION_PARAMETERS = 'distributionParameters' # TODO + if DISTRIBUTION_PARAMETERS in observable \ + and str(observable[DISTRIBUTION_PARAMETERS]): + amici_distrs[observable.name]['parameters'] = \ + str(observable[DISTRIBUTION_PARAMETERS]).split(';') return amici_distrs diff --git a/python/sdist/amici/sbml_import.py b/python/sdist/amici/sbml_import.py index e1f35245e4..0ac169f7f9 100644 --- a/python/sdist/amici/sbml_import.py +++ b/python/sdist/amici/sbml_import.py @@ -47,6 +47,7 @@ class SBMLException(Exception): } ConservationLaw = Dict[str, Union[str, sp.Expr]] +NoiseDistDict = Dict[str, Union[str, list]] logger = get_logger(__name__, logging.ERROR) @@ -215,7 +216,8 @@ def sbml2amici( constant_parameters: Iterable[str] = None, sigmas: Dict[str, Union[str, float]] = None, event_sigmas: Dict[str, Union[str, float]] = None, - noise_distributions: Dict[str, Union[str, Callable]] = None, + noise_distributions: Optional[ + Dict[str, Union[NoiseDistDict, Callable]]] = None, event_noise_distributions: Dict[str, Union[str, Callable]] = None, verbose: Union[int, bool] = logging.ERROR, assume_pow_positivity: bool = False, @@ -374,7 +376,8 @@ def _build_ode_model( constant_parameters: Iterable[str] = None, sigmas: Dict[str, Union[str, float]] = None, event_sigmas: Dict[str, Union[str, float]] = None, - noise_distributions: Dict[str, Union[str, Callable]] = None, + noise_distributions: Optional[ + Dict[str, Union[NoiseDistDict, Callable]]] = None, event_noise_distributions: Dict[str, Union[str, Callable]] = None, verbose: Union[int, bool] = logging.ERROR, compute_conservation_laws: bool = True, @@ -1216,7 +1219,7 @@ def _process_observables( self, observables: Union[Dict[str, Dict[str, str]], None], sigmas: Dict[str, Union[str, float]], - noise_distributions: Dict[str, str] + noise_distributions: Dict[str, NoiseDistDict] ) -> None: """ Perform symbolic computations required for observable and objective @@ -1224,7 +1227,7 @@ def _process_observables( :param observables: dictionary(observableId: {'name':observableName - (optional), 'formula':formulaString)}) + (optional), 'formula':formulaString}) to be added to the model :param sigmas: @@ -1232,7 +1235,8 @@ def _process_observables( parameter name) :param noise_distributions: - dictionary(observableId: noise type) + dictionary(observableId: {'type': noise type, 'parameters': + list with parameter values}) See :py:func:`sbml2amici`. """ @@ -1254,7 +1258,8 @@ def _process_observables( ), 'transformation': noise_distribution_to_observable_transformation( - noise_distributions.get(obs, 'normal') + noise_distributions.get(obs, {'type': 'normal', + 'parameters': []}) ) } for iobs, (obs, definition) in enumerate(observables.items()) @@ -1325,7 +1330,7 @@ def _process_event_observables( 'event': sp.Symbol(definition.get('event')), 'transformation': noise_distribution_to_observable_transformation( - event_noise_distributions.get(obs, 'normal') + event_noise_distributions.get(obs, 'normal') # TODO ) } for iobs, (obs, definition) in @@ -1383,7 +1388,7 @@ def _generate_default_observables(self): def _process_log_likelihood(self, sigmas: Dict[str, Union[str, float]], - noise_distributions: Dict[str, str], + noise_distributions: Dict[str, NoiseDistDict], events: bool = False, event_reg: bool = False): """ @@ -1444,7 +1449,8 @@ def _process_log_likelihood(self, self.symbols[sigma_symbol].items() ): symbol = symbol_with_assumptions(f'J{obs_id}') - dist = noise_distributions.get(str(obs_id), 'normal') + dist = noise_distributions.get( + str(obs_id), {'type': 'normal', 'parameters': []}) cost_fun = noise_distribution_to_cost_function(dist)(obs_id) value = sp.sympify(cost_fun, locals=dict(zip( _get_str_symbol_identifiers(obs_id), @@ -1456,7 +1462,7 @@ def _process_log_likelihood(self, self.symbols[llh_symbol][symbol] = { 'name': f'J{obs["name"]}', 'value': value, - 'dist': dist, + 'dist': dist['type'], } @log_execution_time('processing SBML initial assignments', logger) @@ -2325,7 +2331,7 @@ def _parse_special_functions_sbml(sym: sp.Expr, def _validate_observables( observables: Union[Dict[str, Dict[str, str]], None], sigmas: Dict[str, Union[str, float]], - noise_distributions: Dict[str, str], + noise_distributions: Dict[str, NoiseDistDict], events: bool = False ) -> None: From e6db1007e88cb3fba3c7e669a1963179eb4d3852 Mon Sep 17 00:00:00 2001 From: Polina Lakrisenko Date: Thu, 12 Jan 2023 15:35:42 +0100 Subject: [PATCH 2/5] add log-left-truncated-normal --- python/sdist/amici/import_utils.py | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/python/sdist/amici/import_utils.py b/python/sdist/amici/import_utils.py index 17ff0a01f2..a4e6315a04 100644 --- a/python/sdist/amici/import_utils.py +++ b/python/sdist/amici/import_utils.py @@ -186,7 +186,7 @@ def noise_distribution_to_cost_function( if isinstance(noise_distribution, Callable): return noise_distribution - noise_distribution_type = noise_distribution['dist_type'] + noise_distribution_type = noise_distribution['type'] if noise_distribution_type in ['normal', 'lin-normal']: y_string = '0.5*log(2*pi*{sigma}**2) + 0.5*(({y} - {m}) / {sigma})**2' @@ -196,13 +196,20 @@ def noise_distribution_to_cost_function( elif noise_distribution_type == 'log10-normal': y_string = '0.5*log(2*pi*{sigma}**2*{m}**2*log(10)**2) ' \ '+ 0.5*((log({y}, 10) - log({m}, 10)) / {sigma})**2' - elif noise_distribution_type == ['left-truncated-normal', + elif noise_distribution_type in ['left-truncated-normal', 'lin-left-truncated-normal']: # truncated at a - a = noise_distribution['dist_params'][0] # TODO + a = noise_distribution['parameters'][0] # TODO y_string = f'log(1-0.5*(1+erf({a}-{{y}}/(sqrt(2)*{{sigma}}))))' \ ' + 0.5*log(2*pi*{sigma}**2)' \ ' + 0.5*(({y} - {m}) / {sigma})**2' + elif noise_distribution_type == 'log-left-truncated-normal': + # truncated at a + a = noise_distribution['parameters'][0] # TODO + # TODO a>0 + y_string = f'log(1-0.5*(1+erf(log({a})-log({{y}})/(sqrt(2)*{{sigma}}))))' \ + ' + 0.5*log(2*pi*{sigma}**2*{m}**2)' \ + ' + 0.5*((log({y}) - log({m})) / {sigma})**2' elif noise_distribution_type in ['laplace', 'lin-laplace']: y_string = 'log(2*{sigma}) + Abs({y} - {m}) / {sigma}' elif noise_distribution_type == 'log-laplace': @@ -227,7 +234,7 @@ def noise_distribution_to_cost_function( f'- {{m}} * log({{sigma}})' elif noise_distribution_type == 'left-censored-normal': # left-censored at v (detection limit) - v = noise_distribution['dist_params'][0] # TODO + v = noise_distribution['parameters'][0] # TODO y_string = f'log(0.5*(1+erf(({v}-{{y}})/(sqrt(2)*{{sigma}}))))' else: raise ValueError( From 641c851a337da74eeef85585f561f07714a10f27 Mon Sep 17 00:00:00 2001 From: Polina Lakrisenko Date: Mon, 23 Jan 2023 17:56:57 +0100 Subject: [PATCH 3/5] fix --- python/sdist/amici/import_utils.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/python/sdist/amici/import_utils.py b/python/sdist/amici/import_utils.py index a4e6315a04..72ce902d54 100644 --- a/python/sdist/amici/import_utils.py +++ b/python/sdist/amici/import_utils.py @@ -200,7 +200,7 @@ def noise_distribution_to_cost_function( 'lin-left-truncated-normal']: # truncated at a a = noise_distribution['parameters'][0] # TODO - y_string = f'log(1-0.5*(1+erf({a}-{{y}}/(sqrt(2)*{{sigma}}))))' \ + y_string = f'log(1-0.5*(1+erf(({a}-{{y}})/(sqrt(2)*{{sigma}}))))' \ ' + 0.5*log(2*pi*{sigma}**2)' \ ' + 0.5*(({y} - {m}) / {sigma})**2' elif noise_distribution_type == 'log-left-truncated-normal': @@ -232,7 +232,8 @@ def noise_distribution_to_cost_function( y_string = f'- loggamma({{m}}+{r}) + loggamma({{m}}+1) ' \ f'+ loggamma({r}) - {r} * log(1-{{sigma}}) ' \ f'- {{m}} * log({{sigma}})' - elif noise_distribution_type == 'left-censored-normal': + elif noise_distribution_type in ['lin-left-censored-normal', + 'left-censored-normal']: # left-censored at v (detection limit) v = noise_distribution['parameters'][0] # TODO y_string = f'log(0.5*(1+erf(({v}-{{y}})/(sqrt(2)*{{sigma}}))))' From c3506fb92685eeba39ddd39f503386e5e94f19b5 Mon Sep 17 00:00:00 2001 From: Polina Lakrisenko Date: Thu, 26 Jan 2023 13:20:28 +0100 Subject: [PATCH 4/5] add right-truncated and right-censored --- python/sdist/amici/import_utils.py | 84 ++++++++++++++++++++++++++++-- 1 file changed, 80 insertions(+), 4 deletions(-) diff --git a/python/sdist/amici/import_utils.py b/python/sdist/amici/import_utils.py index 72ce902d54..f075099b9b 100644 --- a/python/sdist/amici/import_utils.py +++ b/python/sdist/amici/import_utils.py @@ -105,9 +105,24 @@ def noise_distribution_to_cost_function( \\pi(m|y,\\sigma) = \\frac{1}{\\sqrt{2\\pi}\\sigma m \\log(10)}\\ exp\\left(-\\frac{(\\log_{10} m - \\log_{10} y)^2}{2\\sigma^2}\\right) - - `'truncated-normal'`: A truncated normal distribution: + - `'left-truncated-normal'`: A truncated normal distribution: .. math:: + \\pi(m|y,\\sigma) = \\frac{1}{1-\\Phi(\\frac{a-y}{\\sigma})} + \\frac{1}{\\sqrt{2\\pi}\\sigma}\\ + exp\\left(-\\frac{(m-y)^2}{2\\sigma^2}\\right) + + where \\Phi is the standard normal cumulative distribution function + + .. math:: + \\Phi(x) = \\frac{1}{2}(1+\\erf(\\frac{x}{\\sqrt{2}})) + + - `'right-truncated-normal'`: A truncated normal distribution: + + .. math:: + \\pi(m|y,\\sigma) = \\frac{1}{\\Phi(\\frac{b-y}{\\sigma})} + \\frac{1}{\\sqrt{2\\pi}\\sigma}\\ + exp\\left(-\\frac{(m-y)^2}{2\\sigma^2}\\right) - `'laplace'`, `'lin-laplace'`: A laplace distribution: @@ -152,6 +167,44 @@ def noise_distribution_to_cost_function( .. math:: r = \\frac{1-\\sigma}{\\sigma} y + - `'left-censored-normal'`: A left-censored normal distribution with + threshold v: + + - if m <= v: + + .. math:: + \\Phi(\\frac{v-y}{\\sigma}) + + where \\Phi is the standard normal cumulative distribution function + + .. math:: + \\Phi(x) = \\frac{1}{2}(1+\\erf(\\frac{x}{\\sqrt{2}})) + + - if m > v: + + .. math:: + \\pi(m|y,\\sigma) = \\frac{1}{\\sqrt{2\\pi}\\sigma}\\ + exp\\left(-\\frac{(m-y)^2}{2\\sigma^2}\\right) + + - `'right-censored-normal'`: A right-censored normal distribution with + threshold v: + + - if m < v: + + .. math:: + \\pi(m|y,\\sigma) = \\frac{1}{\\sqrt{2\\pi}\\sigma}\\ + exp\\left(-\\frac{(m-y)^2}{2\\sigma^2}\\right) + + - if m >= v: + + .. math:: + 1-\\Phi(\\frac{v-y}{\\sigma}) + + where \\Phi is the standard normal cumulative distribution function + + .. math:: + \\Phi(x) = \\frac{1}{2}(1+\\erf(\\frac{x}{\\sqrt{2}})) + The distributions above are for a single data point. For a collection :math:`D=\\{m_i\\}_i` of data points and corresponding simulations :math:`Y=\\{y_i\\}_i` and noise parameters @@ -198,16 +251,24 @@ def noise_distribution_to_cost_function( '+ 0.5*((log({y}, 10) - log({m}, 10)) / {sigma})**2' elif noise_distribution_type in ['left-truncated-normal', 'lin-left-truncated-normal']: - # truncated at a + # left-truncated at a a = noise_distribution['parameters'][0] # TODO y_string = f'log(1-0.5*(1+erf(({a}-{{y}})/(sqrt(2)*{{sigma}}))))' \ ' + 0.5*log(2*pi*{sigma}**2)' \ ' + 0.5*(({y} - {m}) / {sigma})**2' + elif noise_distribution_type in ['right-truncated-normal', + 'lin-right-truncated-normal']: + # right-truncated at b + b = noise_distribution['parameters'][0] # TODO + y_string = f'log(0.5*(1+erf(({b}-{{y}})/(sqrt(2)*{{sigma}}))))' \ + ' + 0.5*log(2*pi*{sigma}**2)' \ + ' + 0.5*(({y} - {m}) / {sigma})**2' elif noise_distribution_type == 'log-left-truncated-normal': # truncated at a a = noise_distribution['parameters'][0] # TODO # TODO a>0 - y_string = f'log(1-0.5*(1+erf(log({a})-log({{y}})/(sqrt(2)*{{sigma}}))))' \ + y_string = f'log(1-0.5*(1+erf(log({a})-log({{y}})/' \ + f'(sqrt(2)*{{sigma}}))))' \ ' + 0.5*log(2*pi*{sigma}**2*{m}**2)' \ ' + 0.5*((log({y}) - log({m})) / {sigma})**2' elif noise_distribution_type in ['laplace', 'lin-laplace']: @@ -236,7 +297,22 @@ def noise_distribution_to_cost_function( 'left-censored-normal']: # left-censored at v (detection limit) v = noise_distribution['parameters'][0] # TODO - y_string = f'log(0.5*(1+erf(({v}-{{y}})/(sqrt(2)*{{sigma}}))))' + y_string = f'Piecewise((-log(0.5*(1+erf(({v}-{{y}})/' \ + f'(sqrt(2)*{{sigma}})))), ' \ + f'{{m}}<={v}), ' \ + f'(0.5*log(2*pi*{{sigma}}**2) + ' \ + f'0.5*(({{y}} - {{m}}) / {{sigma}})**2, ' \ + f'{{m}}>{v}))' + elif noise_distribution_type in ['lin-right-censored-normal', + 'right-censored-normal']: + # right-censored at v (detection limit) + v = noise_distribution['parameters'][0] # TODO + y_string = f'Piecewise((-log(1-0.5*(1+erf(({v}-{{y}})/' \ + f'(sqrt(2)*{{sigma}})))), ' \ + f'{{m}}>={v}), ' \ + f'(0.5*log(2*pi*{{sigma}}**2) + ' \ + f'0.5*(({{y}} - {{m}}) / {{sigma}})**2, ' \ + f'{{m}}<{v}))' else: raise ValueError( f"Cost identifier {noise_distribution_type} not recognized.") From bbd9cc6a21177fa35fb34628b5b60043ce15f689 Mon Sep 17 00:00:00 2001 From: Polina Lakrisenko Date: Mon, 20 Feb 2023 14:28:17 +0100 Subject: [PATCH 5/5] add censored laplace --- python/sdist/amici/import_utils.py | 60 +++++++++++++++++++++++++++++- 1 file changed, 58 insertions(+), 2 deletions(-) diff --git a/python/sdist/amici/import_utils.py b/python/sdist/amici/import_utils.py index f075099b9b..3b12a5b6db 100644 --- a/python/sdist/amici/import_utils.py +++ b/python/sdist/amici/import_utils.py @@ -173,7 +173,7 @@ def noise_distribution_to_cost_function( - if m <= v: .. math:: - \\Phi(\\frac{v-y}{\\sigma}) + \\pi(m|y,\\sigma) = \\Phi(\\frac{v-y}{\\sigma}) where \\Phi is the standard normal cumulative distribution function @@ -198,13 +198,51 @@ def noise_distribution_to_cost_function( - if m >= v: .. math:: - 1-\\Phi(\\frac{v-y}{\\sigma}) + \\pi(m|y,\\sigma) = 1-\\Phi(\\frac{v-y}{\\sigma}) where \\Phi is the standard normal cumulative distribution function .. math:: \\Phi(x) = \\frac{1}{2}(1+\\erf(\\frac{x}{\\sqrt{2}})) + - `'left-censored-laplace'`: A left-censored laplace distribution with + threshold v: + + - if m <= v: + + .. math:: + \\pi(m|y,\\sigma) = \\F(\\frac{v-y}{\\sigma}) + + where \\Phi is the standard laplace cumulative distribution function + + .. math:: + \\F(x) = \\frac{1}{2}(1+\\sign(x)(1-\\exp(-|x|))) + + - if m > v: + + .. math:: + \\pi(m|y,\\sigma) = \\frac{1}{2\\sigma} + \\exp\\left(-\\frac{|m-y|}{\\sigma}\\right) + + - `'right-censored-laplace'`: A right-censored laplace distribution with + threshold v: + + - if m < v: + + .. math:: + \\pi(m|y,\\sigma) = \\frac{1}{2\\sigma} + \\exp\\left(-\\frac{|m-y|}{\\sigma}\\right) + + - if m >= v: + + .. math:: + \\pi(m|y,\\sigma) = 1-\\Phi(\\frac{v-y}{\\sigma}) + + where \\Phi is the standard laplace cumulative distribution function + + .. math:: + \\Phi(x) = \\frac{1}{2}(1+\\sign(x)(1-\\exp(-|x|))) + The distributions above are for a single data point. For a collection :math:`D=\\{m_i\\}_i` of data points and corresponding simulations :math:`Y=\\{y_i\\}_i` and noise parameters @@ -313,6 +351,24 @@ def noise_distribution_to_cost_function( f'(0.5*log(2*pi*{{sigma}}**2) + ' \ f'0.5*(({{y}} - {{m}}) / {{sigma}})**2, ' \ f'{{m}}<{v}))' + elif noise_distribution_type in ['lin-left-censored-laplace', + 'left-censored-laplace']: + # left-censored at v (detection limit) + v = noise_distribution['parameters'][0] # TODO + y_string = f'Piecewise((-log(0.5*(1+sign({v}-{{y}})*' \ + f'(1-exp(-Abs({v}-{{y}})/{{sigma}}))), ' \ + f'{{m}}<={v}), ' \ + f'(log(2*{{sigma}}) + Abs({{y}} - {{m}}) / {{sigma}}, ' \ + f'{{m}}>{v}))' + elif noise_distribution_type in ['lin-right-censored-laplace', + 'right-censored-laplace']: + # right-censored at v (detection limit) + v = noise_distribution['parameters'][0] # TODO + y_string = f'Piecewise((-log(1-0.5*(1+sign({v}-{{y}})*' \ + f'(1-exp(-Abs({v}-{{y}})/{{sigma}}))), ' \ + f'{{m}}>={v}), ' \ + f'(log(2*{{sigma}}) + Abs({{y}} - {{m}}) / {{sigma}}, ' \ + f'{{m}}<{v}))' else: raise ValueError( f"Cost identifier {noise_distribution_type} not recognized.")