Description

Some tests on simulating brownian motion evolution

brownian motion model:

dX(t) = sigma d B(t)

  • sigma = standard deviation
  • B = random noise
  • t = time

In [2]:
import dendropy
from scipy.stats import norm

In [12]:
def brownian(x0, n, dt, delta):
    for i in xrange(n):
        x0 += np.random.normal(scale=delta**2*dt)
    return x0

brownian(0.5, 10, 0.5, 0.25)


Out[12]:
0.44491847561971076

In [4]:
import random
import dendropy

def process_node(node, start=1.0):
    if node.parent_node is None:
        node.value = start
    else:
        node.value = random.gauss(node.parent_node.value, node.edge.length)
    for child in node.child_nodes():
        process_node(child)
    if node.taxon is not None:
        print("%s : %s" % (node.taxon, node.value))

        
mle = dendropy.treesim.birth_death(birth_rate=1, death_rate=0.5, ntax=10)
process_node(mle.seed_node)


T1 : -0.541148641712
T2 : -0.616271210654
T3 : 0.653110263093
T4 : 0.962124941326
T5 : -0.36967857664
T6 : 1.69251796248
T7 : 1.71744530603
T8 : -0.00492935851821
T9 : -0.176923455264
T10 : -0.00337972505324

In [44]:
from math import sqrt
from scipy.stats import norm
import numpy as np
     
     
def brownian(x0, n, dt, delta):
    """\
    Generate an instance of Brownian motion (i.e. the Wiener process):
     
        X(t) = X(0) + N(0, delta**2 * t; 0, t)
     
    where N(a,b; t0, t1) is a normally distributed random variable with mean a and
    variance b.  The parameters t0 and t1 make explicit the statistical
    independence of N on different time intervals; that is, if [t0, t1) and
    [t2, t3) are disjoint intervals, then N(a, b; t0, t1) and N(a, b; t2, t3)
    are independent.
    
    Arguments
    ---------
    x0 : float or numpy array (or something that can be converted to a numpy array
         using numpy.asarray(x0)).
        The initial condition(s) (i.e. position(s)) of the Brownian motion.
    n : int
        The number of steps to take.
    dt : float
        The time step.
    delta : float
        delta determines the "speed" of the Brownian motion.  The random variable
        of the position at time t, X(t), has a normal distribution whose mean is
        the position at time t=0 and whose variance is delta**2*t.

    Returns
    -------
    A numpy array of floats with shape `x0.shape + (n,)`.
    
    Note that the initial value `x0` is not included in the returned array.
    """
     
    x0 = np.asarray(x0)
     
    # For each element of x0, generate a sample of n numbers from a
    # normal distribution.
    r = np.random.normal(size=x0.shape + (n,), scale=delta*np.sqrt(dt))
     
    # If `out` was not given, create an output array.
    if out is None:
        out = np.empty(r.shape)
     
    # This computes the Brownian motion by forming the cumulative sum of
    # the random samples. 
    np.cumsum(r, axis=-1, out=out)
     
    # Add the initial condition.
    out += np.expand_dims(x0, axis=-1)
     
    return out

for i in xrange(10):
    print brownian(0.5, 1, 1, 0.25)


[ 0.88498692]
[ 0.70247526]
[ 0.64399443]
[ 0.84922204]
[ 0.34150368]
[-0.03608391]
[ 0.28334593]
[ 0.60437544]
[ 0.63495074]
[ 0.55608189]

In [50]:
import random
import dendropy

def process_node(node, start=1.0):
    if node.parent_node is None:
        node.value = start
    else:
        x = brownian(node.parent_node.value, 
                                      n = 1,
                                      dt = node.edge.length,
                                      delta = 0.25)
        x = float(x[-1])
        x = x if x >=0 else 0
        x = x if x <= 1 else 1
        node.value = x
    for child in node.child_nodes():
        process_node(child)
    if node.taxon is not None:
        print("%s : %s" % (node.taxon, node.value))

        
mle = dendropy.treesim.birth_death(birth_rate=1, death_rate=0.5, ntax=10)
process_node(mle.seed_node)


T1 : 1
T2 : 0.778752702556
T3 : 0.839397868542
T4 : 0.910064289584
T5 : 1
T6 : 1
T7 : 1
T8 : 0.890143644791
T9 : 1
T10 : 0.664577232122

brownian motion and purely random

  • function composed of both brownian motion and purely random selection
  • Idea from paper: "How to measure and test phylogenetic signal"

  • a ratio parameter determines how much of the brownian motion vs random continuous value is use

    • range of 0-1, 0 = random, 1 = BM
    • BD_value ratio + random_value (1-ratio)

