diff --git a/src/pysatl_core/families/builtins/continuous/exponential.py b/src/pysatl_core/families/builtins/continuous/exponential.py index e0854ad..4794ce3 100644 --- a/src/pysatl_core/families/builtins/continuous/exponential.py +++ b/src/pysatl_core/families/builtins/continuous/exponential.py @@ -165,6 +165,55 @@ def char_func(parameters: Parametrization, t: NumericArray) -> ComplexArray: ) return cast(ComplexArray, result) + def lpdf(parameters: Parametrization, x: NumericArray) -> NumericArray: + """ + Logarithm of the probability density function for exponential distribution. + + Parameters + ---------- + parameters : Parametrization + Distribution parameters object with field: + - lambda_: float (rate parameter) + x : NumericArray + Points at which to evaluate the log-probability density function + + Returns + ------- + NumericArray + Log-probability density values at points x. + For x < 0 returns -np.inf. + """ + parameters = cast(_Rate, parameters) + lambda_ = parameters.lambda_ + return np.where(x >= 0, np.log(lambda_) - lambda_ * x, -np.inf) + + def score(parameters: Parametrization, x: NumericArray) -> np.ndarray: + """ + Score function (gradient of log-pdf w.r.t. lambda) for exponential distribution. + + Parameters + ---------- + parameters : Parametrization + Distribution parameters object with field: + - lambda_: float (rate parameter) + x : NumericArray + Points at which to evaluate the score. Can be a scalar or 1D array. + + Returns + ------- + np.ndarray + Array of shape (..., 1). For scalar input, returns a 1D array of length 1. + For 1D array input of length N, returns a 2D array of shape (N, 1), + where each row corresponds to a data point and column corresponds to + gradient w.r.t. lambda. + """ + parameters = cast(_Rate, parameters) + lambda_ = parameters.lambda_ + # derivative of log pdf w.r.t lambda: 1/lambda - x (for x >= 0) + grad = np.where(x >= 0, 1.0 / lambda_ - x, 0.0) + # add trailing dimension to get (..., 1) + return grad[..., np.newaxis] + def mean_func(parameters: Parametrization, _: Any) -> float: """Mean of exponential distribution.""" parameters = cast(_Rate, parameters) @@ -213,6 +262,8 @@ def _support(_: Parametrization) -> ContinuousSupport: CharacteristicName.CDF: cdf, CharacteristicName.PPF: ppf, CharacteristicName.CF: char_func, + CharacteristicName.LPDF: lpdf, + CharacteristicName.SCORE: score, CharacteristicName.MEAN: mean_func, CharacteristicName.VAR: var_func, CharacteristicName.SKEW: skew_func, diff --git a/src/pysatl_core/families/builtins/continuous/normal.py b/src/pysatl_core/families/builtins/continuous/normal.py index bf86be1..2634852 100644 --- a/src/pysatl_core/families/builtins/continuous/normal.py +++ b/src/pysatl_core/families/builtins/continuous/normal.py @@ -168,6 +168,66 @@ def char_func(parameters: Parametrization, t: NumericArray) -> ComplexArray: mu = parameters.mu return cast(ComplexArray, np.exp(1j * mu * t - 0.5 * (sigma**2) * (t**2))) + def lpdf(parameters: Parametrization, x: NumericArray) -> NumericArray: + """ + Logarithm of the probability density function for normal distribution. + + Parameters + ---------- + parameters : Parametrization + Distribution parameters object with fields: + - mu: float (mean) + - sigma: float (standard deviation) + x : NumericArray + Points at which to evaluate the log-probability density function + + Returns + ------- + NumericArray + Log-probability density values at points x + """ + parameters = cast(_MeanStd, parameters) + sigma = parameters.sigma + mu = parameters.mu + + # log pdf = -0.5*log(2π) - log(σ) - (x-μ)²/(2σ²) + return cast( + NumericArray, + -0.5 * np.log(2 * np.pi) - np.log(sigma) - ((x - mu) ** 2) / (2 * sigma**2), + ) + + def score(parameters: Parametrization, x: NumericArray) -> np.ndarray: + """ + Score function (gradient of log-pdf) for normal distribution. + + Parameters + ---------- + parameters : Parametrization + Distribution parameters object with fields: + - mu: float (mean) + - sigma: float (standard deviation) + x : NumericArray + Points at which to evaluate the score. Can be a scalar or 1D array. + + Returns + ------- + np.ndarray + Array of shape (..., 2). For scalar input, returns a 1D array of length 2. + For 1D array input of length N, returns a 2D array of shape (N, 2), + where each row corresponds to a data point and columns correspond to + gradients w.r.t. mu and sigma (in that order). + """ + parameters = cast(_MeanStd, parameters) + mu = parameters.mu + sigma = parameters.sigma + + x_minus_mu = x - mu + grad_mu = x_minus_mu / sigma**2 + grad_sigma = -1.0 / sigma + (x_minus_mu**2) / sigma**3 + + # Stack along last axis to get (..., 2) + return np.stack([grad_mu, grad_sigma], axis=-1) + def mean_func(parameters: Parametrization, _: Any) -> float: """Mean of normal distribution.""" parameters = cast(_MeanStd, parameters) @@ -216,6 +276,8 @@ def _support(_: Parametrization) -> ContinuousSupport: CharacteristicName.CDF: cdf, CharacteristicName.PPF: ppf, CharacteristicName.CF: char_func, + CharacteristicName.LPDF: lpdf, + CharacteristicName.SCORE: score, CharacteristicName.MEAN: mean_func, CharacteristicName.VAR: var_func, CharacteristicName.SKEW: skew_func, diff --git a/src/pysatl_core/families/builtins/continuous/uniform.py b/src/pysatl_core/families/builtins/continuous/uniform.py index 6d73503..f854272 100644 --- a/src/pysatl_core/families/builtins/continuous/uniform.py +++ b/src/pysatl_core/families/builtins/continuous/uniform.py @@ -191,6 +191,65 @@ def char_func(parameters: Parametrization, t: NumericArray) -> ComplexArray: return cast(ComplexArray, sinc_val * np.exp(1j * center * t_arr)) + def lpdf(parameters: Parametrization, x: NumericArray) -> NumericArray: + """ + Logarithm of the probability density function for uniform distribution. + + Parameters + ---------- + parameters : Parametrization + Distribution parameters object with fields: + - lower_bound: float (lower bound) + - upper_bound: float (upper bound) + x : NumericArray + Points at which to evaluate the log-probability density function + + Returns + ------- + NumericArray + Log-probability density values at points x + For x outside [lower_bound, upper_bound] returns -np.inf + """ + parameters = cast(_Standard, parameters) + a = parameters.lower_bound + b = parameters.upper_bound + width = b - a + return np.where((x >= a) & (x <= b), -np.log(width), -np.inf) + + def score(parameters: Parametrization, x: NumericArray) -> np.ndarray: + """ + Score function (gradient of log-pdf) for uniform distribution. + + Parameters + ---------- + parameters : Parametrization + Distribution parameters object with fields: + - lower_bound: float (lower bound) + - upper_bound: float (upper bound) + x : NumericArray + Points at which to evaluate the score. Can be a scalar or 1D array. + + Returns + ------- + np.ndarray + Array of shape (..., 2). For scalar input, returns a 1D array of length 2. + For 1D array input of length N, returns a 2D array of shape (N, 2), + where each row corresponds to a data point and columns correspond to + gradients w.r.t. lower_bound and upper_bound (in that order). + """ + parameters = cast(_Standard, parameters) + a = parameters.lower_bound + b = parameters.upper_bound + width = b - a + inside = (x >= a) & (x <= b) + + # Derivatives: + # d/da log pdf = 1/width, d/db log pdf = -1/width for points inside support + grad_a = np.where(inside, 1.0 / width, 0.0) + grad_b = np.where(inside, -1.0 / width, 0.0) + + return np.stack([grad_a, grad_b], axis=-1) + def mean_func(parameters: Parametrization, _: Any) -> float: """Mean of uniform distribution.""" parameters = cast(_Standard, parameters) @@ -248,6 +307,8 @@ def _support(parameters: Parametrization) -> ContinuousSupport: CharacteristicName.CDF: cdf, CharacteristicName.PPF: ppf, CharacteristicName.CF: char_func, + CharacteristicName.LPDF: lpdf, + CharacteristicName.SCORE: score, CharacteristicName.MEAN: mean_func, CharacteristicName.VAR: var_func, CharacteristicName.SKEW: skew_func, diff --git a/src/pysatl_core/types.py b/src/pysatl_core/types.py index c63702d..f36f869 100644 --- a/src/pysatl_core/types.py +++ b/src/pysatl_core/types.py @@ -261,6 +261,8 @@ class CharacteristicName(StrEnum): CDF = "cdf" PPF = "ppf" PMF = "pmf" + LPDF = "lpdf" # unimplemented in graph yet + SCORE = "score" # unimplemented in graph yet CF = "cf" # unimplemented in graph yet SF = "sf" # unimplemented in graph yet MEAN = "mean" # unimplemented in graph yet diff --git a/tests/unit/families/builtins/continuous/test_exponential.py b/tests/unit/families/builtins/continuous/test_exponential.py index 6a6566b..cd28b35 100644 --- a/tests/unit/families/builtins/continuous/test_exponential.py +++ b/tests/unit/families/builtins/continuous/test_exponential.py @@ -121,6 +121,8 @@ def test_analytical_computations_availability(self): CharacteristicName.CDF, CharacteristicName.PPF, CharacteristicName.CF, + CharacteristicName.LPDF, + CharacteristicName.SCORE, CharacteristicName.MEAN, CharacteristicName.VAR, CharacteristicName.SKEW, @@ -131,14 +133,30 @@ def test_analytical_computations_availability(self): @pytest.mark.parametrize( "char_name, test_data, scipy_func, scipy_kwargs", [ - (CharacteristicName.PDF, [-1.0, 0.0, 1.0, 2.0, 3.0, 4.0], expon.pdf, {"scale": 2.0}), - (CharacteristicName.CDF, [-1.0, 0.0, 1.0, 2.0, 3.0, 4.0], expon.cdf, {"scale": 2.0}), + ( + CharacteristicName.PDF, + [-1.0, 0.0, 1.0, 2.0, 3.0, 4.0], + expon.pdf, + {"scale": 2.0}, + ), + ( + CharacteristicName.CDF, + [-1.0, 0.0, 1.0, 2.0, 3.0, 4.0], + expon.cdf, + {"scale": 2.0}, + ), ( CharacteristicName.PPF, [0.001, 0.01, 0.1, 0.25, 0.5, 0.75, 0.9, 0.99, 0.999], expon.ppf, {"scale": 2.0}, ), + ( + CharacteristicName.LPDF, + [-1.0, 0.0, 1.0, 2.0, 3.0, 4.0], + expon.logpdf, + {"scale": 2.0}, + ), ], ) def test_array_input_for_characteristics(self, char_name, test_data, scipy_func, scipy_kwargs): @@ -155,6 +173,33 @@ def test_array_input_for_characteristics(self, char_name, test_data, scipy_func, self.assert_arrays_almost_equal(result_array, expected_array) + def test_score_function(self): + """Test score function (gradient of log-pdf) for exponential distribution.""" + dist = self.exponential_dist_example + score_func = dist.query_method(CharacteristicName.SCORE) + + lambda_ = 0.5 + + x_scalar = 2.0 + grad_scalar = score_func(x_scalar) + assert grad_scalar.shape == (1,) + + expected_grad = 1.0 / lambda_ - x_scalar # для x>=0 + assert abs(grad_scalar[0] - expected_grad) < self.CALCULATION_PRECISION + + x_out = -1.0 + grad_out = score_func(x_out) + assert grad_out.shape == (1,) + assert grad_out[0] == 0.0 + + x_array = np.array([-1.0, 0.0, 1.0, 2.0, 3.0, 4.0]) + grad_array = score_func(x_array) + assert grad_array.shape == (len(x_array), 1) + + for idx, x_val in enumerate(x_array): + expected = 1.0 / lambda_ - x_val if x_val >= 0 else 0.0 + assert abs(grad_array[idx, 0] - expected) < self.CALCULATION_PRECISION + def test_characteristic_function_array_input(self): """Test characteristic function calculation with array input.""" char_func = self.exponential_dist_example.query_method(CharacteristicName.CF) diff --git a/tests/unit/families/builtins/continuous/test_normal.py b/tests/unit/families/builtins/continuous/test_normal.py index 2ba8151..64babbe 100644 --- a/tests/unit/families/builtins/continuous/test_normal.py +++ b/tests/unit/families/builtins/continuous/test_normal.py @@ -149,6 +149,8 @@ def test_analytical_computations_availability(self): CharacteristicName.CDF, CharacteristicName.PPF, CharacteristicName.CF, + CharacteristicName.LPDF, + CharacteristicName.SCORE, CharacteristicName.MEAN, CharacteristicName.VAR, CharacteristicName.SKEW, @@ -177,6 +179,12 @@ def test_analytical_computations_availability(self): norm.ppf, {"loc": 2.0, "scale": 1.5}, ), + ( + CharacteristicName.LPDF, + [-1.0, 0.0, 1.0, 2.0, 3.0, 4.0], + norm.logpdf, + {"loc": 2.0, "scale": 1.5}, + ), ], ) def test_array_input_for_characteristics(self, char_name, test_data, scipy_func, scipy_kwargs): @@ -193,6 +201,32 @@ def test_array_input_for_characteristics(self, char_name, test_data, scipy_func, self.assert_arrays_almost_equal(result_array, expected_array) + def test_score_function(self): + """Test score function (gradient of log-pdf) for normal distribution.""" + dist = self.normal_dist_example + score_func = dist.query_method(CharacteristicName.SCORE) + + mu, sigma = 2.0, 1.5 + + x_scalar = 2.5 + grad_scalar = score_func(x_scalar) + assert grad_scalar.shape == (2,) + + expected_grad_mu = (x_scalar - mu) / sigma**2 + expected_grad_sigma = -1.0 / sigma + (x_scalar - mu) ** 2 / sigma**3 + assert abs(grad_scalar[0] - expected_grad_mu) < self.CALCULATION_PRECISION + assert abs(grad_scalar[1] - expected_grad_sigma) < self.CALCULATION_PRECISION + + x_array = np.array([-1.0, 0.0, 1.0, 2.0, 3.0, 4.0]) + grad_array = score_func(x_array) + assert grad_array.shape == (len(x_array), 2) + + for idx, x_val in enumerate(x_array): + expected_mu = (x_val - mu) / sigma**2 + expected_sigma = -1.0 / sigma + (x_val - mu) ** 2 / sigma**3 + assert abs(grad_array[idx, 0] - expected_mu) < self.CALCULATION_PRECISION + assert abs(grad_array[idx, 1] - expected_sigma) < self.CALCULATION_PRECISION + def test_characteristic_function_array_input(self): """Test characteristic function calculation with array input.""" char_func = self.normal_dist_example.query_method(CharacteristicName.CF) diff --git a/tests/unit/families/builtins/continuous/test_uniform.py b/tests/unit/families/builtins/continuous/test_uniform.py index 953c35e..6bcdc27 100644 --- a/tests/unit/families/builtins/continuous/test_uniform.py +++ b/tests/unit/families/builtins/continuous/test_uniform.py @@ -148,6 +148,8 @@ def test_analytical_computations_availability(self): CharacteristicName.CDF, CharacteristicName.PPF, CharacteristicName.CF, + CharacteristicName.LPDF, + CharacteristicName.SCORE, CharacteristicName.MEAN, CharacteristicName.VAR, CharacteristicName.SKEW, @@ -176,6 +178,12 @@ def test_analytical_computations_availability(self): uniform.ppf, {"loc": 2.0, "scale": 3.0}, ), + ( + CharacteristicName.LPDF, + [1.0, 2.0, 2.5, 3.0, 4.0, 5.0, 6.0], + uniform.logpdf, + {"loc": 2.0, "scale": 3.0}, + ), ], ) def test_array_input_for_characteristics(self, char_name, test_data, scipy_func, scipy_kwargs): @@ -192,6 +200,49 @@ def test_array_input_for_characteristics(self, char_name, test_data, scipy_func, self.assert_arrays_almost_equal(result_array, expected_array) + def test_score_function(self): + """Test score function (gradient of log-pdf) for uniform distribution.""" + dist = self.uniform_dist_example + score_func = dist.query_method(CharacteristicName.SCORE) + + a, b = 2.0, 5.0 + width = b - a + + x_scalar = 3.0 # In the interval + grad_scalar = score_func(x_scalar) + assert grad_scalar.shape == (2,) + + expected_grad_a = 1.0 / width + expected_grad_b = -1.0 / width + assert abs(grad_scalar[0] - expected_grad_a) < self.CALCULATION_PRECISION + assert abs(grad_scalar[1] - expected_grad_b) < self.CALCULATION_PRECISION + + x_out = 1.0 + grad_out = score_func(x_out) + assert grad_out.shape == (2,) + assert grad_out[0] == 0.0 + assert grad_out[1] == 0.0 + + x_boundary = a + grad_boundary = score_func(x_boundary) + assert grad_boundary[0] == expected_grad_a + assert grad_boundary[1] == expected_grad_b + + x_array = np.array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0]) + grad_array = score_func(x_array) + assert grad_array.shape == (len(x_array), 2) + + for idx, x_val in enumerate(x_array): + inside = (x_val >= a) & (x_val <= b) + if inside: + expected_a = 1.0 / width + expected_b = -1.0 / width + else: + expected_a = 0.0 + expected_b = 0.0 + assert abs(grad_array[idx, 0] - expected_a) < self.CALCULATION_PRECISION + assert abs(grad_array[idx, 1] - expected_b) < self.CALCULATION_PRECISION + def test_characteristic_function_array_input(self): """Test characteristic function calculation with array input.""" char_func = self.uniform_dist_example.query_method(CharacteristicName.CF)