使用uncertainpy和chaospy的多项式混沌扩展在混沌中寻找秩序

如何使用数学、物理学、Python和数据科学来解决混乱问题

三年前，我从意大利的罗马搬到了美国俄亥俄州辛辛那提市，之后得到了辛辛那提大学的博士学位邀请。我对我的城市有很多的思念：食物、天气、这座永恒之城的美丽之处。我绝对不会怀念的一件事情是我这座城市的疯狂交通。

我的一个好朋友前几天发给我一条短信，说：

“皮耶罗，今天交通太糟糕了，城市完全被混乱包围了。”

现在，显然，我没有纠正他（特别是我知道罗马的交通有多糟糕），但是在数学和物理学中，与我们在日常生活中使用的“混乱”一词相比，“混乱”的意义完全不同。

当我们在数学上提到混乱时，常见的定义是一个问题是由“确定性方程”所调整的，但系统的演变极大地依赖于初始条件。这意味着即使初始条件发生极小变化，系统的演变也可能截然不同。用洛伦兹的话来说，这就是：

“现在决定未来，但近似现在不能近似地决定未来。”

¹ http://mpe.dimacs.rutgers.edu/2013/03/17/chaos-in-an-atmosphere-hanging-on-a-wall/

这意味着我们只能以概率的角度考虑一个状态的演变方式。鉴于一个过程的起点，我们将无法准确预测系统的最终状态，因为它是混乱的，但我们将能够概率地预测它，比如得到平均值和标准偏差。

这种类型的混乱可以用数字方法处理，例如使用Python。在本博客文章中，我们将从抽象的随机游动(Random Walk)开始介绍多项式混沌展开(Polynomial Chaos Expansion)，并在一个真实案例中应用到我们的咖啡温度中☕️

让我们开始吧！

1. 随机游走

对正在阅读的数学家和物理学家来说，随机游走是众所周知的概念。这个模型几乎无所不在，从金融到物理学都有应用，并且非常简单。在文献中它也被称为“布朗运动”，其工作原理如下：

1. 我们从点x = 0开始
2. 以相同的概率，我们可以从x = 0到x = 1或从x = 0到x = -1。我们将这一点定义为x_1
3. 同样，我们可以将x_1的值增加1或减去1。我们将这一点定义为x_2。
4. 我们重复步骤3，使用x_2进行N-2次

有时，我觉得使用伪代码更容易理解，比用文字解释更好

RandomWalk(N):  x = 0  i = 0   while i<N:    p = random(-1,1)    x = x+p    i = i+1  return x

现在让我们来探索一下！

1.2 代码

在这篇文章的这部分，我们将使用Python语言代码来描述随机游走。您需要导入基本的库，如numpy和matplotlib.pyplot。

这是随机游走的代码：

如果我们运行这个程序，假设运行100次，我们将得到以下路径：

非常有趣的是，如果你考虑最后一步，你可以找到高斯分布：

暂时告一段落。我保证我们会用到这个。

2. 微分方程

现在，我们所知道的关于生活的一切，等等，我是说真的一切，都是因为微分方程。

微分方程是物理学用来描述系统演化的工具。我的高中老师曾经这样解释：

“为了描述世界，你需要两样东西：微分和积分。微分很容易，积分却很难。”

例如，让我们考虑松鼠在树上爬的位置y：

假设松鼠的速度为 v(t) = (t/60)**2，其中 t 为秒。所以我们的超级英雄从 v(t=0) = 0 开始，经过两分钟达到速度 v(120) = 2**2 = 4 m/s。

根据这些信息，超级松鼠的位置是什么？

我们需要对速度方程进行积分，得到：

我们如何得到常数 c？我们只需设定在 t = 0 时发生的情况。我们假设松鼠的初始高度为 0，因此 c = 0。

所以松鼠爬树的位置如下：

一般来说，某个确定的解 y 可以视为另一个量 x 的积分，加上一个初始条件。

在这个方案中，我们谈论了：

• 时间（t）是时间变量（从实验开始到结束）
• x 是我们要积分的对象（在上面的例子中，x 是速度）
• 解 y 是通过积分 x 而得到的解（在上面的例子中，y 是超级松鼠的位置）

因此，在这种情况下，我们可以说 x(t) = y(t) 的积分。

但还有更多。你的系统中可能有一些固定的参数，但这些参数可以改变系统的演化。因此：

x(t, 参数列表) = y(t, 参数列表) 的积分

例如，让我们谈谈咖啡。咖啡？是的，咖啡。

2.1 牛顿冷却定律

