diff --git a/cara/monte_carlo/data.py b/cara/monte_carlo/data.py index 2c6268e6..b5e15884 100644 --- a/cara/monte_carlo/data.py +++ b/cara/monte_carlo/data.py @@ -5,7 +5,7 @@ import numpy as np from scipy import special as sp import cara.monte_carlo as mc -from cara.monte_carlo.sampleable import LogNormal,LogCustomKernel,CustomKernel,Uniform +from cara.monte_carlo.sampleable import LogNormal,LogCustomKernel,CustomKernel,Uniform, Beta sqrt2pi = np.sqrt(2.*np.pi) @@ -202,5 +202,5 @@ short_range_expiration_distributions = { } -# Fit from Fig 8 a) "stand-stand" in https://www.mdpi.com/1660-4601/17/4/1445/htm -short_range_distances = LogNormal(-0.269359136417347, 0.4728300188814934) +# Derived from Fig 8 a) "stand-stand" in https://www.mdpi.com/1660-4601/17/4/1445/htm +short_range_distances = Beta(alpha=1.4342766632654418, beta=27.49916410927064) \ No newline at end of file diff --git a/cara/monte_carlo/sampleable.py b/cara/monte_carlo/sampleable.py index 27907e49..5a29de81 100644 --- a/cara/monte_carlo/sampleable.py +++ b/cara/monte_carlo/sampleable.py @@ -2,6 +2,7 @@ import typing import numpy as np from sklearn.neighbors import KernelDensity # type: ignore +from scipy.stats import beta import cara.models @@ -130,6 +131,22 @@ class LogCustomKernel(SampleableDistribution): return 10 ** kde_model.sample(n_samples=size)[:, 0] +class Beta(SampleableDistribution): + """ + Defines a Beta distribution parameterized by two positive shape parameters, + denoted by alpha (α) and beta (β), that appear as exponents of the random + variable and control the shape of the distribution. + """ + + def __init__(self, alpha: float, beta: float): + # these are resp. the alpha and beta of the underlying distribution + self.alpha = alpha + self.beta = beta + + def generate_samples(self, size: int) -> float_array_size_n: + return beta.rvs(a = self.alpha, b = self.beta, loc=0.5, scale=7, size=size) + + _VectorisedFloatOrSampleable = typing.Union[ SampleableDistribution, cara.models._VectorisedFloat, ] diff --git a/cara/tests/models/test_short_range_model.py b/cara/tests/models/test_short_range_model.py index a07a656b..4a191334 100644 --- a/cara/tests/models/test_short_range_model.py +++ b/cara/tests/models/test_short_range_model.py @@ -72,9 +72,9 @@ def test_dilution_factor(activity, expected_dilution): @pytest.mark.parametrize( "time, expected_short_range_concentration", [ [8.5, 0.], - [10.5, 15.24806213], - [10.6, 15.24806213], - [11.0, 15.24806213], + [10.5, 7.491771], + [10.6, 7.491771], + [11.0, 7.491771], [12.0, 0.], ] )