treetime/gtr_site_specific.py
def random(cls, L=1, avg_mu=1.0, alphabet='nuc', pi_dirichlet_alpha=1,
               W_dirichlet_alpha=3.0, mu_gamma_alpha=3.0):
        """
        Creates a random GTR model

        Parameters
        ----------

         mu : float
            Substitution rate

         alphabet : str
            Alphabet name (should be standard: 'nuc', 'nuc_gap', 'aa', 'aa_gap')


        """
        from scipy.stats import gamma
        alphabet=alphabets[alphabet]
        gtr = cls(alphabet=alphabet, seq_len=L)
        n = gtr.alphabet.shape[0]

        if pi_dirichlet_alpha:
            pi = 1.0*gamma.rvs(pi_dirichlet_alpha, size=(n,L))
        else:
            pi = np.ones((n,L))

        pi /= pi.sum(axis=0)
        if W_dirichlet_alpha:
            tmp = 1.0*gamma.rvs(W_dirichlet_alpha, size=(n,n))
        else:
            tmp = np.ones((n,n))
        tmp = np.tril(tmp,k=-1)
        W = tmp + tmp.T

        if mu_gamma_alpha:
            mu = gamma.rvs(mu_gamma_alpha, size=(L,))
        else:
            mu = np.ones(L)

        gtr.assign_rates(mu=mu, pi=pi, W=W)
        gtr.mu *= avg_mu/np.mean(gtr.mu)

        return gtr
def random(cls, L=1, avg_mu=1.0, alphabet='nuc', pi_dirichlet_alpha=1,
               W_dirichlet_alpha=3.0, mu_gamma_alpha=3.0):
        """
        Creates a random GTR model

        Parameters
        ----------

         mu : float
            Substitution rate

         alphabet : str
            Alphabet name (should be standard: 'nuc', 'nuc_gap', 'aa', 'aa_gap')


        """
        from scipy.stats import gamma
        alphabet=alphabets[alphabet]
        gtr = cls(alphabet=alphabet, seq_len=L)
        n = gtr.alphabet.shape[0]

        if pi_dirichlet_alpha:
            pi = 1.0*gamma.rvs(pi_dirichlet_alpha, size=(n,L))
        else:
            pi = np.ones((n,L))

        pi /= pi.sum(axis=0)
        if W_dirichlet_alpha:
            tmp = 1.0*gamma.rvs(W_dirichlet_alpha, size=(n,n))
        else:
            tmp = np.ones((n,n))
        tmp = np.tril(tmp,k=-1)
        W = tmp + tmp.T

        if mu_gamma_alpha:
            mu = gamma.rvs(mu_gamma_alpha, size=(L,))
        else:
            mu = np.ones(L)

        gtr.assign_rates(mu=mu, pi=pi, W=W)
        gtr.mu *= avg_mu/np.mean(gtr.mu)

        return gtr
def random ( cls , L = 1 , avg_mu = 1.0 , alphabet = 'nuc' , pi_dirichlet_alpha = 1 , W_dirichlet_alpha = 3.0 , mu_gamma_alpha = 3.0 ) : from scipy . stats import gamma alphabet = alphabets [ alphabet ] gtr = cls ( alphabet = alphabet , seq_len = L ) n = gtr . alphabet . shape [ 0 ] if pi_dirichlet_alpha : pi = 1.0 * gamma . rvs ( pi_dirichlet_alpha , size = ( n , L ) ) else : pi = np . ones ( ( n , L ) ) pi /= pi . sum ( axis = 0 ) if W_dirichlet_alpha : tmp = 1.0 * gamma . rvs ( W_dirichlet_alpha , size = ( n , n ) ) else : tmp = np . ones ( ( n , n ) ) tmp = np . tril ( tmp , k = - 1 ) W = tmp + tmp . T if mu_gamma_alpha : mu = gamma . rvs ( mu_gamma_alpha , size = ( L , ) ) else : mu = np . ones ( L ) gtr . assign_rates ( mu = mu , pi = pi , W = W ) gtr . mu *= avg_mu / np . mean ( gtr . mu ) return gtr
Creates a random GTR model

        Parameters
        ----------

         mu : float
            Substitution rate

         alphabet : str
            Alphabet name (should be standard: 'nuc', 'nuc_gap', 'aa', 'aa_gap')
Creates a random GTR model