一个对物理学相当不错（哈哈）的家伙名叫艾萨克·牛顿，在他留给我们的众多礼物中，解释了如何描述热量从一个热体传递过程。换句话说，他告诉我们物体如何冷却下来。

牛顿提出的冷却定律指出，体内热量向外界传递的速率与一个常数k成比例，该常数取决于表面积、热传导系数和时间t时的温度T与环境温度T_env之间的差异。

如果我们想得到温度（T），我们需要对热量的速率（dT/dt）进行积分。这是方程：

要得到温度T，给定T_env和k（记住这个！），我们需要对dT/dt进行积分。

3. 维纳混沌！

很少情况下，我们能够对微分方程定义解析解。这就是我的高中教授说它很难积分的原因。我们更有可能需要进行数值积分，这意味着以数值方式（即用算法）解决微分方程。

对于数值积分，有一些非常知名的方法（算法），比如梯形法或黎曼求和法。它们都能工作，并且具有各自的优缺点，效率也很高。它们不是问题所在。

问题在于微分方程的参数（例如T_env和kappa）。我来解释一下。

还记得上面的方程中的T_env和k吗？我们不知道它们可能是什么，而且它们完全可以改变系统的演变。

诺伯特·维纳美妙的思维给微分方程增加了额外的随机参数。特别是，现在我们的讨论变得有意义了，带有随机参数的微分方程被定义为混沌的，可以使用随机游走（啊哈！）作为多项式来描述。通过这样做，我们能够以概率的方式理解解T(t)。

我知道可能有些困惑：我们一步一步来 ：）

3.1 设置

首先，我们需要在Python中定义我们的微分方程：

如你所见，这不仅涉及到T（我们的变量），还涉及到kappa和T_env。

这是我们需要进行积分的函数。在这之前，让我们先引入一些好朋友🦸‍♂️

你可能会遇到错误，因为你没有chaospy和uncertainpy。它们就是我们的魔法巫师：它们实现了多项式混沌展开方法。安装它们非常简单：

pip install uncertainpy

3.2 积分函数

使用梯形法则来设置要进行积分的函数：

所以：

• 我们设置了咖啡的起始温度，假设为T_0 = 95
• 我们设置了问题的时间步长，假设为500个时间步长
• 我们使用梯形法则来进行积分
• 我们返回时间和温度

3.3 关于uncertainpy

现在，我必须说：uncertainpy太棒了。你可以用它做很多事情，我真的建议你在这里花点时间了解一下。

我们要做的是:

• 我们设定可能的kappa分布。例如，kappa从给定的均值和方差的正态分布中采样。
• 我们对T_env做同样设定。
• 我们应用uncertainpy 并从输入分布中提取可能温度值的分布。

游戏很简单：如果我们知道参数的可能分布，得益于Wiener的混沌，我们就能知道输出的分布。

听起来可能有些混乱，但我发誓，看完代码后会更清楚：

使用coffee_cup定义模型，它是我们的微分方程。然后，我们使用chaospy定义参数分布，并定义相应的参数字典。

现在，对于每个kappa和T_env的值，我们有一个具有不同参数和不同温度T(t)的微分方程，这是积分的结果。由于使用了神奇的chaospy，解决方案变成了具有平均值和标准偏差的分布。

让我们看一下：

就是这样！（就像Biggie）非常简单。

3.4 全部内容

以下是完整的代码块：

这不是很美吗？我们能够将输入参数的分布转化为输出结果的分布。在t=0时，温度是T= T_0 = 95。随着时间的增加，参数的不确定性变得越来越明显。在时间=200分钟时，我们有很大的不确定性（可以是从5到30），具体取决于k和T_env，可能会很冷或稍微有点热。

4. 结果

在这篇博文中，我们介绍了美丽的chaospy和uncertainpy库。这些库使我们能够处理使用随机行走来定义一种多项式混沌形式的Wiener混沌问题。这种多项式混沌用于处理包含分布而非参数的微分方程。我们按照以下顺序进行了处理：

• 在第1章中，我们描述了随机行走。
• 在第2章中，我们描述了微分方程，特别是牛顿冷却定律。
• 在第3章中，我们描述了Wiener的混沌并应用了多项式混沌。

5. 结论

如果你喜欢这篇文章，想了解更多关于机器学习的内容，或者你有什么问题想问我，你可以：

A. 在Linkedin上关注我，我会发布我的故事B. 订阅我的新闻简报，它会及时通知你新的故事，并给你机会发短信给我，以获得你可能有的所有更正或疑问。C. 成为推荐会员，这样你就不会有“每月最大故事数量”的限制了，你可以阅读我（以及其他成千上万的机器学习和数据科学顶级作家）写的关于最新技术的内容。