In [158]:
import numpy as np
import scipy.stats as stats
import dendropy

def sim_trait(node, start=0, sigma=0.1, ratio=0.5, verbose=False):
    if node.parent_node is None:
        node.value = start
    else:
        BM = np.random.normal(loc=node.parent_node.value, scale=sigma)
        #rnd = np.random.uniform(minVal, maxVal) 
        rnd = np.random.normal(loc=start, scale=sigma) 
        node.value = BM * ratio + rnd * (1 - ratio)
        #print([BM, rnd, node.value])
        #node.value = node.value if node.value >= minVal else minVal
        #node.value = node.value if node.value <= maxVal else maxVal
    for child in node.child_nodes():
        sim_trait(child, start=start, sigma=sigma,
                  ratio=ratio, verbose=verbose)
    if verbose and node.taxon is not None:
        print('{} : {}'.format(node.taxon, node.value))

        
mle = dendropy.treesim.birth_death(birth_rate=1, death_rate=0.5, ntax=10)
sim_trait(mle.seed_node, verbose=True)
mle.print_plot(display_width=70)


T1 : 0.00578562556579
T2 : 0.158996438024
T3 : 0.189514631499
T4 : 0.186688554744
T5 : -0.0387722622727
T6 : 0.0052849461139
T7 : 0.0860671259267
T8 : 0.0135340232601
T9 : -0.0110602888484
T10 : -0.099397459099
                          /--------------------------------------- T1 
/-------------------------+                                           
|                         |            /-------------------------- T2 
|                         \------------+                              
|                                      |            /------------- T3 
|                                      \------------+                 
|                                                   \------------- T4 
+                                                                     
|                                      /-------------------------- T5 
|                         /------------+                              
|                         |            |            /------------- T6 
|            /------------+            \------------+                 
|            |            |                         \------------- T7 
|            |            |                                           
\------------+            \--------------------------------------- T8 
             |                                                        
             |                                      /------------- T9 
             \--------------------------------------+                 
                                                    \------------- T10
                                                                      
                                                                      

2nd attempt


In [ ]:
from math import sqrt
from scipy.stats import norm
import numpy as np
     
     
def brownian(x0, n, dt, delta):
    """\
    Generate an instance of Brownian motion (i.e. the Wiener process):
     
        X(t) = X(0) + N(0, delta**2 * t; 0, t)
     
    where N(a,b; t0, t1) is a normally distributed random variable with mean a and
    variance b.  The parameters t0 and t1 make explicit the statistical
    independence of N on different time intervals; that is, if [t0, t1) and
    [t2, t3) are disjoint intervals, then N(a, b; t0, t1) and N(a, b; t2, t3)
    are independent.
    
    Arguments
    ---------
    x0 : float or numpy array (or something that can be converted to a numpy array
         using numpy.asarray(x0)).
        The initial condition(s) (i.e. position(s)) of the Brownian motion.
    n : int
        The number of steps to take.
    dt : float
        The time step.
    delta : float
        delta determines the "speed" of the Brownian motion.  The random variable
        of the position at time t, X(t), has a normal distribution whose mean is
        the position at time t=0 and whose variance is delta**2*t.

    Returns
    -------
    A numpy array of floats with shape `x0.shape + (n,)`.
    
    Note that the initial value `x0` is not included in the returned array.
    """
     
    x0 = np.asarray(x0)
     
    # For each element of x0, generate a sample of n numbers from a
    # normal distribution.
    r = np.random.normal(size=x0.shape + (n,), scale=delta*np.sqrt(dt))
     
    # If `out` was not given, create an output array.
    if out is None:
        out = np.empty(r.shape)
     
    # This computes the Brownian motion by forming the cumulative sum of
    # the random samples. 
    np.cumsum(r, axis=-1, out=out)
     
    # Add the initial condition.
    out += np.expand_dims(x0, axis=-1)
     
    return out

for i in xrange(10):
    print brownian(0.5, 1, 1, 0.25)

In [167]:
import numpy as np
import dendropy

