Skip to content

Commit

Permalink
Remove possibility of 0 as a bound
Browse files Browse the repository at this point in the history
  • Loading branch information
zmbc committed Nov 29, 2023
1 parent 0304872 commit 3a3b80c
Show file tree
Hide file tree
Showing 2 changed files with 18 additions and 19 deletions.
33 changes: 16 additions & 17 deletions integration_tests/conftest.py
Original file line number Diff line number Diff line change
Expand Up @@ -201,35 +201,34 @@ def _calculate_bayes_factor(

@cache
def _fit_beta_distribution_to_ui(self, lb: float, ub: float) -> Tuple[float, float]:
assert lb > 0 or ub < 1
assert lb > 0 and ub < 1
# Inspired by https://stats.stackexchange.com/a/112671/
def objective(x):
# np.exp ensures they are always positive
a, b = np.exp(x)
dist = scipy.stats.beta(a=a, b=b)

squared_error_lower = (
self._quantile_squared_error(dist, lb, 0.025) if lb > 0 else 0
)
squared_error_upper = (
self._quantile_squared_error(dist, ub, 0.975) if ub < 1 else 0
)
squared_error_lower = self._quantile_squared_error(dist, lb, 0.025)
squared_error_upper = self._quantile_squared_error(dist, ub, 0.975)

return squared_error_lower + squared_error_upper

# It is quite important to start with a reasonable guess.
ui_midpoint = (lb + ub) / 2
# TODO: Is this reasonable!?
first_guess_concentration = 10
optimization_result = scipy.optimize.minimize(
objective,
x0=[
np.log(ui_midpoint * first_guess_concentration),
np.log((1 - ui_midpoint) * first_guess_concentration),
],
)
# Sometimes it warns that it may not have found a good solution,
# but the solution is very accurate.
for first_guess_concentration in [1_000, 100, 10]:
optimization_result = scipy.optimize.minimize(
objective,
x0=[
np.log(ui_midpoint * first_guess_concentration),
np.log((1 - ui_midpoint) * first_guess_concentration),
],
)
# Sometimes it warns that it may not have found a good solution,
# but the solution is very accurate.
if optimization_result.success or optimization_result.fun < 1e-05:
break

assert optimization_result.success or optimization_result.fun < 1e-05

result = np.exp(optimization_result.x)
Expand Down
4 changes: 2 additions & 2 deletions integration_tests/test_domestic_migration.py
Original file line number Diff line number Diff line change
Expand Up @@ -114,7 +114,7 @@ def test_individuals_move_into_new_households(
"Domestic migration joining recently created households",
# This case -- people who joined a just-created household -- should be exceedingly rare.
joining_recently_created_households,
target_value_lb=0,
target_value_lb=(1 / 1_000_000_000),
target_value_ub=0.001,
name_addl=f"Time step {time_step}",
)
Expand Down Expand Up @@ -160,7 +160,7 @@ def test_individuals_move_into_new_households(
"Domestic migration joining recently created households",
# This case -- people who joined a just-created household -- should be exceedingly rare.
all_time_joining_recently_created_households,
target_value_lb=0,
target_value_lb=(1 / 1_000_000_000),
target_value_ub=0.001,
name_addl="All time steps",
)
Expand Down

0 comments on commit 3a3b80c

Please # to comment.