A note and implementation by python for Barton's coalescent in plane¶
Nicholas H. Barton, Jerome Kelleher, Alison M. Etheridge, A NEW MODEL FOR EXTINCTION AND RECOLONIZATION IN TWO DIMENSIONS: QUANTIFYING PHYLOGEOGRAPHY, Evolution, Volume 64, Issue 9, 1 September 2010, Pages 2701–2715, https://doi.org/10.1111/j.1558-5646.2010.01019.x
By: Longxiao 13.03.2025
Forward process¶
We firstly go forward step by step:
Initialize parameters, include population density $\rho$, rate of event $\Lambda$, total death rate $u_0$, variance of events distribution $\theta$, and relative strength of parent seeking to extinction $\alpha$.
With out loss of generality, we set the region as $1 \times 1$ square. Random seed is set as $123$. And we use two list $\texttt(population\_ history)$ and $\texttt(events)$ to record the population changes and events.
In [1]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.cm as cm
#initialize parameters
rho=100
Lambda=1
u0=0.8
theta=0.1
alpha=1
space_size=(1, 1)
seed=123
Poisson distributed initial population¶
Here we initialize the population. As said in the article, individuals are scattered over a continuous two-dimensional range in a uniform Poisson distribution with density $\rho$. As a result, the initial population size follows a Poisson distribution of $\rho \times area$. And individuals are distributed randomly.
In [ ]:
'''initialize population'''
population_history = []
events = []
area = space_size[0] * space_size[1]
num_individuals = np.random.poisson(rho * area)
x = np.random.uniform(0, space_size[0], num_individuals)
y = np.random.uniform(0, space_size[1], num_individuals)
population = pd.DataFrame({
"id": np.arange(1, num_individuals + 1),
"x": x,
"y": y,
"parent_id": None,
"ancestor_id": np.arange(1, num_individuals + 1)
})
population_history.append(population.copy())
population
Out[ ]:
| id | x | y | parent_id | ancestor_id | |
|---|---|---|---|---|---|
| 0 | 1 | 0.904860 | 0.693196 | None | 1 |
| 1 | 2 | 0.455613 | 0.019816 | None | 2 |
| 2 | 3 | 0.955966 | 0.090190 | None | 3 |
| 3 | 4 | 0.590997 | 0.877293 | None | 4 |
| 4 | 5 | 0.199621 | 0.463968 | None | 5 |
| ... | ... | ... | ... | ... | ... |
| 88 | 89 | 0.667488 | 0.330640 | None | 89 |
| 89 | 90 | 0.606192 | 0.910983 | None | 90 |
| 90 | 91 | 0.103685 | 0.166975 | None | 91 |
| 91 | 92 | 0.517172 | 0.273896 | None | 92 |
| 92 | 93 | 0.014645 | 0.172948 | None | 93 |
93 rows × 5 columns
Visualize the initial population distribution:
In [77]:
plt.figure(figsize=(4, 4))
plt.xlim(0, space_size[0])
plt.ylim(0, space_size[1])
plt.scatter(population["x"], population["y"],
s=40, c='blue', alpha=0.7,
edgecolors='white', linewidths=0.5,
label='Existing Individuals'
)
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0)
plt.title("initial")
plt.show()
Randomly happend event¶
Then, the events happen with the rate $\Lambda$, and the center of a event is draw randomly. We use teh function event_gen to generate time interval and the center of a event:
In [78]:
'''a event happens'''
def event_gen():
dt = np.random.exponential(1 / Lambda)
z_x = np.random.uniform(0, space_size[0])
z_y = np.random.uniform(0, space_size[1])
return dt, z_x, z_y
dt, z_x, z_y = event_gen()
print(dt, z_x, z_y)
0.23909670589208998 0.37200613069215527 0.5477870709778582
Visualize the center of an event
In [79]:
plt.figure(figsize=(4, 4))
plt.xlim(0, space_size[0])
plt.ylim(0, space_size[1])
plt.scatter(population["x"], population["y"],
s=40, c='blue', alpha=0.7,
edgecolors='white', linewidths=0.5,
label='Existing Individuals'
)
plt.scatter(z_x, z_y,
s=20, marker='x', c='red',
linewidths=1, zorder=1,
label='Event Center'
)
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0)
plt.title("event")
plt.show
Out[79]:
<function matplotlib.pyplot.show(close=None, block=None)>
Death and birth happend around the event center¶
Death caused by the event is happened with a probability related to the distance between individual and event center. Specifically, it follows $$ u(z,x) = u_0 e^{-\frac{|z-x|^2}{2\theta^2}} $$ where $z$ is the center and $x$ is the position of a individual. Function ind_death is used to draw individuals that would be dead in the certain event happened in $(z_x, z_y)$.
In [80]:
def ind_death(z_x, z_y,population):
dx = population["x"] - z_x
dy = population["y"] - z_y
squared_dist = dx**2 + dy**2
death_probs = u0 * np.exp(-squared_dist / (2 * theta**2))
dead_mask = np.random.rand(len(population)) < death_probs
survivors = population[~dead_mask].copy()
survivors["parent_id"] = survivors["id"]
death = population[dead_mask]
return survivors, death
survivors, death = ind_death(z_x, z_y, population)
Visualize dead individuals
In [81]:
plt.figure(figsize=(4, 4))
plt.xlim(0, space_size[0])
plt.ylim(0, space_size[1])
plt.scatter(population["x"], population["y"],
s=40, c='blue', alpha=0.7,
edgecolors='white', linewidths=0.5,
label='Existing Individuals'
)
plt.scatter(z_x, z_y,
s=20, marker='x', c='red',
linewidths=1, zorder=1,
label='Event Center'
)
plt.scatter(death['x'], death['y'],
linestyle=':', c='white',
edgecolor='black', linewidth=1.2
)
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0)
plt.title("event")
plt.show
Out[81]:
<function matplotlib.pyplot.show(close=None, block=None)>
The new parent, or the individual reproduce in the event is choosen also follows a normal distribution centered by the event center, but with a variance $(\alpha \cdot \theta)^2$. Basically, it is: $$ v(z,y) = u_0 e^{-\frac{|z-y|^2}{2 (\alpha \cdot \theta)^2}} $$ where $y$ is the position of the new parent.
In [82]:
def ind_parent(z_x, z_y, population):
s_x = population["x"].values
s_y = population["y"].values
s_dx = s_x - z_x
s_dy = s_y - z_y
s_dist = s_dx**2 + s_dy**2
parent_weights = np.exp(-s_dist / (2 * (alpha * theta)**2))
parent_weights /= parent_weights.sum()
parent_idx = np.random.choice(len(population), p=parent_weights)
parent = population.iloc[parent_idx]
return parent
parent = ind_parent(z_x, z_y, population)
Visualize the choosen parent individual
In [83]:
plt.figure(figsize=(4, 4))
plt.xlim(0, space_size[0])
plt.ylim(0, space_size[1])
plt.scatter(population["x"], population["y"],
s=40, c='blue', alpha=0.7,
edgecolors='white', linewidths=0.5,
label='Existing Individuals'
)
plt.scatter(z_x, z_y,
s=20, marker='x', c='red',
linewidths=1, zorder=1,
label='Event Center'
)
plt.scatter(parent['x'], parent['y'],
c='red',
edgecolor='white', linewidth=1.2,
label='parent'
)
plt.scatter(death['x'], death['y'],
linestyle=':', c='white', alpha=0.5,
edgecolor='black', linewidth=1.5,
label='death'
)
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0)
plt.title("event")
plt.show
Out[83]:
<function matplotlib.pyplot.show(close=None, block=None)>
Then, new individuals are shows up with a density $\rho u(z,x)$. That means, the number of new individuals follows Possion distribution with parameter $\rho \cdot \int u(z,x) {\rm d}x = 2 \rho u_0 \pi \theta ^2$, and sitributed around the event center follows normal distribution as used above.
In [86]:
def ind_new(z_x, z_y):
event_area = 2 * np.pi * theta**2
expected_new = rho * event_area * u0
num_new = np.random.poisson(expected_new)
new_x = np.random.normal(z_x, theta, num_new)
new_y = np.random.normal(z_y, theta, num_new)
new_x = new_x[(new_x >= 0)& (new_x <= space_size[0])]
new_y = new_x[(new_x >= 0)& (new_x <= space_size[0])]
return new_x, new_y
new_x, new_y = ind_new(z_x, z_y)
Visualize new individuals
In [87]:
plt.figure(figsize=(4, 4))
plt.xlim(0, space_size[0])
plt.ylim(0, space_size[1])
plt.scatter(population["x"], population["y"],
s=40, c='blue', alpha=0.7,
edgecolors='white', linewidths=0.5,
label='Existing Individuals'
)
plt.scatter(z_x, z_y,
s=20, marker='x', c='red',
linewidths=1, zorder=1,
label='Event Center'
)
plt.scatter(parent['x'], parent['y'],
c='red',
edgecolor='white', linewidth=0.5,
label='parent'
)
plt.scatter(death['x'], death['y'],
linestyle=':', c='white', alpha=0.5,
edgecolor='black', linewidth=1.5,
label='death'
)
plt.scatter(new_x, new_y,
c='yellow', alpha=0.5,
edgecolor='white', linewidth=0.5,
label='new_ind'
)
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0)
plt.title("event")
Out[87]:
Text(0.5, 1.0, 'event')
An event finished. We move the dead individuals and add new individuals to the population, and record population historical state.
In [88]:
def new_population(new_x, new_y, parent, survivors):
new_individuals = pd.DataFrame({
"id": None,
"x": new_x,
"y": new_y,
"parent_id": parent['id'],
"ancestor_id": parent['ancestor_id']
})
population = pd.concat([survivors, new_individuals], ignore_index=True)
population['id'] = np.arange(1, len(population) + 1)
return population
population = new_population(new_x, new_y, parent, survivors)
In [89]:
plt.figure(figsize=(4, 4))
plt.xlim(0, space_size[0])
plt.ylim(0, space_size[1])
plt.scatter(population["x"], population["y"],
s=40, c='blue', alpha=0.7,
edgecolors='white', linewidths=0.5,
label='Existing Individuals'
)
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0)
plt.title("1-step")
plt.show()
Longer simulation¶
For simulation, we now define a function event, to iterate on population
In [ ]:
def event(population):
dt, z_x, z_y = event_gen()
survivors, death = ind_death(z_x, z_y, population)
parent = ind_parent(z_x, z_y,population)
new_x, new_y = ind_new(z_x,z_y)
population = new_population(new_x, new_y, parent, survivors)
events.append({
"time": dt,
"z_x": z_x,
"z_y": z_y,
"num_dead": len(death),
"num_new": len(new_x)
})
population_history.append(population.copy())
return population
population = event(population)
population
Out[ ]:
| id | x | y | parent_id | ancestor_id | |
|---|---|---|---|---|---|
| 0 | 1 | 0.904860 | 0.693196 | 1.0 | 1.0 |
| 1 | 2 | 0.455613 | 0.019816 | 2.0 | 2.0 |
| 2 | 3 | 0.955966 | 0.090190 | 3.0 | 3.0 |
| 3 | 4 | 0.590997 | 0.877293 | 4.0 | 4.0 |
| 4 | 5 | 0.935941 | 0.727350 | 5.0 | 6.0 |
| ... | ... | ... | ... | ... | ... |
| 95 | 96 | 0.392644 | 0.392644 | 29.0 | 32.0 |
| 96 | 97 | 0.348024 | 0.348024 | 29.0 | 32.0 |
| 97 | 98 | 0.322529 | 0.322529 | 29.0 | 32.0 |
| 98 | 99 | 0.377790 | 0.377790 | 29.0 | 32.0 |
| 99 | 100 | 0.342678 | 0.342678 | 29.0 | 32.0 |
100 rows × 5 columns
population density kept¶
Then, we could simulate the population in limited time scale (for example, total time is 10). The population history and events will be recorded by the two list given before.
In [96]:
def simulate(population, total_time=100):
"""simulate until time limitation"""
current_time = 0
event_id = 0
while current_time < total_time:
population = event(population)
current_time += events[event_id]['time']
event_id += 1
simulate(population)
len(population_history),len(events)
Out[96]:
(477, 476)
In [97]:
pop_size = np.zeros(len(population_history))
for i in range(len(population_history)):
pop_size[i] = len(population_history[i])
print(pop_size)
plt.plot(pop_size/rho)
[ 93. 100. 96. 95. 91. 91. 89. 84. 78. 80. 87. 88. 93. 95. 93. 95. 94. 97. 96. 101. 105. 107. 95. 99. 89. 91. 90. 91. 95. 94. 96. 96. 103. 107. 112. 113. 113. 116. 121. 129. 114. 118. 109. 114. 117. 121. 124. 132. 128. 131. 130. 125. 127. 102. 84. 83. 84. 95. 97. 98. 92. 94. 93. 95. 85. 85. 84. 84. 86. 85. 87. 90. 88. 94. 96. 98. 103. 103. 99. 104. 107. 101. 98. 103. 107. 108. 111. 116. 120. 123. 131. 134. 138. 143. 145. 124. 123. 104. 104. 110. 108. 105. 101. 102. 104. 102. 89. 94. 96. 96. 102. 103. 101. 98. 106. 109. 110. 103. 96. 99. 101. 102. 107. 105. 102. 102. 101. 99. 87. 87. 86. 88. 96. 96. 98. 101. 106. 102. 108. 105. 106. 112. 114. 110. 116. 107. 110. 116. 118. 121. 124. 115. 107. 105. 96. 98. 100. 107. 116. 103. 104. 110. 115. 120. 120. 124. 123. 129. 131. 132. 127. 136. 131. 133. 130. 135. 122. 115. 105. 112. 120. 109. 111. 101. 100. 104. 105. 107. 109. 103. 89. 95. 98. 99. 98. 97. 104. 108. 110. 103. 108. 106. 97. 98. 103. 99. 104. 105. 108. 102. 98. 89. 98. 102. 101. 103. 110. 107. 111. 114. 108. 96. 97. 98. 92. 94. 95. 94. 101. 105. 107. 113. 114. 112. 107. 108. 98. 98. 85. 84. 88. 87. 83. 88. 95. 98. 93. 95. 89. 99. 101. 93. 90. 93. 88. 89. 87. 98. 107. 107. 106. 102. 101. 104. 108. 104. 100. 107. 117. 119. 115. 117. 119. 126. 133. 138. 140. 135. 110. 113. 117. 114. 117. 124. 108. 99. 98. 102. 102. 105. 103. 103. 107. 102. 103. 104. 104. 105. 107. 107. 106. 105. 102. 94. 96. 95. 97. 103. 91. 91. 93. 93. 94. 84. 88. 89. 88. 80. 80. 75. 74. 77. 80. 84. 88. 94. 93. 99. 94. 95. 96. 100. 97. 101. 110. 115. 117. 113. 109. 113. 108. 110. 109. 86. 87. 91. 89. 86. 71. 71. 69. 61. 62. 63. 65. 65. 68. 78. 81. 83. 85. 90. 92. 92. 88. 92. 93. 95. 99. 102. 105. 105. 106. 107. 112. 98. 97. 101. 95. 97. 92. 94. 106. 107. 107. 110. 116. 119. 121. 119. 120. 117. 116. 112. 113. 108. 105. 105. 104. 106. 99. 92. 94. 98. 99. 96. 97. 96. 104. 104. 94. 93. 81. 83. 84. 87. 83. 82. 83. 88. 90. 92. 96. 99. 98. 101. 107. 97. 105. 110. 105. 108. 104. 104. 102. 104. 106. 104. 93. 97. 92. 92. 79. 74. 73. 80. 76. 70. 74. 73. 74. 77. 71. 73. 75. 76. 81. 85. 92. 80. 84. 75. 72. 72. 74. 76. 79. 85. 88. 90. 94. 98. 101. 101. 105. 106. 101.]
Out[97]:
[<matplotlib.lines.Line2D at 0x19661e0a220>]
Visualize in 3-D¶
In [599]:
rho = 20
population_history = []
events = []
area = space_size[0] * space_size[1]
num_individuals = np.random.poisson(rho * area)
x = np.random.uniform(0, space_size[0], num_individuals)
y = np.random.uniform(0, space_size[1], num_individuals)
population = pd.DataFrame({
"id": np.arange(1, num_individuals + 1),
"x": x,
"y": y,
"parent_id": None,
"ancestor_id": np.arange(1, num_individuals + 1)
})
population_history.append(population.copy())
In [600]:
simulate(population,total_time=5)
In [597]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
class PopulationHistoryVisualizer:
def __init__(self, events, population_history, space_size=(1, 1)):
"""
:param events: 事件列表，每个元素需包含 "time" 字段（时间间隔，非绝对时间）
:param population_history: 种群历史列表，每个元素为DataFrame，包含 "id", "x", "y", "parent_id"
:param space_size: 模拟区域大小 (默认1x1)
"""
self.events = events
self.population_history = population_history
self.space_size = space_size
self.fig = plt.figure(figsize=(6, 6))
self.ax = self.fig.add_subplot(111, projection='3d')
# 计算绝对时间
self.absolute_times = self._calculate_absolute_times()
def _calculate_absolute_times(self):
"""将事件时间间隔转换为绝对时间"""
absolute_times = [0.0] # 初始时间为0
for event in self.events:
absolute_times.append(absolute_times[-1] + event["time"])
return absolute_times
def plot_3d_history(self):
"""三维时空演化可视化"""
# 创建颜色映射
colors = plt.cm.viridis(np.linspace(0, 1, len(self.absolute_times)))
# 绘制每个时间平面
for i, (t, df) in enumerate(zip(self.absolute_times, self.population_history)):
# 绘制时间平面
self._draw_time_plane(t, color=colors[i], alpha=0.05)
# 绘制当前时间点个体
self.ax.scatter(
df["x"], df["y"], t,
c=[colors[i]]*len(df),
s=20, depthshade=False,
label=f't={t:.2f}'
)
# 绘制亲代连线（跳过初始时间）
if i > 0:
self._draw_parent_connections(
prev_df=self.population_history[i-1],
current_df=df,
prev_time=self.absolute_times[i-1],
current_time=t
)
# 设置坐标轴
self.ax.set_xlim(0, self.space_size[0])
self.ax.set_ylim(0, self.space_size[1])
self.ax.set_zlim(0, max(self.absolute_times))
self.ax.set_xlabel('X')
self.ax.set_ylabel('Y')
self.ax.set_zlabel('Time')
plt.legend()
plt.show()
def _draw_time_plane(self, z, color, alpha=0.05):
"""绘制半透明时间平面（逻辑不变）"""
xx, yy = np.meshgrid(
np.linspace(0, self.space_size[0], 2),
np.linspace(0, self.space_size[1], 2)
)
zz = np.full_like(xx, z)
self.ax.plot_surface(
xx, yy, zz,
color=color,
alpha=alpha,
edgecolor='none'
)
def _draw_parent_connections(self, prev_df, current_df, prev_time, current_time):
"""绘制父子代连线（逻辑不变）"""
parent_map = {
row["id"]: (row["x"], row["y"], prev_time)
for _, row in prev_df.iterrows()
}
for _, row in current_df.iterrows():
parent_id = row["parent_id"]
if parent_id in parent_map:
x_prev, y_prev, t_prev = parent_map[parent_id]
x_curr, y_curr = row["x"], row["y"]
self.ax.plot(
[x_prev, x_curr], [y_prev, y_curr], [t_prev, current_time],
color='gray', alpha=1, lw=0.5
)
In [601]:
visualizer = PopulationHistoryVisualizer(events, population_history)
visualizer.plot_3d_history()
Backward Process (not finished yet)¶
In [54]:
import numpy as np
import pandas as pd
theta = 0.1
u0 = 0.8
Lambda = 1
alpha = 1.3
rho = 100
space_size = (1, 1)
population_history = []
events = []
In [55]:
def initialize_backward_population(sample_size=10):
"""初始化采样种群（无父代信息）"""
return pd.DataFrame({
"x": np.random.uniform(0, space_size[0], sample_size),
"y": np.random.uniform(0, space_size[1], sample_size),
"parent_x": None,
"parent_y": None
})
population = initialize_backward_population(sample_size=20)
population_history = [population.copy()]
population
Out[55]:
| x | y | parent_x | parent_y | |
|---|---|---|---|---|
| 0 | 0.868069 | 0.085209 | None | None |
| 1 | 0.364261 | 0.704067 | None | None |
| 2 | 0.563183 | 0.223180 | None | None |
| 3 | 0.173847 | 0.802975 | None | None |
| 4 | 0.847671 | 0.013925 | None | None |
| 5 | 0.831519 | 0.136961 | None | None |
| 6 | 0.134561 | 0.622288 | None | None |
| 7 | 0.899021 | 0.330266 | None | None |
| 8 | 0.552974 | 0.473884 | None | None |
| 9 | 0.103192 | 0.376543 | None | None |
| 10 | 0.612892 | 0.785135 | None | None |
| 11 | 0.012091 | 0.319185 | None | None |
| 12 | 0.148687 | 0.236002 | None | None |
| 13 | 0.191178 | 0.708069 | None | None |
| 14 | 0.591035 | 0.296511 | None | None |
| 15 | 0.156748 | 0.602910 | None | None |
| 16 | 0.832683 | 0.117243 | None | None |
| 17 | 0.686724 | 0.060198 | None | None |
| 18 | 0.115537 | 0.420695 | None | None |
| 19 | 0.676425 | 0.567785 | None | None |
In [56]:
plt.figure(figsize=(4, 4))
plt.xlim(0, space_size[0])
plt.ylim(0, space_size[1])
plt.scatter(population["x"], population["y"],
s=40, c='blue', alpha=0.7,
edgecolors='white', linewidths=0.5,
label='Existing Individuals'
)
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0)
plt.title("initial")
plt.show()
In [58]:
def backward_event_gen(current_time):
rate = Lambda * np.pi * (theta**2) * u0
dt = np.random.exponential(1/rate)
z_x, z_y = np.random.uniform(0, space_size[0]), np.random.uniform(0, space_size[1])
return current_time + dt, (z_x, z_y)
event_time, event_z = backward_event_gen(0)
event_z,event_time
Out[58]:
((0.5082680356856656, 0.43755032261881033), 11.742827337203758)
In [59]:
plt.figure(figsize=(4, 4))
plt.xlim(0, space_size[0])
plt.ylim(0, space_size[1])
plt.scatter(population["x"], population["y"],
s=40, c='blue', alpha=0.7,
edgecolors='white', linewidths=0.5,
label='Existing Individuals'
)
plt.scatter(event_z[0], event_z[1],
s=20, marker='x', c='red',
linewidths=1, zorder=1,
label='Event Center'
)
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0)
plt.title("event")
plt.show
Out[59]:
<function matplotlib.pyplot.show(close=None, block=None)>
In [ ]:
def apply_backward_event(population, event_z):
"""应用单个回溯事件"""
# 计算击中概率
dx = population.x - event_z[0]
dy = population.y - event_z[1]
dist_sq = dx**2 + dy**2
hit_probs = u0 * np.exp(-dist_sq/(2*theta**2))
# 选择被击中的个体
hit_mask = np.random.rand(len(population)) < hit_probs
if not hit_mask.any():
return population
# 生成父代位置（方差θ²(1+α²)）
sigma_parent = theta * np.sqrt(1 + alpha**2)
parent_x = np.random.normal(event_z[0], sigma_parent)
parent_y = np.random.normal(event_z[1], sigma_parent)
if parent_x*parent_y == 0:
return population
# 更新被击中个体的父代信息
population.loc[hit_mask, "parent_x"] = parent_x
population.loc[hit_mask, "parent_y"] = parent_y
# 添加父代个体
new_parents = pd.DataFrame({
"x": [parent_x],
"y": [parent_y],
"parent_x": None,
"parent_y": None
})
population_history.append(population.copy())
# 合并种群：未被击中的个体 + 唯一父代
survivors = population[~hit_mask]
return population[hit_mask], pd.concat([survivors, new_parents], ignore_index=True)
In [73]:
dx = population.x - event_z[0]
dy = population.y - event_z[1]
dist_sq = dx**2 + dy**2
hit_probs = u0 * np.exp(-dist_sq/(2*theta**2))
hit_mask = np.random.rand(len(population)) < hit_probs
hit_pop = population[hit_mask]
sigma_parent = theta * np.sqrt(1 + alpha**2)
parent_x = np.random.normal(event_z[0], sigma_parent)
parent_y = np.random.normal(event_z[1], sigma_parent)
In [74]:
plt.figure(figsize=(4, 4))
plt.xlim(0, space_size[0])
plt.ylim(0, space_size[1])
plt.scatter(population["x"], population["y"],
s=40, c='blue', alpha=0.7,
edgecolors='white', linewidths=0.5,
label='Existing Individuals'
)
plt.scatter(event_z[0], event_z[1],
s=20, marker='x', c='red',
linewidths=1, zorder=1,
label='Event Center'
)
plt.scatter(parent_x, parent_y,
c='red',
edgecolor='white', linewidth=0.5,
label='parent'
)
plt.scatter(hit_pop[0], hit_pop[1],
linestyle=':', c='white', alpha=0.5,
edgecolor='black', linewidth=1.5,
label='death'
)
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0)
plt.title("event")
--------------------------------------------------------------------------- KeyError Traceback (most recent call last) File c:\Users\muliy\.conda\envs\py_transformer\lib\site-packages\pandas\core\indexes\base.py:3361, in Index.get_loc(self, key, method, tolerance) 3360 try: -> 3361 return self._engine.get_loc(casted_key) 3362 except KeyError as err: File c:\Users\muliy\.conda\envs\py_transformer\lib\site-packages\pandas\_libs\index.pyx:76, in pandas._libs.index.IndexEngine.get_loc() File c:\Users\muliy\.conda\envs\py_transformer\lib\site-packages\pandas\_libs\index.pyx:108, in pandas._libs.index.IndexEngine.get_loc() File pandas\_libs\hashtable_class_helper.pxi:5198, in pandas._libs.hashtable.PyObjectHashTable.get_item() File pandas\_libs\hashtable_class_helper.pxi:5206, in pandas._libs.hashtable.PyObjectHashTable.get_item() KeyError: 0 The above exception was the direct cause of the following exception: KeyError Traceback (most recent call last) Cell In[74], line 19 9 plt.scatter(event_z[0], event_z[1], 10 s=20, marker='x', c='red', 11 linewidths=1, zorder=1, 12 label='Event Center' 13 ) 14 plt.scatter(parent_x, parent_y, 15 c='red', 16 edgecolor='white', linewidth=0.5, 17 label='parent' 18 ) ---> 19 plt.scatter(hit_pop[0], hit_pop[1], 20 linestyle=':', c='white', alpha=0.5, 21 edgecolor='black', linewidth=1.5, 22 label='death' 23 ) 24 plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0) 25 plt.title("event") File c:\Users\muliy\.conda\envs\py_transformer\lib\site-packages\pandas\core\frame.py:3458, in DataFrame.__getitem__(self, key) 3456 if self.columns.nlevels > 1: 3457 return self._getitem_multilevel(key) -> 3458 indexer = self.columns.get_loc(key) 3459 if is_integer(indexer): 3460 indexer = [indexer] File c:\Users\muliy\.conda\envs\py_transformer\lib\site-packages\pandas\core\indexes\base.py:3363, in Index.get_loc(self, key, method, tolerance) 3361 return self._engine.get_loc(casted_key) 3362 except KeyError as err: -> 3363 raise KeyError(key) from err 3365 if is_scalar(key) and isna(key) and not self.hasnans: 3366 raise KeyError(key) KeyError: 0
In [974]:
def simulate_backward(population, total_events=100):
"""主模拟函数"""
# 初始化
current_time = 0
for i in range(total_events):
# 生成事件
event_time, event_z = backward_event_gen(current_time)
events.append({"time": event_time, "z_x": event_z[0], "z_y": event_z[1]})
# 应用事件
population = apply_backward_event(population, event_z)
# 终止条件：只剩一个个体
if len(population) == 10:
break
current_time = event_time
return population_history, events
population_history, events = simulate_backward(population, total_events=100)
len(population_history)
Out[974]:
234
Appendix: Random walk¶
random walk OR: Brownian motion
"why Kingman coalescent not works for population distributed in plane space?"
Random walk in 1 dimension¶
In [522]:
import numpy as np
import matplotlib.pyplot as plt
# 参数设定
T = 1000000 # 总时间步数
x = np.zeros(T) # 记录每个时间步的位置信息
step_size = 1 # 每步移动的距离
# 生成布朗运动（随机游走）
for t in range(1, T):
x[t] = x[t-1] + np.random.choice([-step_size, step_size]) # 随机向左或向右移动
# 绘图
plt.figure(figsize=(6, 6))
plt.plot(x, np.arange(T), linestyle='-', markersize=3, color='blue')
plt.axvline(x=0, color='red', linestyle='--', linewidth=1.5, label="x=0 Reference")
plt.xlabel("Position (X)")
plt.ylabel("Time (T)")
plt.title("1D Brownian Motion")
plt.show()
In [516]:
import numpy as np
import matplotlib.pyplot as plt
# 参数设定
T = 100000 # 总时间步数
sigma = 1 # 运动的标准差，控制步长的大小
# 生成布朗运动
x = np.cumsum(np.random.normal(loc=0, scale=sigma, size=T)) # 累积求和，生成连续轨迹
# 绘图
plt.figure(figsize=(6, 6))
plt.plot(x, np.arange(T), linestyle='-', markersize=3, color='blue')
plt.axvline(x=0, color='red', linestyle='--', linewidth=1.5, label="x=0 Reference")
plt.xlabel("Position (X)")
plt.ylabel("Time (T)")
plt.title("1D Brownian Motion with Gaussian Steps")
plt.show()
Random walk in 2 dimension¶
In [539]:
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# 参数设定
T = 10000 # 总时间步数
sigma = 1 # 标准差，控制步长大小
# 生成二维布朗运动轨迹
x = np.cumsum(np.random.normal(loc=0, scale=sigma, size=T)) # X方向步长累积
y = np.cumsum(np.random.normal(loc=0, scale=sigma, size=T)) # Y方向步长累积
t = np.arange(T) # 时间步
In [540]:
# 创建 3D 图形
fig = plt.figure(figsize=(6, 6))
ax = fig.add_subplot(111, projection='3d')
# 绘制布朗运动轨迹
ax.plot(0, 0, 0, marker = 'o', linestyle='-', markersize=5, color='red')
ax.plot(x, y, t, linestyle='-', markersize=3, color='blue', label="2D Brownian Motion")
# 添加 z=0 的参考平面
ax.plot(x, y, np.zeros_like(t), linestyle="--", color='red', alpha=0.5, label="Projection on XY-plane")
# 设置轴标签
ax.set_xlabel("X Position")
ax.set_ylabel("Y Position")
ax.set_zlabel("Time (T)")
ax.set_title("2D Brownian Motion in 3D Space")
# 添加图例
ax.legend()
# 显示图像
plt.show()
In [541]:
# 绘制二维轨迹
plt.figure(figsize=(6, 6))
plt.plot(x, y, linestyle='-', color='blue', alpha=0.7, label="Brownian Motion Path")
# 标记起点和终点
plt.scatter(x[0], y[0], color='green', s=100, edgecolors='black')
# 添加原点参考线
plt.axhline(y=0, color='gray', linestyle='--', alpha=0.5)
plt.axvline(x=0, color='gray', linestyle='--', alpha=0.5)
# 设置标签和标题
plt.xlabel("X Position")
plt.ylabel("Y Position")
plt.title("2D Brownian Motion")
# 添加图例
plt.legend()
# 显示图像
plt.show()
