The log of the descendent rate, ln(Rd), is ~ N(mu, bl*ac), where
    the variance = bl*ac = branch_length * autocorrel_param, and mu is set
    so that E[Rd] = Rp:
    E[X] where ln(X) ~ N(mu, sigma^2) = exp(mu+(1/2)*sigma_sq)
    so Rp = exp(mu+(1/2)*bl*ac),
    ln(Rp) = mu + (1/2)*bl*ac,
    ln(Rp) - (1/2)*bl*ac = mu,
    so ln(Rd) ~ N(ln(Rp) - (1/2)*bl*ac, bl*ac)
    (NB: Var[Rd] = Rp^2 * (exp(bl*ac)-1),
         Std[Rd] = Rp * sqrt(exp(bl*ac)-1)

    See: H Kishino, J L Thorne, and W J Bruno (2001)

def logn_correlated_rate(parent_rate, branch_length, autocorrel_param, size=1):
    """
    The log of the descendent rate, ln(Rd), is ~ N(mu, bl*ac), where
    the variance = bl*ac = branch_length * autocorrel_param, and mu is set
    so that E[Rd] = Rp:
    E[X] where ln(X) ~ N(mu, sigma^2) = exp(mu+(1/2)*sigma_sq)
    so Rp = exp(mu+(1/2)*bl*ac),
    ln(Rp) = mu + (1/2)*bl*ac,
    ln(Rp) - (1/2)*bl*ac = mu,
    so ln(Rd) ~ N(ln(Rp) - (1/2)*bl*ac, bl*ac)
    (NB: Var[Rd] = Rp^2 * (exp(bl*ac)-1),
         Std[Rd] = Rp * sqrt(exp(bl*ac)-1)

    See: H Kishino, J L Thorne, and W J Bruno (2001)
    """
    if autocorrel_param <= 0:
        raise Exception('Autocorrelation parameter must be greater than 0')

    variance = branch_length * autocorrel_param
    stdev = np.sqrt(variance)
    ln_descendant_rate = np.random.normal(np.log(parent_rate) - 0.5 * variance,
                                          scale=stdev, size=size)
    descendant_rate = np.exp(ln_descendant_rate)
    return float(descendant_rate) if size == 1 else descendant_rate