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)

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)

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]:

``````