Fix bug in likelihood caused by #278.

.pre-commit-config.yaml

@@ -31,7 +31,7 @@ repos:
         additional_dependencies: [
             flake8-bugbear, flake8-builtins, flake8-comprehensions, flake8-docstrings,
             flake8-eradicate, flake8-print, flake8-rst-docstrings, flake8-todo,
-            flake8-unused-arguments, pep8-naming, pydocstyle<4.0,
+            flake8-unused-arguments, pep8-naming, pydocstyle<4.0, Pygments,
         ]
 -   repo: https://github.com/PyCQA/doc8
     rev: 0.8.1rc2

development/testing/regression.py

@@ -8,7 +8,7 @@
 import respy as rp
 from development.testing.notifications import send_notification
 from respy.config import TEST_RESOURCES_DIR
-from respy.config import TOL
+from respy.config import TOL_REGRESSION_TESTS
 from respy.tests.random_model import generate_random_model
 from respy.tests.random_model import simulate_truncated_data
 
@@ -102,7 +102,9 @@ def investigate_regression_test(idx):
 
     crit_val = calc_crit_val(params, options)
 
-    assert np.isclose(crit_val, exp_val, rtol=TOL, atol=TOL)
+    assert np.isclose(
+        crit_val, exp_val, rtol=TOL_REGRESSION_TESTS, atol=TOL_REGRESSION_TESTS
+    )
 
 
 def _check_single(test, strict):
@@ -111,7 +113,9 @@ def _check_single(test, strict):
 
     crit_val = calc_crit_val(params, options)
 
-    is_success = np.isclose(crit_val, exp_val, rtol=TOL, atol=TOL)
+    is_success = np.isclose(
+        crit_val, exp_val, rtol=TOL_REGRESSION_TESTS, atol=TOL_REGRESSION_TESTS
+    )
 
     if strict is True:
         assert is_success, "Failed regression test."

respy/config.py

@@ -12,14 +12,11 @@
 TINY_FLOAT = 1e-8
 PRINT_FLOAT = 1e10
 
-# Number of decimals that are compared for tests This is currently only used in
-# regression tests.
-DECIMALS = 6
 # Some assert functions take rtol instead of decimals
-TOL = 10 ** -DECIMALS
+TOL_REGRESSION_TESTS = 1e-10
 
 # Interpolation
-INADMISSIBILITY_PENALTY = -400000
+INADMISSIBILITY_PENALTY = -400_000
 
 SEED_STARTUP_ITERATION_GAP = 100
 

respy/likelihood.py

@@ -3,8 +3,7 @@
 
 import numba as nb
 import numpy as np
-from scipy.special import logsumexp
-from scipy.special import softmax
+from scipy import special
 
 from respy.conditional_draws import create_draws_and_log_prob_wages
 from respy.config import HUGE_FLOAT
