docs: def calc_bhhh_hessian_approximation(orig_nest_coefs,
                                      index_coefs,
                                      design,
                                      choice_vec,
                                      rows_to_obs,
                                      rows_to_nests,
                                      ridge=None,
                                      weights=None,
                                      use_jacobian=True,
                                      *args,
                                      **kwargs):
      """
      Parameters
      ----------
      orig_nest_coefs : 1D or 2D ndarray.
          All elements should by ints, floats, or longs. If 1D, should have 1
          element for each nesting coefficient being estimated. If 2D, should
          have 1 column for each set of nesting coefficients being used to
          predict the probabilities of each alternative being chosen. There
          should be one row per nesting coefficient. Elements denote the inverse
          of the scale coefficients for each of the lower level nests.
      index_coefs : 1D or 2D ndarray.
          All elements should by ints, floats, or longs. If 1D, should have 1
          element for each utility coefficient being estimated
          (i.e. num_features). If 2D, should have 1 column for each set of
          coefficients being used to predict the probabilities of choosing each
          alternative. There should be one row per index coefficient.
      design : 2D ndarray.
          There should be one row per observation per available alternative.
          There should be one column per utility coefficient being estimated. All
          elements should be ints, floats, or longs.
      choice_vec : 1D ndarray.
          All elements should by ints, floats, or longs. Each element represents
          whether the individual associated with the given row chose the
          alternative associated with the given row.
      rows_to_obs : 2D scipy sparse array.
          There should be one row per observation per available alternative and
          one column per observation. This matrix maps the rows of the design
          matrix to the unique observations (on the columns).
      rows_to_nests : 2D scipy sparse array.
          There should be one row per observation per available alternative and
          one column per nest. This matrix maps the rows of the design matrix to
          the unique nests (on the columns).
      ridge : int, float, long, or None, optional.
          Determines whether or not ridge regression is performed. If an int,
          float or long is passed, then that scalar determines the ridge penalty
          for the optimization. Default == None. Note that if this parameter is
          passed, the values of the BHHH matrix MAY BE INCORRECT since it is not
          100% clear how penalization affects the information matrix.
      use_jacobian : bool, optional.
          Determines whether or not the jacobian will be used when calculating
          the gradient. When performing model estimation, `use_jacobian` should
          be `True` if the values that are actually being estimated are the
          logit of the nest coefficients. Default `== False`.
  
      Returns
      -------
      bhhh_matrix : 2D ndarray.
         The negative of the sum of the outer products of the gradient of the
         log-likelihood function for each observation.
      """
      # Calculate the weights for the sample
      if weights is None:
          weights = np.ones(design.shape[0])
      weights_per_obs = np.max(rows_to_obs.toarray() * weights[:, None], axis=0)
  
      # Transform the nest coefficients into their "always positive" versions
      nest_coefs = naturalize_nest_coefs(orig_nest_coefs)
  
      # Get the vectors and matrices needed to calculate the gradient
      vector_dict = prep_vectors_for_gradient(nest_coefs,
                                              index_coefs,
                                              design,
                                              choice_vec,
                                              rows_to_obs,
                                              rows_to_nests)
  
      # Calculate the index for each alternative for each person
      sys_utility = design.dot(index_coefs)
  
      # Calculate w_ij
      long_w = sys_utility / vector_dict["long_nest_params"]
      # Guard against overflow
      inf_index = np.isposinf(long_w)
      long_w[inf_index] = max_comp_value
  
      ##########
      # Calculate d_log_likelihood_d_nest_params
      ##########
      # Calculate the term that only depends on nest level values
      log_exp_sums = np.log(vector_dict["ind_sums_per_nest"])
      # Guard against overflow
      log_exp_sums[np.isneginf(log_exp_sums)] = -1 * max_comp_value
  
      # Calculate the first term of the derivative of the log-liikelihood
      # with respect to the nest parameters. Note we do not sum this object since
      # we want the values at the 'individual' level, which they already are.
      nest_gradient_term_1 = ((vector_dict["obs_to_chosen_nests"] -
                               vector_dict["nest_choice_probs"]) *
                              log_exp_sums)
  
      # Calculate the second term of the derivative of the log-liikelihood
      # with respect to the nest parameters
      half_deriv = ((vector_dict["long_probs"] -
                     vector_dict["long_chosen_nest"] *
                     vector_dict["prob_given_nest"]) *
                    long_w)[:, None]
      # "Spread out" the second term across the appropriate nests
      spread_half_deriv = rows_to_nests.multiply(half_deriv)
      # Aggregate the spread out half-derivatives to the individual level
      # This object should have shape (num_obs, num_nests)
      nest_gradient_term_2 = rows_to_obs.transpose().dot(spread_half_deriv).A
  
      # Calculate the third term of the derivative of the log-likelihood
      # with respect to the nest parameters
      nest_gradient_term_3a = (choice_vec -
                               vector_dict["long_chosen_nest"] *
                               vector_dict["prob_given_nest"])
  
      nest_gradient_term_3b = ((-1 * nest_gradient_term_3a * long_w) /
                               vector_dict["long_nest_params"])
  
      # Guard against overflow
      inf_idx = np.isposinf(nest_gradient_term_3b)
      nest_gradient_term_3b[inf_idx] = max_comp_value
  
      neg_inf_idx = np.isneginf(nest_gradient_term_3b)
      nest_gradient_term_3b[neg_inf_idx] = -1 * max_comp_value
  
      # Get the nest-wide version of this piece of the gradient
      spread_out_term_3b = rows_to_nests.multiply(nest_gradient_term_3b[:, None])
      nest_gradient_term_3 = rows_to_obs.transpose().dot(spread_out_term_3b).A
  
      # Combine the terms. Note the "nest_coefs * (1 - nest_coefs)" is due to the
      # fact that we're estimating the logit of the nest coefficients instead of
      # the nest coefficient itself. We therefore need to multiply by
      # d_nest_coef_d_estimated_variable to get the correct gradient.
      # d_nest_coef_d_estimated_variable == nest_coefs * (1 - nest_coefs).
      # As with the various nest_gradient_terms, the nest_gradient should be of
      # shape (num_obs, num_nests)
      if use_jacobian:
          jacobian = (nest_coefs * (1.0 - nest_coefs))[None, :]
      else:
          jacobian = 1
  
      nest_gradient = ((nest_gradient_term_1 +
                        nest_gradient_term_2 +
                        nest_gradient_term_3) *
                       jacobian)
  
      ##########
      # Calculate d_loglikelihood_d_beta
      ##########
      beta_gradient_term_1 = (vector_dict["scaled_y"] -
                              vector_dict["p_tilde_given_nest"] +
                              vector_dict["p_tilde_given_nest"] *
                              vector_dict["long_nest_params"] -
                              vector_dict["long_probs"])[:, None]
      #####
      # Calculate the derivative with respect to beta
      #####
      beta_gradient = rows_to_obs.T.dot(beta_gradient_term_1 * design)
  
      #####
      # Combine the gradient pieces
      #####
      gradient_matrix = np.concatenate((nest_gradient, beta_gradient), axis=1)
  
      #####
      # Compute and return the outer product of each row of the gradient
      # with itself. Then sum these individual matrices together. The line below
      # does the same computation just with less memory and time.
      bhhh_matrix =\
          gradient_matrix.T.dot(weights_per_obs[:, None] * gradient_matrix)
  
      if ridge is not None:
          # The rational behind adding 2 * ridge is that the information
          # matrix should approximate the hessian and in the hessian we subtract
          # 2 * ridge at the end. We add 2 * ridge here, since we will multiply
          # by negative one afterwards. I don't know if this is the correct way
          # to calculate the Fisher Information Matrix in ridge regression
          # models.
          bhhh_matrix += 2 * ridge
  
      # Note the "-1" is because we are approximating the Fisher information
      # matrix which has a negative one in the front of it?
      # Note that if we were using the bhhh_matrix to calculate the robust
      # covariance matrix, then we would not multiply by negative one here.
      return -1 * bhhh_matrix