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