@@ -303,12 +302,12 @@ def _internal_log_like_obs(
     per_individual_loglikes = np.add.reduceat(per_period_loglikes, idx_indiv_first_obs)
     if n_types >= 2:
         z = np.dot(type_covariates, optim_paras["type_prob"].T)
-        type_probabilities = softmax(z, axis=1)
+        type_probabilities = special.softmax(z, axis=1)
 
         log_type_probabilities = np.log(type_probabilities)
         weighted_loglikes = per_individual_loglikes + log_type_probabilities
 
-        contribs = logsumexp(weighted_loglikes, axis=1)
+        contribs = special.logsumexp(weighted_loglikes, axis=1)
     else:
         contribs = per_individual_loglikes.flatten()
 
@@ -318,64 +317,37 @@ def _internal_log_like_obs(
 
 
 @nb.njit
-def log_softmax_i(x, i, tau=1):
-    """Calculate log probability of a soft maximum for index ``i``.
+def logsumexp(x):
+    """Compute `logsumexp(x)` of `x`.
 
-    The log softmax function is essentially
+    The function does the same as the following code, but faster.
 
-    .. math::
-
-        log_softmax_i = log(softmax(x_i))
-
-    This function is superior to the naive implementation as takes advantage of the log
-    and uses the fact that the softmax function is shift-invariant. Integrating the log
-    in the softmax function yields
-
-    .. math::
+    .. code-block:: python
 
-        log_softmax_i = x_i - max(x) - log(sum(exp(x - max(x))))
+        max_x = np.max(x)
+        differences = x - max_x
+        log_sum_exp = max_x + np.log(np.sum(np.exp(differences)))
 
-    The second property ensures that overflows and underflows in the exponential
-    function cannot happen as ``exp(x - max(x)) = exp(0) = 1`` at its highest value and
-    ``exp(-inf) = 0`` at its lowest value. Only infinite inputs can cause an invalid
-    output.
+    The subtraction of the maximum prevents overflows and mitigates the impact of
+    underflows.
 
-    The temperature parameter ``tau`` controls the smoothness of the function. For
-    ``tau`` close to 1, the function is more similar to ``max(x)`` and the derivatives
-    of the functions grow to infinity. As ``tau`` grows, the smoothed maximum is bigger
-    than the actual maximum and derivatives diminish. See [1]_ for more information on
-    the smoothing factor ``tau``. Note that the post covers the inverse of the smoothing
-    factor in this function. It is also covered in [2]_ as the kernel-smoothing choice
-    probability simulator. In reinforcement learning and statistical mechanics the
-    function is known as the Boltzmann exploration.
-
-    .. [1] https://www.johndcook.com/blog/2010/01/13/soft-maximum
-    .. [2] McFadden, D. (1989). A method of simulated moments for estimation of discrete
-           response models without numerical integration. Econometrica: Journal of the
-           Econometric Society, 995-1026.
-
-    Parameters
-    ----------
-    x : np.ndarray
-        Array with shape (n,) containing the values for which a smoothed maximum should
-        be computed.
-    i : int
-        Index for which the log probability should be computed.
-    tau : float
-        Smoothing parameter to control the size of derivatives.
+    """
+    # Search maximum.
+    max_x = None
+    length = len(x)
+    for i in range(length):
+        if max_x is None or x[i] > max_x:
+            max_x = x[i]
 
-    Returns
-    -------
-    log_probability : float
-        Log probability for input value ``i``.
+    # Calculate sum of exponential differences.
+    sum_exp = 0
+    for i in range(length):
+        diff = x[i] - max_x
+        sum_exp += np.exp(diff)
 
-    """
-    max_x = np.max(x)
-    smoothed_differences = (x - max_x) / tau
-    log_sum_exp = np.log(np.sum(np.exp(smoothed_differences)))
-    log_probability = smoothed_differences[i] - log_sum_exp
+    log_sum_exp = max_x + np.log(sum_exp)
 
-    return log_probability
+    return log_sum_exp
 
 
 @nb.guvectorize(
@@ -394,14 +366,27 @@ def simulate_log_probability_of_individuals_observed_choice(
     is_inadmissible,
     choice,
     tau,
-    log_prob_choice,
+    smoothed_log_probability,
 ):
     """Simulate the probability of observing the agent's choice.
 
     The probability is simulated by iterating over a distribution of unobservables.
     First, the utility of each choice is computed. Then, the probability of observing
     the choice of the agent given the maximum utility from all choices is computed.
 
+    The naive implementation calculates the log probability for choice `i` with the
+    softmax function.
+
+    .. math::
+
+        \\log(\text{softmax}(x)_i) = \\log\\left(
+            \frac{e^{x_i}}{\\sum^J e^{x_j}}
+        \right)
+
+    The following function is numerically more robust. The derivation with the two
+    consecutive `logsumexp` functions is included in `#278
+    <https://github.com/OpenSourceEconomics/respy/pull/288>`_.
+
     Parameters
     ----------
     wages : numpy.ndarray
@@ -424,15 +409,15 @@ def simulate_log_probability_of_individuals_observed_choice(
 
     Returns
     -------
-    log_prob_choice : float
-        Smoothed log probability of choice.
+    smoothed_log_probability : float
+        Simulated Smoothed log probability of choice.
 
     """
     n_draws, n_choices = draws.shape
 
     value_functions = np.zeros(n_choices)
 
-    log_prob_choice_ = 0
+    smoothed_log_probabilities = np.empty(n_draws)
 
     for i in range(n_draws):
 
@@ -446,13 +431,15 @@ def simulate_log_probability_of_individuals_observed_choice(
                 is_inadmissible[j],
             )
 
-            value_functions[j] = value_function
+            value_functions[j] = value_function / tau
 
-        log_prob_choice_ += log_softmax_i(value_functions, choice, tau)
+        smoothed_log_probabilities[i] = value_functions[choice] - logsumexp(
+            value_functions
+        )
 
-    log_prob_choice_ /= n_draws
+    smoothed_log_prob = logsumexp(smoothed_log_probabilities) - np.log(n_draws)
 
-    log_prob_choice[0] = log_prob_choice_
+    smoothed_log_probability[0] = smoothed_log_prob
 
 
 def _convert_choice_variables_from_categorical_to_codes(df, options):

respy/tests/test_regression.py

@@ -3,7 +3,7 @@
 import pytest
 
 from development.testing.regression import calc_crit_val
-from respy.config import TOL
+from respy.config import TOL_REGRESSION_TESTS
 
 
 @pytest.mark.parametrize("index", range(10))
@@ -12,4 +12,6 @@ def test_single_regression(regression_vault, index):
     params, options, exp_val = regression_vault[index]
     crit_val = calc_crit_val(params, options)
 
-    assert np.isclose(crit_val, exp_val, rtol=TOL, atol=TOL)
+    assert np.isclose(
+        crit_val, exp_val, rtol=TOL_REGRESSION_TESTS, atol=TOL_REGRESSION_TESTS
+    )