docs: def get_param_list_for_prediction(model_obj, replicates):
      """
      Create the `param_list` argument for use with `model_obj.predict`.
  
      Parameters
      ----------
      model_obj : an instance of an MNDC object.
          Should have the following attributes:
          `['ind_var_names', 'intercept_names', 'shape_names', 'nest_names']`.
          This model should have already undergone a complete estimation process.
          I.e. its `fit_mle` method should have been called without
          `just_point=True`.
      replicates : 2D ndarray.
          Should represent the set of parameter values that we now wish to
          partition for use with the `model_obj.predict` method.
  
      Returns
      -------
      param_list : list.
          Contains four elements, each being a numpy array. Either all of the
          arrays should be 1D or all of the arrays should be 2D. If 2D, the
          arrays should have the same number of columns. Each column being a
          particular set of parameter values that one wants to predict with.
          The first element in the list should be the index coefficients. The
          second element should contain the 'outside' intercept parameters if
          there are any, or None otherwise. The third element should contain
          the shape parameters if there are any or None otherwise. The fourth
          element should contain the nest coefficients if there are any or
          None otherwise. Default == None.
      """
      # Check the validity of the passed arguments
      ensure_samples_is_ndim_ndarray(replicates, ndim=2, name='replicates')
      # Determine the number of index coefficients, outside intercepts,
      # shape parameters, and nest parameters
      num_idx_coefs = len(model_obj.ind_var_names)
  
      intercept_names = model_obj.intercept_names
      num_outside_intercepts =\
          0 if intercept_names is None else len(intercept_names)
  
      shape_names = model_obj.shape_names
      num_shapes = 0 if shape_names is None else len(shape_names)
  
      nest_names = model_obj.nest_names
      num_nests = 0 if nest_names is None else len(nest_names)
  
      parameter_numbers =\
          [num_nests, num_shapes, num_outside_intercepts, num_idx_coefs]
      current_idx = 0
      param_list = []
      for param_num in parameter_numbers:
          if param_num == 0:
              param_list.insert(0, None)
              continue
          upper_idx = current_idx + param_num
          param_list.insert(0, replicates[:, current_idx:upper_idx].T)
          current_idx += param_num
      return param_list