Skip to content

Commit

Permalink
correct empricial_cdf for case with weights and include tests
Browse files Browse the repository at this point in the history
  • Loading branch information
aronsho committed Jul 3, 2024
1 parent 3639b21 commit 7b19cc5
Show file tree
Hide file tree
Showing 2 changed files with 22 additions and 8 deletions.
15 changes: 7 additions & 8 deletions seismostats/analysis/estimate_mc.py
Original file line number Diff line number Diff line change
Expand Up @@ -95,8 +95,8 @@ def empirical_cdf(
if mc is None:
mc = np.min(sample)

idx = np.argsort(sample)
x = sample[idx]
idx1 = np.argsort(sample)
x = sample[idx1]
x, y_count = np.unique(x, return_counts=True)

# add empty bins
Expand All @@ -106,24 +106,23 @@ def empirical_cdf(
if mag_bin not in x:
x = np.append(x, mag_bin)
y_count = np.append(y_count, 0)
idx = np.argsort(x)
x = x[idx]
y_count = y_count[idx]
idx2 = np.argsort(x)
x = x[idx2]
y_count = y_count[idx2]

# estimate the CDF
if weights is None:
y = np.cumsum(y_count) / len(sample)
else:
weights_sorted = weights[idx]
weights_sorted = weights[idx1]
y = np.cumsum(weights_sorted) / weights_sorted.sum()

# make sure that y is zero if there are no samples in the first bins
for ii, y_loop in enumerate(y_count):
if y_loop > 0:
break
leading_zeros = np.zeros(ii)

y = leading_zeros.append(y[np.cumsum(y_count[ii:]) - 1])
y = np.append(leading_zeros, y[np.cumsum(y_count[ii:]) - 1])

return x, y

Expand Down
15 changes: 15 additions & 0 deletions seismostats/analysis/tests/test_estimate_mc.py
Original file line number Diff line number Diff line change
Expand Up @@ -128,6 +128,21 @@ def test_empirical_cdf(setup_magnitudes, delta_m=0.1):
assert_equal(len(x), len(y))
assert_equal(y[0], 0.06)

# test that weights function the way that they should
# 1. with equal weights
magnitudes = np.array([0.5, 0.6, 0.7, 0.8, 0.9])
weights = np.array([1, 1, 1, 1, 1])
x, y = empirical_cdf(magnitudes, mc=0, delta_m=0.1, weights=weights)
assert_almost_equal(x, np.arange(0, 1, 0.1))
assert_almost_equal(y, [0, 0, 0, 0, 0, 0.2, 0.4, 0.6, 0.8, 1])

# 2. with different weights
magnitudes = np.array([0.5, 0.6, 0.7, 0.8, 0.9])
weights = np.array([0.5, 0.5, 0.5, 0.5, 3])
x, y = empirical_cdf(magnitudes, mc=0, delta_m=0.1, weights=weights)
assert_almost_equal(x, np.arange(0, 1, 0.1))
assert_almost_equal(y, [0, 0, 0, 0, 0, 0.1, 0.2, 0.3, 0.4, 1])


@pytest.fixture
def setup_ks_dists():
Expand Down

0 comments on commit 7b19cc5

Please sign in to comment.