1

我需要求解一组延迟微分方程，并且我想在这些方程中使用粉红噪声。

我尝试使用 Python 的 Pydelay 包来实现，但问题是我需要生成噪声然后将其传递给仿真，或者在仿真过程中生成噪声。第一个选项不起作用，因为在此求解器中指定的参数需要保持不变，并且噪声会随着时间而变化。第二个选项不起作用，因为 Pydelay 只支持生成每个样本都独立于其他样本的噪声，而粉红噪声则不是这样。

有人知道怎么做吗？

这是我的代码（我需要 I1 和 I2 是粉红噪声，而不是像我的代码中那样恒定）：

# -*- coding: utf-8 -*-
"""
Created on Tue Aug  2 15:35:15 2016

@author: kasienka
"""
#!python
import pydelay
from time import time
import math
import numpy as np
from numpy import fft
from scipy import integrate
import matplotlib as mpl
mpl.use('Agg')
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.interpolate import griddata
import sys
import pylab as pl
from pydelay import dde23

# define the equations
eqns = {
    'y1' : 'I1 - y1(t-tau) + epsilon * pow(y2(t-tau), 1.1)',
    'y2' : 'I2 - y2(t-tau) + epsilon * pow(y1(t-tau), 1.1)'
    }
#define the parameters
params = {
    'I2' : 0.2,
    'I1' : 0.4,
    'tau': 1.31,
    'epsilon': 0.2
    }

# Initialise the solver
dde = dde23(eqns=eqns, params=params)

# set the simulation parameters
dde.set_sim_params(tfinal=50, dtmax=0.0001)

histdic = {
    'y1': lambda t: 0.2,
    'y2': lambda t: 0.4
    }
dde.hist_from_funcs(histdic, 1000)

# run the simulator
dde.run()
p1 = []
p2 = []
x1 = []
x2 = []
czasy = np.linspace(0, 50, 10000)
for el in times:
    x10 = dde.sol_spl(el)['y1']
    x20 = dde.sol_spl(el)['y2']
    prob1 = 1.0 / (1 + 3.14**(-(x10 - x20)) )
    prob2 = 1.0 / (1 + 3.14**(-(x20 - x10)) )
    p1.append(prob1)
    p2.append(prob2)
    x1.append(x10)
    x2.append(x20)
fig = pl.figure()
pl.plot(times, p1, label = '$p_1$')
pl.plot(times, p2, label = '$p_2$')
pl.xlabel('$time$')
pl.ylabel('$probabilities$')
pl.legend()
pl.savefig(sys.argv[1])
4

2 回答 2

1

我假设您想要动态噪声¹。在这种情况下，您首先应该意识到有几种类型的微分方程具有根本不同的求解器：

  • 延迟微分方程 (DDE)通常使用以下任一方法求解：

    • 对过去进行插值的嵌入式多步 Runge-Kutta 方法（这就是 Pydelay 所做的），

    • 单步积分器，积分步除所有延迟。

  • 随机微分方程 (SDE)使用简单的单步法求解（均基于欧拉法）。多步方法仍然是一个热门话题，嵌入式方法是最近才提出的。

    我所知道的所有关于 SDE 的论文（理论和方法；但诚然它们并不多），只考虑白噪声（维纳过程）；意识到粉红噪声本身就是一个问题。通过快速搜索，我只能找到一篇关于使用白噪声 SDE 模拟白噪声的论文。

  • 随机延迟微分方程 (SDDE)再次需要特殊的求解器。我对它们不是很熟悉，但显然你至少继承了 DDE 和 SDE 的所有限制。因此，最好的求解器将是单步方法，其中积分步骤划分所有延迟。快速浏览一下文献，确实是这样做的。请注意，这并不能解决粉红噪声的问题。

如前所述，Pydelay 使用带插值的多步方法。这不是为随机性而设计的，也无法正确处理随机性（如果是的话，它已经是 SDE 的一个了不起的求解器，即没有延迟）。相反，它所做的是将多步方法视为单步方法，然后添加噪声。正如 Pydelay 的作者自己所说的那样，这“非常粗糙”。此外，实际使用粉红噪声（即使使用粗略的方法）将需要您重新编程软件，并可能导致积分器的误差估计器存在固有问题。我强烈建议不要这样做。即使你成功了，使用高级 DDE 求解器的所有优点（例如，自适应步骤）无论如何都会失去，因此从头开始编写一个新的、更简单的积分器会更容易。

如果您确实需要这样做，我建议您先了解如何解决带有粉红噪声的 SDE，然后将该方法扩展到 SDDE（希望这相当简单）。

¹如果您想要观察噪音，这很简单：只需将其添加到您的解决方案中即可。

于 2017-01-13T13:24:02.727 回答
0

那么问题是为数据集中的所有元素生成粉红噪声？

def voss(nrows, ncols=16):
    """Generates pink noise using the Voss-McCartney algorithm.

    nrows: number of values to generate
    rcols: number of random sources to add

    returns: NumPy array
    """
    array = np.empty((nrows, ncols))
    array.fill(np.nan)
    array[0, :] = np.random.random(ncols)
    array[:, 0] = np.random.random(nrows)

    # the total number of changes is nrows
    n = nrows
    cols = np.random.geometric(0.5, n)
    cols[cols >= ncols] = 0
    rows = np.random.randint(nrows, size=n)
    array[rows, cols] = np.random.random(n)

    df = pd.DataFrame(array)
    df.fillna(method='ffill', axis=0, inplace=True)
    total = df.sum(axis=1)

    return total.values

资料来源：ThinkDSP

于 2017-01-13T12:40:30.933 回答