docs: def _convert_eta_to_c(eta, ref_position):
      """
      Parameters
      ----------
      eta : 1D or 2D ndarray.
          The elements of the array should be this model's 'transformed' shape
          parameters, i.e. the natural log of (the corresponding shape parameter
          divided by the reference shape parameter). This array's elements will
          be real valued. If `eta` is 2D, then its shape should be
          (num_estimated_shapes, num_parameter_samples).
      ref_position : int.
          Specifies the position in the resulting array of shape ==
          `(eta.shape[0] + 1,)` that should be equal to 1 - the sum of the other
          elements in the resulting array.
  
      Returns
      -------
      c_vector : 1D or 2D ndarray based on `eta`.
          If `eta` is 1D then `c_vector` should have shape
          `(eta.shape[0] + 1, )`. If `eta` is 2D then `c_vector` should have
          shape `(eta.shape[0] + 1, eta.shape[1])`. The returned array will
          contains the 'natural' shape parameters that correspond to `eta`.
      """
      # Exponentiate eta
      exp_eta = np.exp(eta)
  
      # Guard against overflow
      exp_eta[np.isposinf(exp_eta)] = max_comp_value
      # Guard against underflow
      exp_eta[exp_eta == 0] = min_comp_value
  
      # Calculate the denominator in a logistic transformation
      # Note the +1 is for the reference alternative which has been
      # constrained so that its corresponding eta = 0 and exp(0) = 1
      denom = exp_eta.sum(axis=0) + 1
  
      # Get a list of all the indices (or row indices) corresponding to the
      # alternatives whose shape parameters are being estimated.
      replace_list = list(range(eta.shape[0] + 1))
      replace_list.remove(ref_position)
  
      # Initialize an array for the vector of shape parameters, c
      if len(eta.shape) > 1 and eta.shape[1] > 1:
          # Get an array of zeros with shape
          # (num_possible_alternatives, num_parameter_samples). This is used when
          # working with samples from a Bayesian posterior distribution
          c_vector = np.zeros((eta.shape[0] + 1,
                               eta.shape[1]))
  
          # Calculate the natural shape parameters
          c_vector[replace_list, :] = exp_eta / denom
          c_vector[ref_position, :] = 1.0 / denom
      else:
          # Get an array of zeros with shape (num_possible_alternatives,)
          c_vector = np.zeros(eta.shape[0] + 1)
  
          # Calculate the natural shape parameters
          c_vector[replace_list] = exp_eta / denom
          c_vector[ref_position] = 1.0 / denom
  
      return c_vector