def sim_traits(tree, start=0, sigma=0.1, weight=0.5, verbose=False):
    """Trait simulation as detailed in:
    author = {Münkemüller, Tamara and Lavergne, Sébastien and Bzeznik,
              Bruno and Dray, Stéphane and Jombart,
              Thibaut and Schiffers, Katja and Thuiller, Wilfried},
    title = {How to measure and test phylogenetic signal},
    journal = {Methods in Ecology and Evolution}
    
    Args:
    tree -- dendropy tree object
    start -- starting value for continuous character evolution
    sigma -- sigma use for drawing from a normal distribution
    weight -- weight parameter for random vs Brownian motion
              range: 0-1; 0 = purely random; 1 = purely Brownian
    verbose -- verbose output
    """
    ntaxa = len(tree.nodes())
    # simulate brownian motion
    BM = np.random.normal(loc=0, scale=sigma, size=ntaxa)
    BM = np.cumsum(BM) + start
    # random values
    rnd = np.random.permutation(BM)
    # making weighted sums
    ws = weight * BM + (1-weight) * rnd
    # z-scaling weighted sums
    ws = (ws - np.mean(ws)) / np.std(ws)
    
    for i, node in enumerate(tree.preorder_node_iter()):
        node.value = ws[i]
        if verbose and node.taxon is not None:
            print('{} : {}'.format(node.taxon, node.value))
            
    

mle = dendropy.treesim.birth_death(birth_rate=1, death_rate=0.5, ntax=10)
sim_traits(mle, verbose=True)
mle.print_plot(display_width=70)


T1 : -0.891054451744
T2 : -0.322280465379
T3 : -2.44657691116
T4 : -0.0740668194434
T5 : 0.997935100731
T6 : 1.81915911621
T7 : 0.884005888766
T8 : 0.0804171223293
T9 : -0.0614338765143
T10 : 0.636976007026
                                                       /---------- T1 
                      /--------------------------------+              
                      |                                \---------- T2 
                      |                                               
           /----------+                                /---------- T3 
           |          |          /---------------------+              
           |          |          |                     \---------- T4 
           |          \----------+                                    
           |                     |          /--------------------- T5 
/----------+                     \----------+                         
|          |                                |          /---------- T6 
|          |                                \----------+              
|          |                                           \---------- T7 
+          |                                                          
|          \------------------------------------------------------ T8 
|                                                                     
|                                                      /---------- T9 
\------------------------------------------------------+              
                                                       \---------- T10
                                                                      
                                                                      

In [169]:
mle = dendropy.treesim.birth_death(birth_rate=1, death_rate=0.5, ntax=10)
sim_traits(mle, weight=1, verbose=True)
mle.print_plot(display_width=70)


T1 : 1.50565341718
T2 : 0.508000165942
T3 : 0.144948439186
T4 : 0.84560056451
T5 : 0.39228645008
T6 : -0.0694699478639
T7 : -1.25830579824
T8 : -1.63654304486
T9 : -0.610682714374
T10 : -1.32751346635
                /------------------------------------------------- T1 
/---------------+                                                     
|               |                /-------------------------------- T2 
|               \----------------+                                    
|                                |               /---------------- T3 
|                                \---------------+                    
+                                                \---------------- T4 
|                                                                     
|                                                /---------------- T5 
|               /--------------------------------+                    
|               |                                \---------------- T6 
|               |                                                     
\---------------+                                /---------------- T7 
                |                /---------------+                    
                |                |               \---------------- T8 
                \----------------+                                    
                                 |               /---------------- T9 
                                 \---------------+                    
                                                 \---------------- T10
                                                                      
                                                                      

In [170]:
mle = dendropy.treesim.birth_death(birth_rate=1, death_rate=0.5, ntax=10)
sim_traits(mle, weight=0, verbose=True)
mle.print_plot(display_width=70)


T1 : -1.35914169871
T2 : -0.154150467806
T3 : -1.40942068381
T4 : -0.586737000775
T5 : 1.2507582239
T6 : -1.17633229963
T7 : -0.223550833593
T8 : 0.442030123667
T9 : 0.556943323431
T10 : -2.51103996866
                      /------------------------------------------- T1 
                      |                                               
                      |                                /---------- T2 
           /----------+                     /----------+              
           |          |          /----------+          \---------- T3 
           |          |          |          |                         
           |          \----------+          \--------------------- T4 
/----------+                     |                                    
|          |                     |                     /---------- T5 
|          |                     \---------------------+              
|          |                                           \---------- T6 
+          |                                                          
|          \------------------------------------------------------ T7 
|                                                                     
|                                                      /---------- T8 
|                                           /----------+              
\-------------------------------------------+          \---------- T9 
                                            |                         
                                            \--------------------- T10
                                                                      
                                                                      

Dendropy sandbox


In [171]:
tree = dendropy.treesim.birth_death(birth_rate=1, death_rate=0.5, ntax=10)
tree.find_node

In [178]:
chrMtx = dendropy.ContinuousCharacterMatrix()
chrMtx.extend?

In [177]: