docs: def get_pij_method(model=F81, frequencies=None, kappa=None):
      """
      Returns a function for calculation of probability matrix of substitutions i->j over time t.
  
      :param kappa: kappa parameter for HKY model
      :type kappa: float
      :param frequencies: array of state frequencies \pi_i
      :type frequencies: numpy.array
      :param model: model of character evolution
      :type model: str
      :return: probability matrix
      :rtype: function
      """
      if is_f81_like(model):
          mu = get_mu(frequencies)
          return lambda t: get_f81_pij(t, frequencies, mu)
      if JTT == model:
          return get_jtt_pij
      if HKY == model:
          return lambda t: get_hky_pij(t, frequencies, kappa)