Use estimagic for NFXP (#42)

Use estimagic to wrap scipy optimizers for the likelihood optimization in the NFXP.

Author: Pascalheid

.conda/meta.yaml

@@ -25,3 +25,4 @@ requirements:
     - pytest
     - pytest-xdist
     - scipy
+    - estimagic >=0.0.30

.gitignore

@@ -21,4 +21,5 @@ figures
 ruspy.rst
 modules.rst
 _generated
+*.db
 

environment.yml

@@ -1,5 +1,7 @@
 name: ruspy
-
+channels:
+  - opensourceeconomics
+  - conda-forge
 dependencies:
   - python=3.7
   - numpy
@@ -16,9 +18,10 @@ dependencies:
   - numba
   - mpmath
   - pip
+  - estimagic>=0.0.30
   - pip:
       - flake8
       - black
       - pre-commit
       - jupyter_contrib_nbextensions
-      - sphinx_rtd_theme
+      - sphinx_rtd_theme

ruspy/estimation/bootstrapping.py

@@ -58,6 +58,7 @@ def loglike_cost_params(
 
 
     """
+    params = params["value"].to_numpy()
     costs = calc_obs_costs(num_states, maint_func, params, scale)
 
     ev = get_ev(params, trans_mat, costs, disc_fac, alg_details)
@@ -109,6 +110,7 @@ def derivative_loglike_cost_params(
 
 
     """
+    params = params["value"].to_numpy()
     dev = np.zeros_like(params)
     obs_costs = calc_obs_costs(num_states, maint_func, params, scale)
 

ruspy/estimation/est_cost_params.py

@@ -2,9 +2,8 @@
 This module contains the main function for the estimation process.
 """
 import numpy as np
-import scipy.optimize as opt
+from estimagic.optimization.optimize import minimize
 
-from ruspy.estimation.bootstrapping import bootstrapp
 from ruspy.estimation.est_cost_params import create_state_matrix
 from ruspy.estimation.est_cost_params import loglike_cost_params
 from ruspy.estimation.estimation_interface import select_model_parameters
@@ -60,32 +59,30 @@ def estimate(init_dict, df):
 
     alg_details = {} if "alg_details" not in init_dict else init_dict["alg_details"]
 
+    kwargs = {
+        "maint_func": maint_func,
+        "maint_func_dev": maint_func_dev,
+        "num_states": num_states,
+        "trans_mat": trans_mat,
+        "state_mat": state_mat,
+        "decision_mat": decision_mat,
+        "disc_fac": disc_fac,
+        "scale": scale,
+        "alg_details": alg_details,
+    }
+
     result_cost_params = {}
 
-    min_result = opt.minimize(
+    min_result = minimize(
         loglike_cost_params,
-        args=(
-            maint_func,
-            maint_func_dev,
-            num_states,
-            trans_mat,
-            state_mat,
-            decision_mat,
-            disc_fac,
-            scale,
-            alg_details,
-        ),
-        **optimizer_options
+        criterion_kwargs=kwargs,
+        gradient_kwargs=kwargs,
+        **optimizer_options,
     )
-    result_cost_params["x"] = min_result["x"]
-    result_cost_params["fun"] = min_result["fun"]
-    result_cost_params["message"] = min_result["message"]
-    result_cost_params["jac"] = min_result["jac"]
-
-    if "hess_inv" in min_result:
-        (
-            result_cost_params["95_conf_interv"],
-            result_cost_params["std_errors"],
-        ) = bootstrapp(min_result["x"], min_result["hess_inv"])
+    result_cost_params["x"] = min_result[1]["value"].to_numpy()
+    result_cost_params["fun"] = min_result[0]["fitness"]
+    result_cost_params["status"] = min_result[0]["status"]
+    result_cost_params["message"] = min_result[0]["message"]
+    result_cost_params["jac"] = min_result[0]["jacobian"]
 
     return transition_results, result_cost_params

ruspy/estimation/estimation.py

@@ -1,4 +1,5 @@
 import numpy as np
+import pandas as pd
 
 from ruspy.estimation.est_cost_params import derivative_loglike_cost_params
 from ruspy.model_code.cost_functions import cubic_costs
@@ -97,46 +98,32 @@ def select_optimizer_options(init_dict, num_params_costs):
 
     """
 
-    optimizer_dict = {} if "optimizer" not in init_dict else init_dict["optimizer"]
-    optimizer_options = {}
+    optimizer_options = {} if "optimizer" not in init_dict else init_dict["optimizer"]
 
-    if "optimizer_name" in optimizer_dict:
-        optimizer_options["method"] = optimizer_dict["optimizer_name"]
+    if "algorithm" not in optimizer_options:
+        optimizer_options["algorithm"] = "scipy_L-BFGS-B"
     else:
-        optimizer_options["method"] = "L-BFGS-B"
+        pass
 
-    if "start_values" in optimizer_dict:
-        optimizer_options["x0"] = np.array(optimizer_dict["start_values"])
-    else:
-        optimizer_options["x0"] = np.power(
-            np.full(num_params_costs, 10, dtype=float),
-            np.arange(1, -num_params_costs + 1, -1),
+    if "params" not in optimizer_options:
+        optimizer_options["params"] = pd.DataFrame(
+            np.power(
+                np.full(num_params_costs, 10, dtype=float),
+                np.arange(1, -num_params_costs + 1, -1),
+            ),
+            columns=["value"],
         )
-
-    if "use_search_bounds" in optimizer_dict:
-        if optimizer_dict["use_search_bounds"] == "yes":
-            optimizer_options["bounds"] = [(1e-6, None)] * num_params_costs
-        else:
-            pass
     else:
         pass
 
-    if "search_bounds" in optimizer_dict:
-        optimizer_options["bounds"] = np.array(optimizer_dict["search_bounds"])
+    if "gradient" not in optimizer_options:
+        optimizer_options["gradient"] = derivative_loglike_cost_params
     else:
         pass
 
-    if "use_gradient" in optimizer_dict:
-        if optimizer_dict["use_gradient"] == "yes":
-            optimizer_options["jac"] = derivative_loglike_cost_params
-        else:
-            pass
+    if "logging" not in optimizer_options:
+        optimizer_options["logging"] = False
     else:
         pass
 
-    if "additional_options" in optimizer_dict:
-        optimizer_options["options"] = optimizer_dict["additional_options"]
-    else:
-        optimizer_options["options"] = {}
-
     return optimizer_options

ruspy/estimation/estimation_interface.py

@@ -4,9 +4,8 @@
 """
 import numba
 import numpy as np
-import scipy.optimize as opt
-
-from ruspy.estimation.bootstrapping import bootstrapp
+import pandas as pd
+from estimagic.optimization.optimize import minimize
 
 
 def estimate_transitions(df):
@@ -29,25 +28,34 @@ def estimate_transitions(df):
     usage = df["usage"].to_numpy(dtype=float)
     usage = usage[~np.isnan(usage)].astype(int)
     result_transitions["trans_count"] = transition_count = np.bincount(usage)
-    raw_result_trans = opt.minimize(
-        loglike_trans,
-        args=transition_count,
-        x0=np.full(len(transition_count), 0.1),
-        method="BFGS",
-    )
-    p_raw = raw_result_trans["x"]
-    result_transitions["x"] = reparam_trans(p_raw)
 
-    result_transitions["95_conf_interv"], result_transitions["std_errors"] = bootstrapp(
-        p_raw, raw_result_trans["hess_inv"], reparam=reparam_trans
+    # Prepare DataFrame for estimagic
+    name = ["trans_prob"]
+    number = np.arange(1, len(transition_count) + 1)
+    index = pd.MultiIndex.from_product([name, number], names=["name", "number"])
+    params = pd.DataFrame(
+        np.full(len(transition_count), 1 / len(transition_count)),
+        columns=["value"],
+        index=index,
+    )
+    constr = [{"loc": "trans_prob", "type": "probability"}]
+
+    raw_result_trans = minimize(
+        criterion=loglike_trans,
+        params=params,
+        algorithm="scipy_L-BFGS-B",
+        constraints=constr,
+        criterion_kwargs={"transition_count": transition_count},
+        logging=False,
     )
 
-    result_transitions["fun"] = loglike_trans(p_raw, transition_count)
+    result_transitions["x"] = raw_result_trans[1]["value"].to_numpy()
+    result_transitions["fun"] = raw_result_trans[0]["fitness"]
 
     return result_transitions
 
 
-def loglike_trans(p_raw, transition_count):
+def loglike_trans(params, transition_count):
     """
     Log-likelihood function of transition probability estimation.
 
@@ -65,31 +73,11 @@ def loglike_trans(p_raw, transition_count):
     log_like : numpy.float
         The negative log-likelihood value of the transition probabilities
     """
-    trans_probs = reparam_trans(p_raw)
-    log_like = -np.sum(np.multiply(transition_count, np.log(trans_probs)))
+    p_raw = params.loc["trans_prob", "value"].to_numpy()
+    log_like = -np.sum(np.multiply(transition_count, np.log(p_raw)))
     return log_like
 
 
-def reparam_trans(p_raw):
-    """
-    The reparametrization function for transition probabilities.
-
-    Parameters
-    ----------
-    p_raw : numpy.array
-        The raw values before reparametrization, on which there are no constraints
-        or bounds.
-
-    Returns
-    -------
-    trans_prob : numpy.array
-        The probabilities of an state increase.
-
-    """
-    p = np.exp(p_raw) / np.sum(np.exp(p_raw))
-    return p
-
-
 @numba.jit(nopython=True)
 def create_transition_matrix(num_states, trans_prob):
     """

ruspy/estimation/estimation_transitions.py

@@ -28,11 +28,7 @@ def inputs():
             "maint_cost_func": "cubic",
             "cost_scale": scale,
         },
-        "optimizer": {
-            "optimizer_name": "BFGS",
-            "use_gradient": "yes",
-            "use_search_bounds": "no",
-        },
+        "optimizer": {"algorithm": "scipy_L-BFGS-B"},
     }
     df = pd.read_pickle(TEST_FOLDER + "group_4.pkl")
     result_trans, result_fixp = estimate(init_dict, df)
@@ -45,7 +41,7 @@ def inputs():
     out["disc_fac"] = disc_fac
     out["num_states"] = num_states
     out["scale"] = scale
-    out["message"] = result_fixp["message"]
+    out["status"] = result_fixp["status"]
     return out
 
 
@@ -56,7 +52,7 @@ def outputs():
     out["params_base"] = np.loadtxt(TEST_FOLDER + "repl_params_cubic.txt")
     out["transition_count"] = np.loadtxt(TEST_FOLDER + "transition_count.txt")
     out["trans_ll"] = 3140.570557
-    out["cost_ll"] = 163.624044  # 162.885
+    out["cost_ll"] = 164.632939  # 162.885
     return out
 
 
@@ -85,7 +81,7 @@ def test_ll_params_derivative(inputs, outputs):
     disc_fac = inputs["disc_fac"]
     assert_array_almost_equal(
         derivative_loglike_cost_params(
-            inputs["params_est"],
+            pd.DataFrame(inputs["params_est"], columns=["value"]),
             cubic_costs,
             cubic_costs_dev,
             num_states,
@@ -102,4 +98,4 @@ def test_ll_params_derivative(inputs, outputs):
 
 
 def test_success(inputs):
-    assert inputs["message"] == "Optimization terminated successfully."
+    assert inputs["status"] == "success"