diff --git a/csep/core/binomial_evaluations.py b/csep/core/binomial_evaluations.py index 9ecd61ea..ec187cf4 100644 --- a/csep/core/binomial_evaluations.py +++ b/csep/core/binomial_evaluations.py @@ -88,13 +88,13 @@ def binary_joint_log_likelihood_ndarray(forecast, catalog): It has to be a either zero or positive integer only (No Floating Point) """ #First, we mask the forecast in cells where we could find log=0.0 singularities: - forecast_masked = np.ma.masked_where(forecast.ravel() <= 0.0, forecast.ravel()) + forecast_masked = numpy.ma.masked_where(forecast.ravel() <= 0.0, forecast.ravel()) #Then, we compute the log-likelihood of observing one or more events given a Poisson distribution, i.e., 1 - Pr(0) target_idx = numpy.nonzero(catalog.ravel()) y = numpy.zeros(forecast_masked.ravel().shape) y[target_idx[0]] = 1 - first_term = y * (np.log(1.0 - np.exp(-forecast_masked.ravel()))) + first_term = y * (numpy.log(1.0 - numpy.exp(-forecast_masked.ravel()))) #Also, we estimate the log-likelihood in cells no events are observed: second_term = (1-y) * (-forecast_masked.ravel().data) diff --git a/tests/test_evaluations.py b/tests/test_evaluations.py index df6e22e6..ffc6f3ce 100644 --- a/tests/test_evaluations.py +++ b/tests/test_evaluations.py @@ -3,6 +3,8 @@ import unittest from csep.core.poisson_evaluations import _simulate_catalog, _poisson_likelihood_test +from csep.core.binomial_evaluations import binary_joint_log_likelihood_ndarray + def get_datadir(): root_dir = os.path.dirname(os.path.abspath(__file__)) @@ -71,3 +73,17 @@ def test_likelihood(self): # calculated by hand given the expected data, see explanation in zechar et al., 2010. numpy.testing.assert_allclose(simulated_ll[0], -7.178053830347945) + +class TestBinomialLikelihood(unittest.TestCase): + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + self.forecast_data = numpy.array([[0.1, 0.3, 0.4], [0.2, 0.1, 0.1]]) + self.observed_data = numpy.array([[0, 1, 2], [1, 1, 0]]) + + def test_likelihood(self): + bill = binary_joint_log_likelihood_ndarray(self.forecast_data, self.observed_data) + + numpy.testing.assert_allclose(bill, -6.7197988064) + +if __name__ == '__main__': + unittest.main() \ No newline at end of file