随机动力系统的Python模拟

### 随机动力系统的Python模拟 #### 1. 随机动力系统简介 随机动力系统是受噪声影响的动力系统，噪声带来的随机性考虑了现实世界现象中的变异性。例如，股票价格的演变通常表现出长期趋势，同时伴随着更快、更小幅度的振荡，反映了每日或每小时的变化。 随机系统在数据科学中的应用包括统计推断方法（如马尔可夫链蒙特卡罗方法）以及时间序列或地理空间数据的随机建模。随机离散时间系统包括离散时间马尔可夫链，其马尔可夫性质意味着系统在时间 $n + 1$ 的状态仅取决于其在时间 $n$ 的状态。随机元胞自动机是元胞自动机的随机扩展，是特殊的马尔可夫链。 对于连续时间系统，带有噪声的常微分方程产生随机微分方程（SDEs），带有噪声的偏微分方程产生随机偏微分方程（SPDEs）。点过程是另一种随机过程，用于模拟瞬时事件在时间（如队列中客户的到达或神经系统中的动作电位）或空间（如森林中树木的位置、领土中的城市或天空中的星星）上的随机发生。 数学上，随机动力系统的理论基于概率论和测度论。连续时间随机系统的研究建立在随机微积分的基础上，随机微积分是将微积分（包括导数和积分）扩展到随机过程。 #### 2. 模拟离散时间马尔可夫链 离散时间马尔可夫链是在状态空间中从一个状态过渡到另一个状态的随机过程，每个时间步都会发生过渡。马尔可夫链的特点是无记忆性，即从当前状态过渡到下一个状态的概率仅取决于当前状态，而不取决于先前的状态。这些模型广泛应用于科学和工程领域。 以下是模拟一个简单的马尔可夫链来建模种群演变的步骤： 1. 导入所需的库： ```python import numpy as np import matplotlib.pyplot as plt %matplotlib inline ``` 2. 定义种群的最大规模、出生率和死亡率： ```python N = 100 # 最大种群规模 a = 0.5/N # 出生率 b = 0.5/N # 死亡率 ``` 3. 模拟在有限空间 $\{0, 1, ..., N\}$ 上的马尔可夫链，每个状态代表一个种群规模。设置初始状态为 $x_0 = 25$： ```python nsteps = 1000 x = np.zeros(nsteps) x[0] = 25 ``` 4. 模拟马尔可夫链： ```python for t in range(nsteps - 1): if 0 < x[t] < N-1: # 是否有出生？ birth = np.random.rand() <= a*x[t] # 是否有死亡？ death = np.random.rand() <= b*x[t] # 更新种群规模 x[t+1] = x[t] + 1*birth - 1*death # 如果达到 0 或 N，演变停止 else: x[t+1] = x[t] ``` 5. 查看种群规模的演变： ```python plt.plot(x) ``` 可以看到，在每个时间步，种群规模可以保持稳定、增加或减少 1。 接下来，模拟该马尔可夫链的多个独立试验： 1. 设置试验次数并初始化种群规模： ```python ntrials = 100 x = np.random.randint(size=ntrials, low=0, high=N) ``` 2. 定义模拟函数： ```python def simulate(x, nsteps): """运行模拟""" for _ in range(nsteps - 1): # 哪些试验需要更新？ upd = (0 < x) & (x < N-1) # 哪些试验发生出生？ birth = 1*(np.random.rand(ntrials) <= a*x) # 哪些试验发生死亡？ death = 1*(np.random.rand(ntrials) <= b*x) # 更新所有试验的种群规模 x[upd] += birth[upd] - death[upd] ``` 3. 查看不同时间的种群规模直方图： ```python bins = np.linspace(0, N, 25) nsteps_list = [10, 1000, 10000] for i, nsteps in enumerate(nsteps_list): plt.subplot(1, len(nsteps_list), i + 1) simulate(x, nsteps) plt.hist(x, bins=bins) plt.xlabel("种群规模") if i == 0: plt.ylabel("直方图") plt.title("{0:d} 时间步".format(nsteps)) ``` 最初，种群规模在 0 到 N 之间看起来是均匀分布的，但经过足够长的时间后，它们似乎收敛到 0 或 N。这是因为状态 0 和 N 是吸收状态，一旦到达，链就不能离开这些状态，并且从任何其他状态都可以到达这些状态。 #### 3. 模拟泊松过程 泊松过程是一种特殊的点过程，是表示瞬时事件随机发生的随机模型。大致来说，泊松过程是结构最少或最随机的点过程，是一种特殊的连续时间马尔可夫过程。点过程，特别是泊松过程，可以模拟随机瞬时事件，如队列或服务器中客户的到达、电话呼叫、放射性衰变、神经细胞的动作电位等。 以下是模拟均匀平稳泊松过程的不同方法： 1. 导入所需的库： ```python import numpy as np import matplotlib.pyplot as plt %matplotlib inline ``` 2. 指定速率值，即每秒的平均事件数： ```python rate = 20. # 每秒的平均事件数 ``` 3. 使用 1 毫秒的小时间间隔模拟过程： ```python dt = .001 # 时间步长 n = int(1./dt) # 时间步数 ``` 4. 在每个时间间隔