Python-仿真建模实用指南-二-

布客飞龙
66 阅读1小时+

Python 仿真建模实用指南（二）

原文：Hands-On Simulation Modeling with Python

协议：CC BY-NC-SA 4.0

五、基于仿真的马尔可夫决策过程

马尔可夫决策过程（MDPs）在结果部分随机，部分受决策者控制的情况下，对决策进行建模。MDP 是一个由五个要素构成的随机过程：决策时期、状态、行动、转移概率和奖励。马尔可夫过程的特征要素是系统发现自己的状态以及决策者可以对这些状态执行的可用操作。这些元素确定了两组：可以找到系统的状态集和每个特定状态可用的操作集。决策者选择的行动决定了系统的随机响应，从而使系统进入新状态。这种转变会带来回报，决策者可以用它来评估他们选择的好处。在本章中，我们将学习如何用马尔可夫链处理决策过程。我们将分析马尔可夫过程的基本概念，然后分析一些实际应用，以了解如何为系统的不同状态之间的过渡选择正确的操作。

在本章中，我们将介绍以下主要主题：

  • 马尔可夫过程综述
  • 引入马尔可夫链
  • 马尔可夫链的应用
  • 贝尔曼方程解释
  • 多智能体仿真

技术要求

本章将介绍 MDP。为了处理本章中的主题，您必须具备代数和数学建模的基本知识。

要使用本章中的 Python 代码，您需要以下文件（可在 GitHub 上通过以下 URL 获得：github.com/PacktPublis… ：

  • SimulatingRandomWalk.py
  • WeatherForecasting.py

马尔可夫过程概述

马尔可夫决策过程被定义为一个离散时间随机控制过程。在*第二章*【了解随机性和随机数】中，我们说随机过程是用于根据随机规律仿真系统演化的数值模型。自然现象，无论是自然现象还是观测误差，都具有随机因素的特点。这些因素在系统的观测中引入了一个随机数。这种随机因素决定了观测结果的不确定性，因为不可能确定地预测结果。在这种情况下，我们只能说它将以一定的概率假设许多可能值中的一个。

如果从对系统进行观察的瞬间t开始，则过程的演变将仅取决于t*，而不受先前瞬间的影响。这里，我们可以说随机过程是马尔可夫过程。

重要提示

一个过程被称为马尔可夫，当过程的未来演化仅取决于系统的观察瞬间，而不以任何方式依赖于过去时。

马尔可夫过程的特征要素包括系统发现自己的状态以及决策者可以对这些状态执行的可用操作。这些元素确定了两组：可以找到系统的状态集和每个特定状态可用的操作集。决策者选择的动作决定了系统的随机响应，从而使其进入新状态。这种转变会带来回报，决策者可以用它来评估他们选择的好处。

agent 环境接口

为了实现目标，马尔可夫过程具有两个元素之间相互作用问题的特征。这种相互作用的两个特征元素是主体和环境。agent是必须达到目标的元素，而环境是 agent 必须与之交互的元素。环境对应于代理外部的所有内容。

代理是一个软件，以完全自动和智能的方式执行另一个软件所需的服务。它们被称为智能代理。

代理人的基本特征如下所示：

  • 代理会持续监视环境，此操作会导致环境状态发生变化。
  • 可用操作属于连续或离散的集合。
  • 代理的行为选择取决于环境的状态，这种选择需要一定程度的智能，因为它不是微不足道的。
  • 代理具有所做选择的记忆–智能记忆。

代理的行为以试图实现特定目标为特征。要做到这一点，它会在事先不知道或至少不完全知道的环境中执行操作。这种不确定性通过代理和环境之间的交互作用来填补。在这个阶段，代理通过测量环境来学习了解环境的状态，以这种方式规划其未来的行动。

agent 采用的策略是基于错误理论的原理：证明行为和记忆可能犯的错误，以便反复尝试，直到目标实现。代理不断重复这些操作，导致环境发生变化，从而改变其状态。

重要提示

奖励的概念对代理人未来的选择至关重要，它代表了环境对所采取行动的反应。这一反应与行动在实现目标中所决定的权重成正比：如果行动导致正确的行为，它将是积极的，而如果行动不正确，它将是消极的。

引导代理人实现其目标的决策过程可概括为三个要点：

  • 代理人的目的
  • 与环境的互动
  • 环境的全部或部分不确定性

在此过程中，代理通过传感器进行的测量接收来自环境的刺激。代理根据从环境接收到的刺激来决定要采取什么行动。作为代理人行为的结果，确定环境状态的变化将获得奖励。

决策过程中的关键要素如下图所示：

Figure 5.1 – The agent’s decision-making process

图 5.1——代理人的决策过程

当选择动作时，对环境进行正式描述变得至关重要。此描述必须返回有关环境属性的基本信息，而不是环境的精确表示。

探索 MDP

我们在前一节中讨论过的 agent 与环境的交互，被看作是一个马尔可夫决策过程。这种选择取决于加载问题和计算困难。如马尔可夫过程概述一节所述，马尔可夫决策过程被定义为离散时间随机控制过程。

在这里，我们需要执行一系列操作，每个操作都会导致环境状态的不确定性变化。通过观察环境，我们知道它在执行一个动作后的状态。另一方面，如果无法对环境进行观察，则即使在执行操作后，我们也无法了解状态。在这种情况下，状态是环境所有可能状态的概率分布。在这种情况下，可以将更改过程视为快照序列。

时间t的状态由随机变量st 表示。决策被解释为一个离散时间随机过程。离散时间随机过程是一系列随机变量xt，其中 tN。我们可以定义如下一些元素：

  • 状态空间：随机变量可以假设的值的集合
  • 随机过程的历史（路径）：随机变量序列的实现

环境对某一行为的反应表现为奖励。马尔可夫决策过程中的代理-环境交互可通过下图进行总结：

Figure 5.2 – The agent-environment interaction in MDP

图 5.2–MDP 中的代理环境交互

代理环境交互的基本步骤如上图所示，如下所示：

  1. 代理和环境之间的交互在一段时间内的离散瞬间发生。
  2. 在每一瞬间，代理通过获取其状态stS来监控环境，其中S是可能的状态集。
  3. 代理执行动作aa（st），其中 a（st）是状态st可用的一组可能动作。
  4. 代理根据要实现的目标选择一种可能的行动。
  5. 该选择由策略π（s，a）决定，它表示动作a在s状态下执行的概率。
  6. 在时间t+1 时，代理收到与先前选择的动作对应的数字奖励rt+1r。
  7. 由于选择，环境进入了新的状态。
  8. 代理必须监视环境状态并执行新操作。
  9. 此迭代将重复，直到达到目标。

在我们描述的迭代过程中，st+1 状态取决于之前的状态和所采取的动作。此功能将流程定义为 MDP，可通过以下等式表示：

在前面的等式中，δ表示状态函数。我们可以将 MDP 总结如下：

  1. 代理监控环境的状态并执行一系列操作。

  2. 在离散时间t时，代理检测到当前状态并决定在a执行动作*。*

  3. The environment reacts to this action by returning a reward rt= r (st*, at) and moving to the state st + 1 = δ (st, a*t).

    重要提示

    r 和δ函数是环境的特征，仅取决于当前状态和动作。MDP 的目标是学习一个策略，对于系统的每个状态，该策略为代理提供一个动作，该动作使整个动作序列中累积的总奖励最大化。

现在，让我们分析一下前面介绍的一些术语。它们代表了帮助我们理解马尔可夫过程的关键概念。

报酬函数

奖励函数在马尔可夫过程中识别目标。它将代理检测到的环境状态映射为一个表示奖励的数字。这一过程的目的是使代理人因其选择而获得的长期总回报最大化。然后，奖励函数收集从代理选择的操作中获得的积极和消极结果，并使用它们修改策略。如果根据策略提供的指示选择的操作返回低奖励，则将修改策略以选择其他操作。报酬函数具有两个功能：它刺激决策的效率，并决定代理人的风险规避程度。

政策

策略决定代理在决策方面的行为。它映射了环境的状态和在这些状态下要选择的动作，代表了一组对刺激做出反应的规则或关联。策略是马尔可夫代理的基本组成部分，因为它决定了其行为。在马尔可夫决策模型中，策略提供了一个解决方案，该解决方案将推荐的操作与代理可能实现的每个状态相关联。如果策略在可能的操作中提供了最高的预期效用，则称为最优策略（π*）。这样，代理就不必将以前的选择保存在内存中。要做出决策，代理只需要执行与当前状态关联的策略。

状态值函数

状态值函数为我们提供了评估代理状态质量所需的信息。它返回根据每个状态的策略获得的预期目标的值，该值由总预期奖励表示。代理依赖于策略来选择要执行的操作。

了解折扣累积奖励

MDP 的目标是学习一个策略，该策略指导代理为环境的每个状态选择要执行的动作。该政策旨在最大限度地提高代理人在整个行动过程中获得的总奖励。让我们学习如何最大化这个总回报。通过采用政策获得的总回报计算如下：

在前面的等式中，rT 是将环境带入终端状态sT 的行为的回报。

为了获得最大的总回报，我们可以选择为每个州提供最高回报的行动，从而选择使总回报最大化的最优策略。

重要提示

此解决方案不适用于所有情况；例如，当目标或终端状态未在有限的步骤数内实现时。在这种情况下，rt 和您想要最大化的奖励之和趋于无穷大。

另一种方法是使用折扣累积奖励，尝试最大化以下金额：

在前面的式中，γ被称为折扣因子，代表未来奖励的重要性。折扣系数是 0≤ γ≤ 1 且具有以下条件：

  • γ<1：序列rt收敛为有限值。
  • γT=γ= 0 席 T1：代理不考虑未来的奖励，从而试图最大限度地奖励目前的状态。
  • γ=1：代理人倾向于未来奖励而非即时奖励。

在学习过程中，折扣系数的值可能会有所不同，以考虑特殊的操作或状态。最佳策略可能包括回报较低的个人行动，前提是总回报较高。

勘探开发概念比较

当达到目标时，代理寻找最有回报的行为。要做到这一点，他们必须将每个动作与返回的奖励联系起来。在具有多个状态的复杂环境中，由于存在多个动作-奖励对，这种方法不可行。

重要提示

这就是众所周知的探索-开发困境：对于每一个国家，代理人探索所有可能的行动，利用在实现目标中获得最大回报的行动。

决策需要在两种可用方法之间进行选择：

  • 开发：根据当前信息做出最佳决策
  • 勘探：最好的决策是收集更多的信息

最佳的长期战略可能会带来短期牺牲，因为这种方法需要收集足够的信息以做出最佳决策。

在日常生活中，我们经常发现自己不得不在两种选择中做出选择，至少在理论上，这两种选择会导致相同的结果：这种方法就是探索-开发困境。例如，假设我们需要决定是选择我们已经知道的（开发）还是选择新的（探索）。开发使我们的知识保持不变，而探索使我们更多地了解系统。很明显，探索使我们面临错误选择的风险。

让我们看一个在现实生活中使用这种方法的例子——我们必须选择到达我们信任的餐厅的最佳路径：

  • 利用：选择你已经知道的路径。
  • 探索：尝试一条新路。

在复杂问题中，收敛到最优策略可能太慢。在这些情况下，这个问题的解决办法表现为勘探和开发之间的平衡。

完全基于探索的代理在每个状态下都会随机行为，收敛到实际上不可能的最优策略。相反，如果一个代理完全基于剥削进行操作，它们将始终使用相同的操作，这可能不是最优的。

引入马尔可夫链

马尔可夫链是离散动态系统，表现出可归因于马尔可夫过程的特征。这些是有限状态系统——有限马尔可夫链——其中从一种状态到另一种状态的转换是基于概率而非确定性的。在通用时刻t的链的可用信息由其处于任何状态的概率提供，并且通过指定这些概率如何在时刻t+1从时刻t更新来指定链的时间演化。

重要提示

马尔可夫链是一种随机模型，在这种模型中，系统随着时间的推移以这样一种方式演化：过去只通过现在影响未来：马尔可夫链对过去没有记忆。

给出了一个随机过程，其特征是一系列随机变量X=X0，Xn，其值在一组j0*、j1、n 中。如果该过程的演化仅取决于当前状态，即n*步之后的状态，则该过程是马尔可夫过程。使用条件概率，我们可以用以下等式表示该过程：

如果离散时间随机过程X具有马尔可夫性质，则称为马尔可夫链。如果以下转移概率不依赖于n，且仅依赖于i和j，则称马尔可夫链是齐次的：

在这些假设中，假设我们有以下几点：

所有联合概率都可以通过知道数字 pij 和以下初始分布来计算：

该概率表示零时间过程的分布。概率 pij 称为转移概率，pij 是每个时间阶段从i到j的转移概率。

转移矩阵

通过采用矩阵表示，齐次马尔可夫链的应用很容易。通过这一点，前面的方程所表达的公式变得更具可读性。我们可以通过以下转移矩阵来表示马尔可夫链的结构：

这是一个正矩阵，其中每行元素之和为酉。事实上，i第行的元素是在瞬间 t 处于状态 Si 的链通过S1或 S2的概率。或者下一个瞬间的 Sn。这种转变是相互排斥的，并且是所有可能的。这种单位和线的正矩阵是随机的。我们将每个向量正直线 x 称为随机线，使得 T=[x1 x2…xn]，其中元素之和假设为单位值：

通过执行单个实验，转移矩阵具有从结果 i 传递到结果 j 的位置（i，j）。

过渡图

转移矩阵不是描述马尔可夫链的唯一解。另一种选择是称为过渡图的定向图，其中顶点由状态S1*、S2、*n 标记，当且仅当从 Si 到 Sj 的转移概率为正时，存在一条连接顶点 Si 和顶点 Sj 的直边。

重要提示

转移矩阵和转移图包含表示相同马尔可夫链所需的相同信息。

让我们来看看一个例子：考虑一个具有三种可能状态的马尔可夫链——1, 2 和 3——由下面的转移矩阵表示：

如前所述，转换矩阵包含与转换图相同的信息。让我们学习如何绘制这个图表。有三种可能的状态–1、2 和 3–从每个状态到其他状态的直接边界显示了转变概率 pij。当状态i到状态j没有箭头时，表示 pij=0：

Figure 5.3 – Diagram of the transition matrix

图 5.3–过渡矩阵图

在前面的转换图中，状态中的箭头加起来总是正好为 1，就像转换矩阵中的每一行一样。

马尔可夫链的应用

现在，让我们来看一系列使用马尔可夫链可以实现的实际应用。我们将介绍这个问题，然后分析 Python 代码，让我们仿真它的工作方式。

引入随机游动

随机游动识别一类用于仿真由一系列随机步骤组成的路径的数学模型。模型的复杂性取决于我们想要仿真的系统特征，这些特征由自由度和方向表示。这个词的作者是卡尔·皮尔森，他在 1905 年第一次提到了“随意散步”这个词。在该模型中，每个步骤都有一个随机方向，该方向通过一个随机过程演变，该过程涉及遵循精确统计分布的已知量。随时间追踪的路径不一定是真实运动的描述：它只是返回变量随时间的演化。这就是该模型在所有科学领域广泛使用的原因：化学、物理、生物学、经济学、计算机科学和社会学。

一维随机游动

一维随意行走仿真了一个准时粒子的运动，该粒子必然沿着直线运动，因此只有两个运动：右运动和左运动。每一个移动都与一步的随机移动相关联，向右移动的概率为p，向左移动的概率为q。每一步的长度相同，并且独立于其他步骤。下图显示了准时粒子绑定到的路径，以及允许的方向和两个顶点：

Figure 5.4 – One-dimensional walk

图 5.4–一维步行

n 次通过后，该点的位置将由其横坐标X（n）确定，其特征为随机项。我们的目标是计算在n步之后粒子返回起点的概率。

重要提示

该点将实际返回其起始位置的想法并不确定。为了表示点在直线上的位置，我们将采用 X（n）变量，该变量表示粒子移动 n 步后直线的横坐标：该变量是一个具有二项式分布的离散随机变量。

点粒子的路径可以概括为：对于每个瞬间，粒子根据随机变量*Z**（n）*返回的值向右或向左移动一步。此随机变量仅取两个二分法值：

  • +1，概率 p>0
  • -概率为 q 的 1

这两种概率通过以下等式相互关联：

让我们考虑随机变量 Po.t0，z，Po.t1，n，与 Po.t2，n，t3，t＝1, 2，…假设这些变量是独立的且分布均匀。瞬间n的粒子位置将由以下等式表示：

在前面的公式中，Xn 是行走中的下一个值，Xn-1 是前一个时间段的观察值，Zn 是该步骤中的随机波动。

重要提示

Xn 变量标识马尔可夫链；也就是说，下一时刻粒子处于某个位置的概率仅取决于当前位置，即使我们知道当前位置之前的所有时刻。

仿真一维随机游动

一次随意行走的仿真并不代表一个随机数的小序列，因为当前一步的下一步代表了它的进化。接下来两个步骤之间的依赖性保证了从一个通道到下一个通道的某种一致性。这在一个平庸的独立随机数生成中是不能保证的，相反，它会返回一个数字到另一个数字之间的巨大差异。让我们学习如何通过以下伪代码在简单的随意行走模型中表示要执行的动作序列：

  1. 从 0 位置开始。
  2. 随机选择一个二分法值（-1，1）。
  3. 将此值添加到上一时间步。
  4. 继续重复步骤 2。

这个简单的迭代过程可以在 Python 中通过处理随机行走的 1000 个时间步列表来实现。让我们来看一看：

  1. Let's start by loading the necessary libraries:

    from random import seed
from random import random
from matplotlib import pyplot

    random模块实现各种分布的伪随机数生成器。随机模块基于 Mersenne Twister 算法。Mersenne Twister 是一种伪随机数生成器。最初开发用于生成蒙特卡罗仿真的输入，通过 Mersenne Twister 生成几乎一致的数字，使其适用于广泛的应用。

    从random模块导入两个库：seed和random。在这段代码中，我们将生成随机数。要做到这一点，我们将使用random()函数，它在每次调用时生成不同的值。任何数字重复都有很长的一段时间。这对于生成唯一的值或变体非常有用，但有时，以不同的方式处理相同的数据集非常有用。这对于确保实验的再现性是必要的。为此，我们可以使用seed库中包含的random.seed()函数。此函数用于初始化基本随机数生成器。

    matplotlib库是用于打印高质量图形的 Python 库。使用matplotlib可以通过几个命令生成图形、直方图、条形图、功率谱、误差图、散点图等。这是一组命令行函数，类似于 MATLAB 软件提供的那些函数。

  2. Now, we will investigate the individual operations. Let's start with setting the seed:

    seed(1)

    如果我们希望以不同的方式处理相同的数据集，random.seed()功能非常有用，因为这可以使仿真重现。

    重要提示

    此函数用于初始化基本随机数生成器。如果在两个连续的仿真中使用相同的种子，则始终会得到相同的数字对序列。

  3. Let's move on and initialize the crucial variable of the code:

    RWPath= list()

    RWPath变量表示一个列表，其中包含代表随机游动的值序列。列表是值的有序集合，可以是各种类型。它是一个可编辑的容器，意味着我们可以添加、删除和修改现有值。出于我们的目的，当我们希望通过路径的后续步骤不断更新我们的值时，列表表示最合适的解决方案。list()函数接受一系列值并将其转换为列表。使用前面的命令，我们只是初始化了当前为空的列表，并使用以下代码开始填充它：

    RWPath.append(-1 if random() < 0.5 else 1)

    我们添加到列表中的第一个值是二分法值。这只是决定要添加的值是1还是-1的问题。然而，选择是随机的。在这里，我们使用random()函数生成一个介于 0 和 1 之间的随机数，然后检查它是否为<0.5。如果是，则添加-1；否则，我们添加1。此时，我们将使用一个for循环的迭代循环，该循环将重复 1000 个步骤：

    for i in range(1, 1000):

    在每个步骤中，我们将生成一个随机项，如下所示：

    	ZNValue = -1 if random() < 0.5 else 1

    重要提示

    正如我们选择要添加到列表中的第一个值时所做的那样，我们使用random()函数生成一个随机值，因此如果返回的值低于0.5，则ZNValue变量采用-1值；否则，1。

  4. Now, we can calculate the value of the random walk at the current step:

    XNValue = RWPath[i-1] + ZNValue

    XNValue变量表示当前步长上横坐标的值。它由两个术语组成：第一个术语表示前一状态中横坐标的值，第二个术语表示生成随机值的结果。必须将此值添加到列表中：

    RWPath.append(XNValue)

    对于我们要执行的 1000 个步骤，将重复此过程。在循环结束时，我们将把整个序列存储在列表中。

  5. Finally, we can visualize it through the following piece of code:

    pyplot.plot(RWPath)
pyplot.show()

    pyplot.plot()函数使用x作为索引数组，在y轴上绘制RWPath列表中包含的值，其值如下：0..N-1。plot()函数用途极其广泛，可以接受任意数量的参数。

    最后，pyplot.show()函数显示创建的图形，如下所示：

Figure 5.5 – The trend plot of the random walk path

图 5.5–随机行走路径的趋势图

在前面的图中，我们可以分析随机过程中点粒子跟随的路径。该曲线可以描述一般功能的趋势，不一定与道路路线相关。正如预期的那样，该过程被配置为马尔可夫过程，因为下一步独立于上一步的位置，并且仅取决于当前步骤。散散步是一种广泛应用于金融领域的数学模型。事实上，它被广泛用于仿真来自市场的信息的效率：价格随着新信息的到来而变化，这与我们已经知道的无关。

仿真天气预报

马尔可夫链的另一个潜在应用是开发天气预报模型。让我们学习如何在 Python 中实现该算法。首先，我们可以用一个简化的模型进行工作：我们只考虑两个气候条件/状态，即晴天和雨天。我们的模型将假设明天的天气条件将受到今天天气条件的影响，使该过程具有马尔可夫特征。这两个状态之间的联系将由以下过渡矩阵表示：

转移矩阵返回条件概率 P（A | B），表示事件 A 发生在事件 B 发生之后的概率。因此，该矩阵包含以下条件概率：

在上一个转移矩阵中，每行包含一个完整的分布。因此，所有数字必须是非负的，且总和必须等于 1。气候条件显示出抵抗变化的趋势。因此，在一个晴天之后，另一个晴天–*P（晴天）*的概率大于一个雨天–*P（晴天）*的概率。明天的气候条件与昨天的气候条件没有直接关系；因此，这个过程是马尔可夫过程。前面的转移矩阵等价于以下内容：

Figure 5.6 – Transition diagram

图 5.6–过渡图

我们要详细说明的仿真模型必须计算未来几天下雨的概率。它还必须允许您恢复一段时间内晴天和雨天比例的统计数据。如前所述，该过程是马尔可夫过程，我们在前面章节中分析的工具允许我们获得所需的信息。让我们开始：

  1. Let's see the Python code that alternates sunny and rainy days, starting from a specific initial condition. As always, we will analyze it line by line, starting with loading the necessary libraries:

    import numpy as np
import matplotlib.pyplot as plt

    numpy库是一个 Python 库，它包含许多函数，可以帮助我们管理多维矩阵。此外，它还包含大量高级数学函数，我们可以在这些矩阵上使用这些函数。我们将使用两个函数：random.seed()和random.choose()。

    matplotlib库是用于打印高质量图形的 Python 库。使用matplotlib可以通过几个命令生成图形、直方图、条形图、功率谱、误差图、散点图等。它是一组命令行函数，类似于 MATLAB 软件提供的那些函数。让我们继续说明代码：

    np.random.seed(3)

    random.seed()函数初始化随机数生成器的种子。这样，使用随机数的仿真将是可复制的。实验的再现性将通过以下事实得到保证：将生成的随机数始终相同。

  2. Now, let's define the possible states of the weather conditions:

    StatesData = ['Sunny','Rainy']

    提供两种状态：晴天和雨天。表示天气条件之间转换的转换矩阵将设置如下：

    TransitionStates = [['SuSu','SuRa'],['RaRa','RaSu']]
TransitionMatrix = [[0.80,0.20],[0.25,0.75]]

    转移矩阵返回条件概率P（A|B），表示事件B发生后事件A发生的概率。一行中的所有数字必须是非负的，且总和必须等于 1。让我们继续并设置将包含状态转换列表的变量：

    WeatherForecasting = list()

    WeatherForecasting变量将包含天气预报的结果。此变量将为list类型。

    重要提示

    列表是值的有序集合，可以是各种类型。它是一个可编辑的容器，允许我们添加、删除和修改现有值。

    对于我们的目的，即通过路径的后续步骤不断更新我们的值，列表表示最合适的解决方案。list()函数接受一系列值并将其转换为列表。

  3. 现在，我们决定预测天气状况的天数：

    NumDays = 365

  4. 目前，我们已决定仿真 1 年时间范围内的天气预报；即 365 天。让我们修复一个包含当天预测的变量：

    TodayPrediction = StatesData[0]

  5. 此外，我们还使用包含可能状态的第一个向量值初始化它。该值对应于Sunny条件。我们在屏幕上打印此值：

    print('Weather initial condition =',TodayPrediction)

  6. 此时，我们可以预测NumDays变量设定的每一天的天气条件。为此，我们将使用一个for循环，该循环将执行同一段代码，执行次数等于我们预先设置的天数：

    for i in range(1, NumDays):

  7. 现在，我们将分析整个程序的主要部分。在for循环中，通过一个附加的条件结构if语句来预测连续每一天的时间。从TodayPrediction变量中包含的气象条件开始，我们必须预测第二天的天气。我们有两种情况：晴天和下雨。实际上有两种控制条件，如下代码所示：

    if TodayPrediction == 'Sunny':        
        TransCondition = np.random.choice(TransitionStates[0],replace=True,p=TransitionMatrix[0])
        if TransCondition == 'SuSu':
            pass
        else:
            TodayPrediction = 'Rainy'

 elif TodayPrediction == 'Rainy':
        TransCondition = np.random.choice(TransitionStates[1],replace=True,p=TransitionMatrix[1])
        if TransCondition == 'RaRa':
            pass
        else:
            TodayPrediction = 'Sunny'

  8. If the current state is Sunny, we use the numpy random.choice() function to forecast the weather condition for the next state. A common use for random number generators is to select a random element from a sequence of enumerated values, even if these values are not numbers. The random.choice() function returns a random element of the non-empty sequence passed as an argument. Three arguments are passed:

    TransitionStates[0]：第一行过渡态

    replace=True：该样品为替换品

    p=TransitionMatrix[0]：与通过状态下的每个条目相关联的概率

    random.choise()函数根据TransitionStates矩阵中包含的值返回SuSu、SuRa、RaRa和RaSu类型的随机样本。前两个将从晴天开始返回，其余两个将从雨天开始返回。这些值将存储在TransCondition变量中。

    在每个if语句中，还有一个额外的if语句。这用于确定是更新天气预报的当前值还是保持不变。让我们看看如何：

    if TransCondition == 'SuSu':
            pass
        else:
            TodayPrediction = 'Rainy'

    如果TransCondition变量包含SuSu值，则当天的天气条件保持不变。否则，将其替换为Rainy值。elif条款从降雨条件开始执行类似程序。在for循环的每次迭代结束时，更新天气预报列表，并打印当前预报：

    WeatherForecasting.append(TodayPrediction) 
print(TodayPrediction)

    现在，我们需要预测未来 365 天的天气预报。

  9. Let's draw a graph with the sequence of forecasts for the next 365 days:

    plt.plot(WeatherForecasting)
plt.show()

    打印以下图表：

    Figure 5.7 – Plot of the weather forecast

    图 5.7–天气预报图

    在这里，我们可以看到晴天的天气预报比雨天的天气预报更为普遍。

    重要提示

    顶部的平坦点代表所有的晴天，而中间的凹陷点代表雨天。

  10. To quantify this prevalence, we can draw a histogram. In this way, we will be able to count the occurrences of each condition:

```py

plt.figure()
plt.hist(WeatherForecasting)
plt.show()
```

以下是未来 365 天天气状况的柱状图：

Figure 5.8 – Histogram of the weather forecast

图 5.8–天气预报直方图

有了这个，我们可以确认晴天的普遍性。我们得到的结果来自于转移矩阵。事实上，我们可以看到，太阳状况持续的概率大于降雨的概率。此外，初始条件已设置为阳光条件。我们还可以尝试看看当初始条件设置为降雨条件时会发生什么。

解释了贝尔曼方程

1953 年，理查德·贝尔曼（Richard Bellman）引入了动态规划原理，以有效地解决顺序决策问题。在这类问题中，决策会定期执行，并影响模型的大小。反过来，这些因素会影响未来的决策。Bellman 阐述的优化原则允许通过智能应用程序有效处理决策和模型大小之间交互的复杂性。动态规划技术也从一开始就应用于没有时间或顺序方面的问题。

重要提示

虽然通过提供一个通用的抽象模型，动态规划可以应用于广泛的问题，但从实用的角度来看，许多问题需要这样的维度的模型来排除任何计算方法。这种不便后来被称为“维度诅咒”，并且是对计算复杂性概念的一种预期，用仍然非正式的术语来说。

动态规划的最大成功在于序列决策模型，特别是随机类型的决策模型，如马尔可夫决策过程，但也包括一些组合模型。

动态规划概念

动态规划（DP是一种规划技术，旨在基于以MDP形式的环境完美模型计算最优策略。动态规划的基础是使用状态值和行动值来确定好的政策。

DP方法应用于马尔可夫决策过程，使用两个相互作用的过程，称为政策评估和政策改进：

  • 政策评估是通过迭代过程来完成的，该过程旨在求解贝尔曼方程。k 过程的收敛性→ ∞ 施加近似规则，从而引入停止条件。
  • 策略改进基于当前值改进策略。

在策略迭代技术中，刚才描述的两个阶段是交替的，每个阶段在另一个阶段开始之前结束。

重要提示

评估策略时的迭代过程要求我们通过一个迭代过程在每一步评估策略，该迭代过程的收敛性是未知的，取决于起始策略。为了解决这个问题，我们可以在某个时刻停止评估策略，同时仍然确保我们收敛到最佳值。

最优性原则

动态优化程序的有效性由 Bellman 的最优性原则保证：最优策略的性质是，无论初始状态和初始决策是什么，剩余决策必须构成与第一个决策产生的状态相关的最优策略。

基于这一原则，可以使用动态确定的目标函数值，将问题分为多个阶段，并按顺序解决这些阶段，而不考虑导致这些阶段的决策。这允许我们一次优化一个阶段，将初始问题减少为一系列更小的子问题，因此更容易解决。

贝尔曼方程

贝尔曼方程通过寻找最优策略和价值函数，帮助我们求解 MDP。最优值函数 V*（S）是返回状态最大值的函数。该最大值是对应于在每个状态下使最佳行动的奖励值最大化的行动的最大值。然后，通过递归过程，将贴现因子乘以下一状态的值，再乘以贝尔曼方程。以下是贝尔曼方程的一个示例：

在上一个等式中，我们有以下等式：

  • 为s状态下的功能值。
  • 是我们在s状态下扮演一个角色后得到的奖励。
  • γ是贴现系数。
  • 为下一状态的功能值。

对于一个随机系统，当我们采取行动时，并不是说我们将在以后的状态中结束，而是我们只能指出在那个状态中结束的概率。

多智能体仿真

agent 可以定义为任何能够通过传感器感知环境并通过执行器在其中起作用的事物。人工智能关注的是理性代理的概念，或者总是试图优化适当性能度量的代理。rational agent 可以是人类 agent、机器人 agent 或软件 agent。在下图中，我们可以看到代理和环境之间的交互：

Figure 5.9 – Interaction between the agent and the environment

图 5.9–代理和环境之间的交互

当 agent 能够灵活、独立地选择要采取的行动以实现其目标，而无需不断地求助于外部决策系统的干预时，它被认为是自治的。请注意，在大多数复杂域中，代理只能部分获取信息并在其插入的环境中拥有控制权，因此最多只能对其施加一定的影响。

如果代理具有以下特征，则可以将其视为自治和智能代理：

  • 反应性：它必须感知环境，设法及时适应所发生的变化。
  • 主动性：必须表现出以主动性为目标的行为。
  • 社交技能：它必须能够与其他代理人互动，以实现他们的目标。

在许多情况下，多个代理共存于同一环境中，并以不同的方式相互作用。事实上，代理表示一个孤立的系统是非常罕见的。我们可以将多智能体系统（MAS定义为一组可能相互交互的智能体。MAS 可以是竞争性的，即每个代理都试图完全最大化自己的利益，甚至牺牲其他代理的利益，而不是合作性的，即代理愿意放弃部分目标，试图最大化系统的全球效用。

可能的交互类型如下所示：

  • 协商：当代理必须就分配给某些变量的值寻求一致时，会发生。
  • 合作：这发生在有共同目标的情况下，代理人试图调整和协调他们的行动。
  • 协调：这是一种旨在避免代理之间冲突的交互。

MAS 系统的使用带来了一系列优势：

  • 效率和速度：由于可以并行执行计算。
  • 鲁棒性：系统能够克服单 agent 故障。
  • 灵活性：向系统中添加新代理非常简单。
  • 模块化：在软件设计阶段非常有用，因为可以重用代码。
  • 成本：与整个系统相比，单个单元代理的成本非常低。

基于决策过程处理问题的多智能体系统越来越受到重视，这与它们的一些区别特征有关，例如灵活性和通过相互交互的不同计算单元表示独立实体的可能性。事实上，决策系统的各个利益相关者可以建模为自治代理。一些实际应用程序最近采用了一种基于问题的方法，例如满足和优化分布式约束，以及识别在集中（最优）和非协调（坏）约束之间具有中间效率的规则。

总结

在本章中，我们学习了马尔可夫过程的基本概念。在这一点上，过程的未来演变只取决于对系统的即时观察，而决不取决于过去。我们已经看到了一个代理和周围环境是如何相互作用的，以及描述其行为特征的元素。我们现在了解了决策背后的奖励和政策概念。然后，我们继续通过分析控制其演化的矩阵和转移图来探索马尔可夫链。

然后，我们讨论了一些应用程序，以便将我们学到的概念付诸实践。我们采用一种基于马尔可夫链的方法来处理一次随意散步和天气状况的预测模型。接下来，我们研究了 Bellman 方程作为最优值函数的一致性条件来确定最优策略。最后，我们引入了多智能体系统，它允许我们在决策过程中考虑不同的利益相关者。

在下一章中，我们将了解如何获得总体参数的置信区间和标准误差的稳健估计，以及如何估计统计量的失真和标准误差。然后，我们将了解如何执行统计显著性测试以及如何验证预测模型。*

六、重采样方法

重采样方法是随机仿真和随机数最有趣的推理应用之一。它们在非参数领域特别有用，因为传统的推理方法无法正确应用。它们生成随机数，分配给随机变量或随机样本。它们需要与重复操作增长相关的机器时间。它们的实现非常简单，一旦实现，它们就会自动执行。选择所需元素必须提供代表总体的样本，或者至少可以提供代表总体的样本。为了实现这一点，样本中必须包含人口的所有特征。在本章中，我们将尝试推断从整个人口的代表性样本中获得的结果。鉴于此推断可能出错，有必要评估样本的准确度和得出错误预测的风险。在本章中，我们将学习如何应用重采样方法来近似样本分布的某些特征，以便验证统计模型。我们将分析最常见的重采样方法的基础知识，并通过解决一些实际案例来学习如何使用它们。

在本章中，我们将介绍以下主要主题：

  • 介绍重采样方法
  • jackknife 技术初探
  • 解谜引导
  • 解释置换测试
  • 交叉验证技术探讨

技术要求

在本章中，我们将讨论重采样方法技术。为了处理本章中的主题，您必须具备代数和数学建模的基本知识。要使用本章中的 Python 代码，您需要以下文件（可在 GitHub 上通过以下 URL 获得：github.com/PacktPublis… ：

  • JackknifeEstimator.py
  • BootstrapEstimator.py
  • KfoldCrossValidation.py

介绍重采样方法

重采样方法是一套基于使用数据子集的技术，可以随机提取，也可以根据系统程序提取。这项技术的目的是近似样本分布的某些特征——统计、测试或估计器——以验证统计模型。

重采样方法是随机仿真和随机数生成最有趣的推理应用之一。这些方法在 20 世纪 60 年代开始流行，起源于蒙特卡罗方法的基本概念。蒙特卡罗方法的发展主要发生在 20 世纪 80 年代，随着信息技术的进步和计算机能力的提高。在经典推理方法无法正确应用的情况下，它们的有用性与非参数方法的发展有关。

从重采样方法可以观察到以下细节：

  • 他们多次重复简单的操作。
  • 它们生成随机数，分配给随机变量或随机样本。
  • 随着重复操作数量的增加，它们需要更多的机器时间。
  • 它们的实现非常简单，一旦实现，它们就会自动执行。

随着时间的推移，各种重采样方法已经开发出来，并且可以根据某些特征进行分类。

重要提示

可以在基于随机抽取样本数据子集的方法和根据非随机程序进行重采样的方法之间进行第一分类。

可按如下方式进行进一步分类：

  • bootstrap 方法及其变体（如子抽样）属于随机抽取类别。
  • 刀切和交叉验证等程序属于非随机类别。
  • 统计测试，称为置换或精确测试，也包括在重采样方法家族中。

抽样概念概述

抽样是所有统计研究的基础课题之一。采样生成一组基本单元，即总体的子集，其属性与整个总体相同，至少具有定义的错误风险。

所谓总体，我们指的是所有基本单位的集合，有限的或无限的，这些基本单位具有一定的特征，这表明它们是同质的。

重要提示

例如，这可以是每个地方的温度值总体，时间跨度可以是每天、每月或每年。

抽样理论是统计推断的一个完整的和准备部分，以及由此产生的抽样技术，允许我们识别要分析变量的单位。

统计抽样是一种用于随机选择项目的方法，以使总体中的每个项目都有已知的非零概率包含在样本中。随机选择被认为是构建代表性样本的有力手段，该样本的结构和多样性反映了所考虑的群体。统计抽样使我们能够获得一个客观的样本：元素的选择不取决于出于研究便利性或可用性的原因而定义的标准，也不系统地排除和支持群体中的任何元素组。

在与概率计算相关的随机抽样中，执行以下操作：

  • 通过数学公式外推结果并估计相关误差
  • 控制得出与现实相反结论的风险
  • 通过公式计算获得给定精度和精密度所需的最小样本量

关于抽样的推理

现在，让我们来了解为什么可能更倾向于分析样本的数据而不是整个人群的数据：

  • Consider a case in which the statistical units do not present variability. Here, it is useless to make many measurements because the population parameters are determined with few measurements. For example, if we wanted to determine the average of 1,000 identical statistical units, this value would be equal to that obtained if we only considered 10 units.

    重要提示

    如果并非所有人口要素都可用，则使用抽样。例如，对过去的调查只能根据现有的历史数据进行，而这些数据往往是不完整的。

  • 在获得结果时，如果节省了大量时间，则表示采样。这是因为，即使使用电子计算机，如果调查仅限于总体人口的几个要素，数据输入阶段也会大大缩短。

抽样的利弊

收集信息后，对构成研究人群的所有单位进行调查。当对收集到的信息进行分析时，可以只对构成人口的部分单位进行分析。

抽样的优点如下：

  • 降低成本
  • 缩短时间
  • 减轻组织负担

抽样的缺点如下：

  • 抽样基数并不总是可用或容易知道

在参考群体的组成或大小部分未知的情况下，可通过强制选择进行抽样。抽样并不总能取代完整的调查，例如关于婚姻状况、出生和死亡的调查：所有个案都必须知道。

概率抽样

在概率抽样中，必须提取人口的每个单位的概率是已知的。相比之下，在非概率抽样中，必须提取人口的每个单位的概率是未知的。

让我们来看一个例子。如果我们通过抽签的方式从任何一天在校的大学生中抽取样本，我们无法获得概率样本，原因如下：未入学的学生没有进入的机会，而入学人数最多的学生比随后几年的其他学生更有可能被抽取。

抽样工作原理

抽样程序涉及一系列需要适当遵循的步骤，以便提取能够充分代表总体的数据。抽样工作如下：

  1. 在检测统计中定义目标总体。
  2. 定义采样单位。
  3. 确定样本的大小。
  4. 根据抽样方法，选择一个或多个样本，对其上的载荷进行统计检测。
  5. 最后，对样本的优劣做出判断。

探索刀切技术

此方法用于估计统计数据的失真和标准偏差等特征。这种技术使我们能够获得所需的估计，而无需借助参数假设。Jackknife 是基于计算我们获得的子样本的感兴趣的统计数据，一次忽略一个样本观察。刀切估计对于各种样本统计是一致的，例如均值、方差、相关系数、最大似然估计等。

定义刀切法

Jackknife 方法由 M.H.Quenouille 于 1949 年提出，由于当时的计算能力较低，他创建了一种需要固定数量账户的算法。

重要提示

该方法的主要思想是每次从原始样本中截取不同的观测值，并重新评估感兴趣的参数。该估算值将与在原始样本上计算的同一估算值进行比较。

由于变量的分布未知，估计量的分布也未知。

每次在原始样品外保留一个观察值xi，由此构建刀切样品，如下式所示：

然后，获得尺寸为m=n-1 的n样本。让我们来看一个例子。考虑尺寸 n＝5 的样品，其产生尺寸为 m＝4 的五刀切样品，如下：

伪值在通用ith 样本折刀上重新计算。该程序在每个可用的刀切样品上重复n次：

下图显示了该初步程序：

图 6.1——折刀法的表示

为了计算刀切估算的方差，将使用以下等式：

在前面的等式中，术语定义如下：

计算出的标准偏差将用于建立参数的置信区间。

为了评估并可能减少估计器失真，失真的刀切估计值计算如下：

本质上，刀切法可以减少偏差，并评估估计器的方差。

估计变异系数

为了对不同分布之间的变异性进行比较，我们可以使用变异系数（CV，因为它考虑了分布的平均值。变异系数是色散的相对度量，是一个无量纲量值。它允许我们评估平均值周围值的分散度，而不考虑测量单位。

重要提示

例如，以美元表示的收入样本的标准差与以欧元表示的相同收入的标准差完全不同，而在这两种情况下，分散系数是相同的。

使用以下方程式计算变异系数：

在前面的等式中，我们使用以下参数：

  • 是分布的标准偏差。
  • 是分布平均值的绝对值。

方差是一组数据中每个观测值与数据算术平均值之间的差值平方的平均值：

因此，它代表了我们所犯的平方误差，平均而言，用平均的平均值 T0（T1）代替一般的观测席。标准偏差是方差的平方根，因此代表均方误差的平方根：

CV 可以从平均值和标准偏差开始定义，是比较两个性状变异性的合适指标。当您想要比较具有不同测量单位或不同变化范围的数据分散度时，CV 特别有用。

使用 Python 应用刀切重采样

现在，让我们看一看一些 Python 代码，将分布的 CV 与根据 Jackknife 方法通过重采样获得的 CV 进行比较：

  1. Let's see the code step by step, starting with loading the necessary libraries:

    import random
import statistics 
import matplotlib.pyplot as plt

    random模块实现各种分布的伪随机数生成器。random模块基于 Mersenne Twister 算法。Mersenne Twister 是一种伪随机数生成器。最初开发用于生成蒙特卡罗仿真的输入，通过 Mersenne Twister 生成几乎一致的数字，使其适用于广泛的应用。

    statistics模块包含许多从数值数据计算数理统计的功能。利用本模块中可用的工具，可以计算平均值，并测量中心位置和扩散测量值。

    matplotlib库是用于打印高质量图形的 Python 库。使用matplotlib可以通过几个命令生成图形、直方图、条形图、功率谱、误差图、散点图等。这是一组命令行函数，类似于 MATLAB 软件提供的那些函数。

  2. We now generate a distribution that represents our data population. We will use this data to extract samples using the sampling methods we are studying. To do this, we first create an empty list that will contain such data:

    PopData = list()

    列表是值的有序集合，可以是各种类型。它是一个可编辑的容器–事实上，它允许我们添加、删除和修改现有值。为了不断更新我们的价值观，该列表代表了最合适的解决方案。list()函数接受一系列值并将其转换为列表。使用此命令，我们只需初始化当前为空的列表。

  3. The list will be populated through the generation of random numbers. Then, to make the experiment reproducible, we will fix the seed in advance:

    random.seed(5)

    如果我们想要以不同的方式处理相同的数据集，那么random.seed ()函数非常有用，因为这使得仿真具有可再现性。

    重要提示

    此函数用于初始化基本随机数生成器。如果在两个连续的仿真中使用相同的种子，则始终会得到相同的数字对序列。

  4. Now, we can populate the list with 100 randomly generated values:

    for i in range(100):
    DataElem = 10 * random.random()
    PopData.append(DataElem)

    在前一段代码中，我们使用random()函数生成了 100 个介于 0 和 1 之间的随机数。然后，对于for循环的每个步骤，将该数字乘以 10，以获得 0 到 10 之间的数字分布。

  5. Now, let's define a function that calculates the coefficient of variation, as follows:

    def CVCalc(Dat):
    CVCalc = statistics.stdev(Dat)/statistics.mean(Dat) 
    return CVCalc

    如估算变异系数一节所述，该系数只是标准偏差和平均值之间的比率。为了计算标准偏差，我们使用了statistics.stdev()函数。此函数计算样本标准偏差，它表示样本方差的平方根。为了计算数据的平均值，我们使用了statistics.mean函数。此函数用于计算数据的样本算术平均值。我们可以立即使用新创建的函数来计算我们创建的分布的变异系数：

    CVPopData = CVCalc(PopData)
print(CVPopData)

    返回以下结果：

    0.6569398125747403

    现在，我们省略了这个结果，但我们稍后将使用它来比较我们通过重采样获得的结果。

  6. Now, we can move on and resample according to the Jackknife method. To begin, we fix the variables that we will need in the following calculations:

    N = len(PopData)
JackVal = list()
PseudoVal = list()

    N表示起始分布中存在的样本数。JackVal列表将包含Jackknife样本，PseudoVal列表将包含Jackknife伪值。

  7. The two newly created lists must be initialized to zero to avoid problems in subsequent calculations:

    for i in range(N-1):
    JackVal.append(0)
for i in range(N):
    PseudoVal.append(0)

    JackVal列表的长度为N-1，与我们在定义刀切法一节中讨论的内容相关。

  8. At this point, we have all the tools necessary to apply the Jackknife method. We will use two for loops to extract the samples from the initial distribution by calculating the pseudo value at each step of the external loop:

    for i in range(N):
    for j in range(N):
        if(j < i): 
            JackVal[j] = PopData[j]
        else:
            if(j > i):
                JackVal[j-1]= PopData[j]
    PseudoVal[i] = N*CVCalc(PopData)-
                       (N-1)*CVCalc(JackVal)

    Jackknife samples（Apple T0）是通过在外部环路的每个步骤留下原始样品的观察席（为 To.T2。在外部循环的每个步骤结束时，使用以下方程式计算伪值：

  9. To analyze the distribution of pseudo values, we can draw a histogram:

    plt.hist(PseudoVal)
plt.show()

    打印以下图形：

    图 6.2–伪值分布

  10. Now, let's calculate the average of the pseudo values that we have obtained:

```py

MeanPseudoVal=statistics.mean(PseudoVal)
print(MeanPseudoVal)
```

返回以下结果：

```py
0.6545985339842991
```

正如我们所见，我们得到的值与我们从初始分布得到的值是可比的。现在，我们将计算伪值的方差：

```py
VariancePseudoVal=statistics.variance(PseudoVal)
print(VariancePseudoVal)
```

返回以下结果：

```py
0.2435929299444099
```

最后，让我们评估刀切估计器的方差：

```py
VarJack = statistics.variance(PseudoVal)/N
print(VarJack)
```

返回以下结果：

```py
0.002435929299444099
```

我们将使用这些结果来比较不同的重采样方法。

揭开自举的神秘面纱

最著名的重采样技术是 B.Efron 在 1993 年引入的一种定义为自举的技术。bootstrap 方法的逻辑是构建未观察到的样本，但在统计上与观察到的样本相似。这是通过提取过程对观测序列重新采样实现的，我们在提取过程中重新插入观测值。

引入自举

这个过程就像从 urn 中提取一个数字，然后在下一次提取之前重新插入该数字。一旦选择了统计检验，则对观察到的样本和大量与观察到的样本大小相同的样本进行计算，并通过重新取样获得。检验统计量的N值允许我们定义样本分布；即所选统计数据的经验分布。

重要提示

统计检验是一种区分样本的规则，如果观察到样本，会导致拒绝初始假设，而如果观察到样本，则会导致接受相同假设，直到得到证明。

由于自举样本源于随机抽取过程，并从原始序列重新整合，因此观察序列的任何时间相关结构都不会保留。因此，自举样本具有观察到的样本等属性，但至少大致尊重独立性假设。这使得它们适用于计算测试统计分布，假设没有趋势、变化点或通用系统时间趋势的零假设。

一旦知道零假设下一般检验统计量的样本分布，就可以将观察样本上计算的统计量值与从样本分布中推导出的分位数进行比较，并检查该值是否落入显著性水平为 5%和 10%的关键区域，分别地或者，您可以定义在观察样本上计算的统计值被来自 N 个样本的值超过的时间百分比。该值是观察样本的统计 p 值，并检查该百分比是否远离通常采用的 5%和 10%的含义。

引导定义问题

Bootstrap 是一种统计重采样技术，它通过再入来近似统计的样本分布。因此，它允许我们近似估计值的均值和方差，以便我们可以建立置信区间，并在不知道感兴趣的统计分布时计算测试 p 值。

重要提示

Bootstrap 是基于这样一个事实：唯一可用的样本用于生成更多的样本，并构建理论参考分布。使用原始样本的数据计算统计数据并估计其样本分布，而无需对分布模型进行任何假设。

使用插件原理生成分发；也就是说，θ的估计是通过替换未知总体分布函数的经验等价物来构建的。样本的分布函数是通过构造其在该实验情况下可以假设的所有值的频率分布来获得的。

在简单随机抽样的简单情况下，操作如下。考虑一个观察到的样品，用如下公式描述：

从该分布中，对 m 个数量等于n、例如x1、xm的其他恒定样本进行重新采样。在每次引导提取中，可以多次提取样本第一个元素的数据。提供的每一个都有 1/n 的被提取概率。

假设 E 是我们感兴趣研究的θ的估计量，例如，E（x）=θ。该数量是针对每个引导样本 E（x1），…，E（xm）计算的。通过这种方式，θ的 m 估计可用，由此可以计算 bootstrap 平均值、bootstrap 方差、bootstrap 百分位数等。这些值是相应未知值的近似值，并携带有关 E（x）分布的信息。因此，从这些估计量开始，可以计算置信区间、检验假设等。

使用 Python 进行引导重采样

我们以类似于 Jackknife重采样的方式进行。我们将生成一个随机分布，根据 bootstrap 方法进行重采样，然后比较结果。让我们一步一步地查看代码，以了解过程：

  1. Let's start by importing the necessary libraries:

    import random
import numpy as np 
import matplotlib.pyplot as plt

    NumPy 是一个 Python 库，它包含许多函数，可以帮助我们管理多维矩阵。此外，它包含大量高级数学函数，我们可以在这些矩阵上使用这些函数。

  2. Now, we will generate a distribution that represents our data population. We will use this data to start extracting samples using the sampling methods we have studied. To do this, we will create an empty list that will contain such data:

    PopData = list()

    此列表将通过生成随机数来填充。

  3. To make the experiment reproducible, we will fix the seed in advance:

    random.seed(7)

    如果我们想要以不同的方式处理相同的数据集，那么random.seed``()函数非常有用，因为它可以使仿真重现。

  4. Now, we can populate the list with 1,000 randomly generated values:

    for i in range(1000):
    DataElem = 50 * random.random()
    PopData.append(DataElem)

    在前一段代码中，我们使用函数random()生成了 1000 个介于 0 和 1 之间的随机数。然后，对于for循环的每个步骤，将该数字乘以 50，以获得 0 到 50 之间的数字分布。

  5. At this point, we can start extracting a sample of the initial population. The first sample can be extracted using the random.choices() function, as follows:

    PopSample = random.choices(PopData, k=100)

    此函数从具有替换的总体中提取大小为k元素的样本。我们从 1000 个元素的原始总体中提取了 100 个元素的样本。

  6. Now, we can apply the bootstrap method, as follows:

    PopSampleMean = list()
for i in range(10000):  
    SampleI = random.choices(PopData, k=100)
    PopSampleMean.append(np.mean(SampleI))

    在这段代码中，我们创建了一个包含示例的新列表。在这里，我们使用了 1000 步的for循环。在每个步骤中，使用random.choices ()函数从初始总体中提取 100 个元素的样本。然后，我们得到了这个样本的平均值。然后将该值添加到列表的末尾。

    重要提示

    我们使用替换项对数据进行了重采样，从而使重采样大小与原始数据集的大小相等。

  7. We now print a histogram of the sample we obtained to visualize its distribution:

    plt.hist(PopSampleMean)
plt.show()

    打印以下图表：

    图 6.3–样本分布直方图

    在这里，我们可以看到样本具有正态分布。

  8. We can now calculate the mean of the three distributions that we have generated. Let's start with the bootstrap estimator:

    MeanPopSampleMean = np.mean(PopSampleMean)
print("The mean of the Bootstrap estimator is ",MeanPopSampleMean)

    返回以下结果：

    The mean of the Bootstrap estimator is  24.105354873028915

    然后，我们可以计算初始总体的平均值：

    MeanPopData = np.mean(PopData)
print("The mean of the population is ",MeanPopData)

    返回以下结果：

    The mean of the population is  24.087053989747968

    最后，我们计算从初始总体中提取的简单样本的平均值：

    MeanPopSample = np.mean(PopSample)
print("The mean of the simple random sample is ",MeanPopSample)

    返回以下结果：

    The mean of the simple random sample is  23.140472976536497

我们现在可以比较结果了。在这里，总体和自举样本均值实际上是相同的，而一般样本均值偏离这些值。这告诉我们，与从中提取的一般样本相比，引导样本更能代表初始总体。

比较 Jackknife 和 bootstrap

在本节中，我们将通过强调两种研究的抽样方法的优缺点来比较它们：

  • 引导需要大约 10 倍的计算工作量。至少从理论上讲，折刀可以手工完成。
  • 自举在概念上比折刀简单。Jackknife 需要对 n 个样本进行n次重复，而 bootstrap 需要一定数量的重复。这导致选择一个数字来使用，这并不总是一件容易的事情。一般的经验法则是这个数字是 1000，除非你有足够的计算能力。
  • 引导引入了错误，这是由于完成了重新采样而产生了额外的变化源。请注意，对于大尺寸或仅使用特定引导样本集的情况，此错误会减少。
  • Jackknife 比 bootstrap 更保守，因为它会产生稍大的估计标准误差。
  • 由于复制品之间的微小差异，Jackknife 始终提供相同的结果。另一方面，引导在每次运行时提供不同的结果。
  • Jackknife最适合于估计配对一致性度量的置信度区间。
  • Bootstrap 对于扭曲的分布执行得更好。
  • Jackknife 最适合原始数据的小样本。

解释排列测试

当观察一个属于一组可能结果的现象时，我们会问自己概率定律是什么，我们可以分配给这个集合。统计检验提供了一条规则，允许我们根据样本观察结果决定是否拒绝假设。

参数方法对于实验计划和总体模型非常不确定。如果不遵守这些假设，尤其是当数据法则不符合测试的需要时，参数结果的可靠性就会降低。如果假设不是基于数据分布的知识，并且假设没有得到验证，则使用非参数检验。非参数检验提供了一个非常重要的选择，因为它们需要更少的假设。

置换测试是随机测试的一个特例，它使用从统计推断中得出的一系列随机数。现代计算机的计算能力使其广泛应用成为可能。这些方法不要求满足其关于数据分布的假设。

通过以下步骤执行置换测试：

  1. 定义的统计值与所研究过程或关系的强度成正比。
  2. 定义了一个零假设H0。
  3. 数据集是基于对实际观测数据的置乱而创建的。根据零假设定义混合模式。
  4. 重新计算参考统计数据，并将该值与观察到的值进行比较。
  5. 最后两个步骤重复多次。
  6. 如果观察到的统计数据大于 95%的基于随机洗牌的情况下获得的限值，则拒绝H0。

两个实验使用相同样本空间中的值，分别分布在P1 和P2 下，这两个分布都是未知总体分布的成员。假设同一数据集为x，如果使用相同的检验统计量获得的x上的推理条件相同，则假设在无效假设中满足各组的可交换性。排列测试的重要性在于其健壮性和灵活性。使用这些方法的想法是从数据中生成参考分布，并参考产生的离散规律重新计算数据每个排列的测试统计数据。

接近交叉验证技术

交叉验证是一种基于预测准确性原则的模型选择程序中使用的方法。样本分为两个子集，其中第一个子集（训练集）用于构建和估计，而第二个子集（验证集）用于验证估计模型预测的准确性。通过对重复预测的综合，获得了模型精度的度量。交叉验证方法就像折刀，每次只留下一个观察结果。在另一种称为 K-fold 验证的方法中，样本被划分为 K 个子集，而每个子集又被作为一个验证集忽略。

重要提示

交叉验证可用于估计统计学习技术的均方误差（MSE）测试（或一般而言，任何精度度量），以评估其性能或选择其灵活性水平。

交叉验证可用于回归和分类问题。仿真模型的三种主要验证技术是验证集方法、遗漏交叉验证（LOOCV）和 k-折叠交叉验证方法。在以下部分中，我们将更详细地了解这些概念。

验证集方法

此技术包括将可用数据集随机分成两部分：

  • 训练器材
  • 一个验证集，称为保持集

统计学习模型适用于训练数据，并随后用于使用验证集的数据进行预测。

对产生的测试误差的测量，通常是回归情况下的 MSE，提供了实际测试误差的估计值。事实上，验证集是抽样程序的结果，因此不同的抽样会导致对测试误差的不同估计。

这种验证技术有各种优点和缺点。让我们来看看几个：

  • 该方法具有较高的变异性；也就是说，结果可以随着所选测试集的变化而发生实质性变化。
  • 只有一部分可用单位用于功能估算。这可能导致函数估计的精度降低，并高估测试误差。

LOOCV和 k-fold 交叉验证技术试图克服这些问题。

在交叉验证中遗漏一项

LOOCV 还将观察集分为两部分。但是，我们不创建两个大小相当的子集，而是执行以下操作：

  1. 使用单个观察（x1、y1）进行验证，剩余观察组成训练集。

  2. 该函数基于训练集的n-1观察值进行估计。

  3. The prediction is made using x1. Since (x1, y1) was not used in the function estimate, an estimate of the test error is as follows:

    但是，即使MSE1 对测试误差是公平的，这也是一个很差的估计，因为它是非常可变的。这是因为它是基于单个观察（x1、y1）。

  4. The procedure is repeated by selecting for validation (x2, y2), where a new estimate of the function is made based on the remaining n-1 observations, and calculating the test error again, as follows:

  5. 重复此方法 n 次会产生n测试错误。

  6. The LOOCV estimate for the MSE test is the average of the n MSEs available, as follows:

与验证集方法相比，LOOCV 具有一些优势：

  • 使用n-1单元估计函数的偏差较小。因此，LOOCV 方法不会高估测试误差。
  • 由于测试集的选择不存在随机性，因此相同初始数据集的结果不存在可变性。

LOOCV 可能是计算密集型的，因此对于大型数据集，计算需要很长时间。然而，在线性回归的情况下，存在计算强度低的直接计算公式。

K-折叠交叉验证

在k-折叠交叉验证（k-折叠 CV中），观察集被随机分为大小大致相同的k组或文件夹。将第一个文件夹视为验证集，并对其余 k-1 文件夹的功能进行估计。然后根据所保留文件夹的观测值计算均方误差MSE1。该程序重复 k 次，每次选择不同的文件夹进行验证，从而获得测试误差的k估计值。通过平均这些值计算 k 倍 CV 估计值，如下所示：

如果 k<< n. Furthermore, the k-fold CV tends to have less variability than the LOOCV on different size datasets n，该方法具有计算量小的优点。

选择k是 k-折叠交叉验证的关键。当交叉验证中k发生变化时会发生什么？让我们看看k的极端选择意味着什么：

  • 较高的k值会导致更大的训练集，从而减少偏差。这意味着验证集较小，因此差异较大。
  • 较低的k值会导致较小的训练集，从而导致更大的偏差。这意味着验证集更大，因此方差更低。

使用 Python 进行交叉验证

在本节中，我们将看一个应用交叉验证的示例。首先，我们将创建一个示例数据集，其中包含要识别的简单数据，以便验证算法正在执行的过程。然后，我们将应用 k-折叠交叉验证并分析结果：

  1. As always, we start by importing the necessary libraries:

    import numpy as np
from sklearn.model_selection import KFold

    numpy是一个 Python 库，其中包含许多函数，可帮助我们管理多维矩阵：此外，它还包含大量可用于这些矩阵的高级数学函数。

    Scikit learn 是一个开源 Python 库，它为机器学习提供了多种工具。特别是，它包含许多分类、回归和聚类算法；支持向量机；逻辑回归；还有更多。自 2007 年发布以来，Scikit learn 已成为机器学习领域中使用最多的库之一，无论是有监督的还是无监督的，这要归功于它提供的工具的广泛性，也要归功于它的 API，它是有文档记录的、易于使用的和多功能的。

    重要提示

    应用程序编程接口（API）是应用软件创建和集成的一组定义和协议。它们允许产品或服务与其他产品或服务进行通信，而不知道它们是如何实现的，从而简化了应用程序开发，并实现了时间和金钱的净节约。在创建新工具和产品或管理现有工具和产品时，API 提供了灵活性；简化设计、管理和使用；提供创新机会。

    scikit 学习 API将功能用户界面与众多分类和元分类算法的优化实现结合在一起。它还提供了各种各样的数据预处理、交叉验证、优化和模型评估功能。Scikit learn 在学术研究中特别受欢迎，因为开发人员只需更改几行代码，就可以使用工具来试验不同的算法。

  2. Now, let's generate the starting dataset:

    StartedData=np.arange(10,110,10)
print(StartedData)

    这里，我们生成了一个包含 10 个整数的向量，从值 10 到 100，步长等于 10。为此，我们使用了 NumPyarange()函数。此函数在一定范围内生成等距值。已通过三个参数，如下所示：

    10：间隔的开始。包括该值。如果省略此值，则使用默认值 0。

    110：范围结束。除了浮点数的情况外，此值不包括在范围内。

    10：值之间的间距。这是两个相邻值之间的距离。默认情况下，此值等于 1。

    返回了以下数组：

    [ 10  20  30  40  50  60  70  80  90 100]

  3. Now, we can set the function that will allow us to perform k-fold cross validation:

    kfold = KFold(5, True, 1)

    Sklearn 的KFold()函数通过将数据集划分为 k 个连续折叠来执行 k 折叠交叉验证，默认情况下不进行洗牌。然后，每个折叠被用作一次验证，而剩余的 k-1 折叠形成训练集。通过了三个参数，如下所示：

    5：需要的折叠次数。此数字必须至少为 2。

    True：可选布尔值。如果等于True，则将数据混合后再分批。

    1：随机数生成器使用的种子。

  4. Finally, we can resample the data by using k-fold cross validation:

    for TrainData, TestData in kfold.split(StartedData):
	print("Train Data :", StartedData[TrainData],"Test Data :", StartedData[TestData])

    为此，我们对kfold.split()方法生成的元素使用了一个循环，返回数据集划分的索引。然后，对于每个步骤（等于折叠次数），打印绘制的子集的元素。

    返回以下结果：

    Train Data : [ 10  20  40  50  60  70  80  90] 
Test Data : [ 30 100]
Train Data : [ 10  20  30  40  60  80  90 100] 
Test Data : [50 70]
Train Data : [ 20  30  50  60  70  80  90 100] 
Test Data : [10 40]
Train Data : [ 10  30  40  50  60  70  90 100] 
Test Data : [20 80]
Train Data : [ 10  20  30  40  50  70  80 100] 
Test Data : [60 90]

    这些数据对（Train Data、Test Data将依次用于训练模型和验证模型。这样，您可以避免过度拟合和偏差问题。每次评估模型时，都会使用数据集的提取部分，数据集的其余部分用于训练。

总结

在本章中，我们学习了如何对数据集重新采样。我们分析了几种不同的解决问题的技术。首先，我们分析了抽样的基本概念，了解了促使我们使用从人群中提取的样本的原因。然后我们分析了这一选择的利弊。我们还分析了重采样算法的工作原理。

然后我们解决了第一种重采样方法：折刀法。我们首先定义了该方法背后的概念，然后进入过程，该过程允许我们从原始总体中获取样本。为了将我们学到的概念付诸实践，我们将 Jackknife 重采样应用于一个实际案例。

然后，我们探索了 bootstrap 方法，该方法与观察样本一样，构建了未观察到但具有统计意义的样本。这是通过提取过程对观测序列进行重新采样来实现的，在提取过程中，我们重新插入观测值。在定义了方法之后，我们通过一个示例来强调过程的特征。此外，还对 Jackknife 和 bootstrap 进行了比较。

在分析了置换测试的基本概念之后，我们通过查看各种交叉验证方法来总结本章。通过一个例子加深了我们对 k-折叠交叉验证方法的了解。

在下一章中，我们将学习各种优化技术的基本概念以及如何实现它们。我们将了解数值优化技术和随机优化技术之间的区别，并将学习如何实现随机梯度下降。然后，我们将发现如何估计缺失或潜在变量并优化模型参数。最后，我们将了解如何在实际应用中使用优化方法。

七、利用仿真改进和优化系统

仿真模型允许我们使用很少的资源获得大量的信息。正如生活中经常发生的那样，仿真模型也需要进行改进，以提高其性能。通过优化技术，我们试图修改模型的性能，以在结果方面和尝试利用资源时获得改进。优化问题通常非常复杂，无法通过解析方法确定解。复杂性首先由定义问题大小的变量和约束的数量决定，然后由非线性函数的可能存在决定。为了解决优化问题，有必要使用迭代算法，在给定当前近似解的情况下，通过适当的操作序列确定对该近似的更新。从初始近似开始，确定逐步改进解决方案的近似序列。

在本章中，我们将学习如何使用主要的优化技术来提高仿真模型的性能。我们将学习如何使用梯度下降技术、牛顿-拉斐逊方法和随机梯度下降。我们还将通过实际例子学习如何应用这些技术。

在本章中，我们将介绍以下主要主题：

  • 介绍数值优化技术
  • 探索梯度下降技术
  • 面对 Newton-Raphson 方法
  • 加深对随机梯度下降的认识
  • 在 Python 中发现多元优化应用程序

技术要求

在本章中，我们将学习如何使用仿真模型来改进和优化系统。要处理本章中的主题，您必须具备代数和数学建模的基本知识。要使用本章中的 Python 代码，您需要以下文件（可在 GitHub 上通过以下 URL 获得：github.com/PacktPublis… ：

  • gradient_descent.py
  • newton_raphson.py
  • scipy_optimize.py

引入数值优化技术

在现实生活中，优化意味着在多个可用选项中选择最佳选项。我们每个人都会优化行程以到达目的地，安排我们的一天，如何使用储蓄，等等。在数学中，优化是指确定函数变量的值，使其假定为最小值或最大值。优化是一门处理应用中有用模型的公式化的学科，因此使用有效的方法来确定最佳解决方案。

优化模型对于提供的许多应用具有极大的实际意义。事实上，有许多决策过程需要您确定最小化成本或最大化收益的选择，因此可归因于优化模型。在最优化理论中，数学最优化模型占据着一个相关的位置，对于数学最优化模型，评价函数和描述允许方案的约束通过方程和不等式表示。数学优化模型有不同的形式：

  • 线性优化
  • 整数优化
  • 非线性优化

定义一个优化问题

优化问题包括尝试确定属于集合F的点，其中函数F取尽可能低的值。此问题以以下形式表示：

在这里，我们有以下几点：

  • f称为目标函数

  • F is called the feasible set and contains all the admissible choices for x

    重要提示

    如果你有一个最大化问题，也就是说，如果你必须找到一个函数f具有最高可能值的点，你总是可以回到最小问题，从而改变目标函数的符号。

通过满足上述关系使函数f最小化的元素称为全局最优解，也称为最优解或最小解。目标函数f的对应值称为全局最优值，也称为最优或最小值。

优化问题的复杂性，即其求解的难度，显然取决于目标函数的特性和柔性集的结构。通常，优化问题的特征是向量x的选择是否完全自由。因此，我们可以说有两类问题，如下所示：

  • 无约束极小化问题，若F=Rn；即柔性组F与整组Rn 重合
  • 约束极小化问题，若FRn；也就是说，如果柔性集合F只是集合Rn 的一部分

在解决优化问题时，我们面临的第一个困难是理解值是否正确放置，在这个意义上，函数 F（x）可能不存在一个点 F，其中函数 F（x）取具有最小十进制值的 pi 值。事实上，可能发生以下情况之一：

  • 柔性套件F可能为空。
  • 柔性集合F可能不为空，但目标函数的下限可以等于−∞.
  • 柔性集合F不能为空，目标函数的下限可以等于−∞ 但是，同样在这种情况下，在f上不可能有f的全局最低点。

优化问题中存在全局最小点的一个充分但非必要的条件是 Weierstrass 定理通过以下命题表示的条件：设FRn 为非空紧集。设f为f上的连续函数。如果是，则在f中存在f的全局最低点。

前面的命题只适用于柔性集是紧集的约束问题。为了建立非紧柔性集问题的存在性结果，即在F=Rn 的情况下，有必要尝试刻画包含问题最优解的F的子集。

重要提示

紧空间是拓扑空间，它的每个开覆盖都包含一个完备子覆盖。

一般来说，手头的问题并不总是有一个最优的解决方案，而且，如果存在，也不总是唯一的。

解释局部最优

不幸的是，所有全局最优性条件的应用兴趣有限。事实上，它们与容许集上目标函数的整体行为有关，因此，从计算角度来看，它们必然由复杂条件描述。除了通过定义优化模型引入的全局最优性概念之外，还应该定义局部最优性的概念。

我们可以将局部最优定义为问题在可能解的小邻域内的最佳解。在下图中，我们可以确定函数*f（x）*的四个局部最小条件，因此表示局部最优：

Figure 7.1 – Local minimum conditions for the f(x) function

图 7.1–f（x）函数的局部最小条件

然而，其中只有一个是全局最优，而其他三个仍然是局部最优。从应用的角度来看，局部最优性条件更有用。这些只是必要条件，但总的来说，它们是不够的。这是因为指定点是最小化问题的局部最小点。因此，从理论上看，它们并没有给出优化问题的局部极小值的令人满意的表征，但它们在定义最小化算法中起着重要的作用。

现实生活中面临的许多问题都可以表示为非线性优化问题。从技术和科学的角度来看，这激发了人们对研究和开发解决这类数学难题的方法的兴趣。

探索梯度下降技术

任何仿真算法的目标都是减少模型预测值与数据返回的实际值之间的差异。这是因为实际值和期望值之间的误差较小，表明该算法的仿真效果良好。减少这种差异仅仅意味着最小化所建立模型所基于的目标函数。

确定下降方式

下降法是迭代法，从初始点 x0Rn 开始，生成一系列点{xn}nn，由以下等式定义：

这里，向量是一个搜索方向，标量是一个称为步长的正参数，表示我们在方向上移动的距离。

在下降法中，选择向量和参数，以确保目标f函数在每次迭代时减少，如下所示：

使用向量，我们取下降方向，即线x=xn*+与梯度向量f（xn）形成钝角。这样，只要足够小，就可以保证f*的减少。

根据的选择，有不同的下降方法。最常见的情况如下：

  • 梯度下降法
  • 牛顿-拉斐逊法

让我们从分析梯度下降算法开始。

逼近梯度下降算法

梯度是一个向量值函数，表示函数图切线的斜率，指示函数最大增长率的方向。让我们考虑下面图中所示的凸函数：

Figure 7.2 – The convex function

图 7.2–凸函数

梯度下降算法的目标是达到该函数的最低点。在更专业的术语中，梯度表示表示目标函数的斜率或倾斜度的导数。

为了更好地理解这一点，让我们假设我们在夜间能见度低的山中迷路了。我们只能感觉到脚下地面的坡度。我们的目标是到达山的最低点。要做到这一点，我们必须走几步，向最高坡度的方向移动。我们将迭代地这样做，一步一步地移动，直到我们最终到达山谷。

在数学中，导数是函数在给定点的变化率或斜率。导数的值是特定点的斜率。梯度表示同样的东西，另外它是一个存储偏导数的向量值函数。这意味着梯度是一个向量，它的每个分量都是对特定变量的偏导数。

让我们分析一个函数，f（x，y），即一个双变量函数，x和y。其梯度是一个包含f偏导数的向量：第一个相对于x，第二个相对于y。如果我们计算f的偏导数，我们得到以下结果：

这两个表达式中的第一个称为关于x的偏导数，而第二个偏导数称为关于y。梯度为以下向量：

前面的方程是一个表示二维空间中的点或二维向量的函数。每个分量表示每个函数变量的最陡爬坡方向。因此，梯度指向函数增加最多的方向。

类似地，如果我们有一个五个变量的函数，我们会得到一个五个偏导数的梯度向量。通常，一个带有n变量的函数会产生n维梯度向量，如下所示：

然而，对于梯度下降，我们不希望尽可能快地最大化f。相反，我们想要最小化它，也就是说，找到最小化函数的最小点。

假设我们有一个函数y=f（x）。梯度下降是基于这样的观察，即如果函数 f 在x的邻域中定义且可微，那么如果我们沿着负梯度的方向移动，该函数下降得更快。从x的值开始，我们可以写出以下内容：

在这里，我们有以下几点：

  • 是学习率
  • 是坡度

对于足够小的值，该算法在有限的迭代次数内收敛到函数f的最小值。

重要提示

基本上，如果梯度为负，则该点的目标函数将减小，这意味着参数必须向更大的值移动才能达到最小值。相反，如果梯度为正，则参数会朝较小的值移动，以达到目标函数的较低值。

了解学习率

梯度下降算法通过迭代过程搜索目标函数的最小值。在每一步中，对梯度进行估计，以将下降方向指向使目标函数最小化的方向。在此过程中，学习率参数的选择变得至关重要。此参数决定我们将以多快或多慢的速度移动到目标函数的最佳值：

  • 如果它太小，我们将需要太多的迭代来收敛到最佳值。
  • 如果学习率很高，我们将跳过最优解。

在下图中，您可以看到学习率值施加的两种可能情况：

Figure 7.3 – The scenarios for the learning rate

图 7.3–学习率的场景

因此，必须使用良好的学习率。确定最佳学习率的最佳方法是通过反复试验。

讲解试错法

术语反复试验定义了一种启发式方法，旨在通过尝试并检查是否产生了预期效果来找到问题的解决方案。如果是，则该尝试构成问题的解决方案；否则，我们将继续进行不同的尝试。

让我们分析一下该方法的基本特征：

  • 它面向解决方案：它的目的不是找出一次尝试的原因，而是简单地寻找它。
  • 它是针对所讨论的问题的：它没有权利推广到其他问题。
  • 它不是最优的：它通常局限于寻找一个通常不是最好的解决方案的单一解决方案。
  • 它不需要对它有透彻的了解：它建议找到一个对它知之甚少或一无所知的问题的解决方案。

试错法可以用来找到问题的所有解决方案，或者如果有多个解决方案，则可以找到其中的最佳解决方案。在这种情况下，我们不会在提供期望结果的第一次尝试时停止，而是注意到它并继续尝试，直到找到所有解决方案。最后，根据给定的标准对它们进行比较，从而确定哪一个被认为是最好的。

在 Python 中实现梯度下降

在这一节中，我们将通过完成一个实际示例，将迄今为止所学的应用于梯度下降。我们将定义一个函数，然后使用该方法找到函数的最小点。一如既往，我们将逐行分析代码：

  1. Let's start by importing the necessary libraries:

    import numpy as np
import matplotlib.pyplot as plt

    numpy是一个 Python 库，其中包含许多函数，可帮助我们管理多维矩阵。此外，它还包含大量高级数学函数，我们可以在这些矩阵上使用这些函数。

    matplotlib是用于打印高质量图形的 Python 库。使用 matplotlib，可以通过几个命令生成图形、直方图、条形图、功率谱、误差图、散点图等。它是一组命令行函数，类似于 Matlab 软件提供的那些函数。

  2. Now, let's define the function:

    x = np.linspace(-1,3,100)
y=x**2-2*x+1

    首先，我们为依赖的变量x定义了一个区间。我们只需要将函数可视化并绘制图形。为此，我们使用了linspace()函数。此函数用于创建数字序列。然后，我们传递了三个参数：起点、终点和要生成的点数。接下来，我们定义了一个抛物线函数。

  3. Now, we can draw a graph and display it:

    fig = plt.figure()
axdef = fig.add_subplot(1, 1, 1)
axdef.spines['left'].set_position('center')
axdef.spines['bottom'].set_position('zero')
axdef.spines['right'].set_color('none')
axdef.spines['top'].set_color('none')
axdef.xaxis.set_ticks_position('bottom')
axdef.yaxis.set_ticks_position('left')
plt.plot(x,y, 'r')
plt.show()

    首先，我们定义了一个新图形，然后设置轴，使 x 轴与函数的最小值重合，y 轴与抛物线的中心重合。这将使更容易直观地定位函数的最小点。打印下图：

    Figure 7.4 – The minimum point of the function

    图 7.4–功能的最小点

    如我们所见，y = 0对应的函数的最小值出现在x等于1的值。这将是我们必须通过梯度下降法确定的值。

  4. Now, let's define the gradient function:

    Gradf = lambda x: 2*x-2

    回想一下，函数的梯度是它的导数。在这种情况下，这样做很容易，因为它是一个单变量函数。

  5. Before applying the iterative procedure, it is necessary to initialize a series of variables:

    actual_X  = 3 
learning_rate  = 0.01 
precision_value = 0.000001 
previous_step_size = 1 
max_iteration = 10000 
iteration_counter = 0

    如图所示：

    • actual_X变量将包含自变量x的当前值。首先，我们在x = 3处初始化它，它表示函数在图形中显示范围的最右边的值。
    • learning_rate变量包含学习率。如了解学习率部分所述，设置为0.01。我们可以试着看看如果我们改变这个值会发生什么。
    • precision_value变量将包含定义我们算法精确度的值。作为一个迭代过程，解在每次迭代时都会被细化并趋于收敛。但是，这种收敛可能发生在大量迭代之后，因此为了节省资源，建议在达到足够精度后停止迭代过程。
    • previous_step_size变量将包含该精度的计算，并将初始化为1。
    • max_iteration变量包含我们为算法提供的最大迭代次数。如果该过程不收敛，将使用该值停止该过程。
    • 最后，iteration_counter变量将作为迭代计数器。

  6. Now, we are ready for the iteration procedure:

    while previous_step_size  > precision_value  and iteration_counter  < max_iteration :
    PreviousX = actual_X 
    actual_X  = actual_X  - learning_rate  * Gradf(PreviousX) 
    previous_step_size  = abs(actual_X  - PreviousX) 
    iteration_counter  = iteration_counter +1 
    print('Number of iterations = ',iteration_counter ,'\nActual value of x  is = ',actual_X ) 
    print('X value of f(x) minimum = ', actual_X )

    迭代过程使用一个while循环，该循环将重复自身，直到两个条件都得到验证（TRUE。当两个参数中至少有一个假设为FALSE值时，循环将停止。这两个条件提供了对精度和迭代次数的检查。

    正如定义梯度下降部分中所预期的，该程序要求我们在梯度下降方向上更新x的当前值。我们使用以下等式进行计算：

    在循环的每个步骤中，x的前一个值被存储，以便我们可以计算前一个步骤的精度，作为两个x值之差的绝对值。此外，每一步都会增加步进计数器。在每个步骤结束时，打印迭代次数和x的当前值，如下所示：

    Number of iterations =  520 
Actual value of x  is =  1.0000547758790321
Number of iterations =  521 
Actual value of x  is =  1.0000536803614515
Number of iterations =  522 
Actual value of x  is =  1.0000526067542224
Number of iterations =  523 
Actual value of x  is =  1.000051554619138
Number of iterations =  524 
Actual value of x  is =  1.0000505235267552
Number of iterations =  525 
Actual value of x  is =  1.0000495130562201
Number of iterations =  526 
Actual value of x  is =  1.0000485227950957

    正如我们在每个步骤中所看到的，x的值逐渐接近精确值。这里执行了526迭代。

  7. At the end of the procedure, we can print the result:

    print('X value of f(x) minimum = ', actual_X )

    返回以下结果：

    X value of f(x) minimum =  1.0000485227950957

    我们可以验证，返回值与精确值非常接近，等于1。它与我们作为迭代过程的术语施加的精度值完全不同。

面对牛顿-拉斐逊法

牛顿法是逼近非线性方程根的主要数值方法。该函数在每次迭代时线性近似，以获得更好的零点估计。

使用牛顿-拉斐逊算法寻根

给定一个非线性函数 f 和一个初始近似值x0，牛顿方法通过在x的邻域中为每个k构建函数 f 的线性模型，生成一系列近似值{xk}k>0并用模型本身近似函数。该模型可从泰勒在属于迭代当前点xk 邻域的点x处对函数f的展开开始构建，如下所示：

截断泰勒的一阶发展，我们得到以下线性模型：

前面的方程在xk 的足够小的邻域内仍然有效。

以x0 为初始数据，第一次迭代计算x1 为k=0的前一线性模型的零点，即求解以下标量线性方程：

上一个等式导致下一个迭代的x1，形式如下：

类似地，后续方程迭代x2，其中构造x3，以便我们可以详细说明一般有效性方程，如下所示：

更新方程的形式类似于下降法的通用公式。从几何角度来看，之前的更新方程表示在点（xk，f（xk））处与函数f相切的线。正是由于这个原因，该方法也称为切线方法。

从几何角度来看，我们可以通过以下步骤定义此过程：

  • 函数的切线在起点x0 处绘制。
  • 这条线的截距用 x 轴标识。此点表示新值x1。
  • 重复此过程直到收敛。

下图显示了此过程：

Figure 7.5 – The procedure of finding a tangent

图 7.5–寻找切线的程序

如果每k一个f'（xk*）=0*，则该算法定义良好。关于计算成本，可以注意到，在每次迭代中，需要对函数f及其在点xk 之前的导数进行求值。

逼近牛顿-拉斐逊进行数值优化

牛顿-拉斐逊法也用于解决数值优化问题。在这种情况下，该方法采用牛顿法求函数零点的形式，但适用于函数f的导数。这是因为确定函数f的最小点等同于确定一阶导数*f'*的根。

在这种情况下，更新公式采用以下形式：

在上一个等式中，我们有以下等式：

  • 是函数f的一阶导数

  • is the second derivative of the function f

    重要提示

    牛顿-拉斐逊法由于其速度快，通常优于下降梯度法。然而，它需要了解一阶导数和二阶导数的解析表达式，并且不分青红皂白地收敛到极小值和极大值。

有一些变体使该方法全局收敛，并通过避免使用直接方法确定研究方向来降低计算成本。

应用牛顿-拉斐逊技术

在本节中，我们将通过完成一个实践练习来应用到目前为止所学的牛顿-拉斐逊法。我们将定义一个函数，然后使用该方法找到函数的最小点。一如既往，我们将逐行分析代码：

  1. 让我们从导入必要的库开始：

    import numpy as np
import matplotlib.pyplot as plt

  2. Now, let's define the function:

    x = np.linspace(0,3,100)
y=x**3 -2*x**2 -x + 2

    首先，我们为因变量 x 定义了一个区间。为了绘制图形，我们只需要将函数可视化。为此，我们使用了linspace()函数。此函数用于创建数字序列。然后，我们传递了三个参数：起点、终点和要生成的点数。接下来，我们定义了一个三次函数。

  3. Now, we can draw a graph and display it:

    fig = plt.figure()
axdef = fig.add_subplot(1, 1, 1)
axdef.spines['left'].set_position('center')
axdef.spines['bottom'].set_position('zero')
axdef.spines['right'].set_color('none')
axdef.spines['top'].set_color('none')
axdef.xaxis.set_ticks_position('bottom')
axdef.yaxis.set_ticks_position('left')
plt.plot(x,y, 'r')
plt.show()

    首先，我们定义了一个新图形，然后设置轴，使 x 轴与函数的最小值重合，y 轴与抛物线的中心重合。这将更容易直观地定位函数的最小点。打印以下图表：

    Figure 7.6 – The minimum point of the function

    图 7.6–功能的最小点

    在这里，我们可以看到函数的最小值出现在x值大致等于1.5的情况下。这将是我们必须通过梯度下降法确定的值。但是为了获得精确的值，以便我们可以将其与稍后得到的值进行比较，我们需要提取该值：

    print('Value of x at the minimum of the function', x[np.argmin(y)])

    为了确定这个值，我们使用了 NumPy 的argmin()函数。此函数返回向量中最小元素的位置索引。返回以下结果：

    Value of x at the minimum of the function 1.5454545454545454

  4. 现在，让我们定义一阶和二阶导数函数：

    FirstDerivative = lambda x: 3*x**2-4*x -1 
SecondDerivative = lambda x: 6*x-4

  5. Now, we will initialize some parameters:

    actual_X  = 3 
precision_value  = 0.000001 
previous_step_size  = 1 
max_iteration  = 10000 
iteration_counter  = 0

    这些参数具有以下含义：

    • actual_X变量将包含自变量x的当前值。首先，我们在x = 3处初始化它，它表示函数在图形中显示范围的最右边的值。
    • precision_value变量将包含定义我们算法精确度的值。作为一个迭代过程，解在每次迭代时都会被细化并趋于收敛。但是，这种收敛可能发生在大量迭代之后，因此为了节省资源，建议在达到足够精度后停止迭代过程。
    • previous_step_size变量将包含该精度的计算，并将初始化为1。
    • max_iteration变量包含我们为算法提供的最大迭代次数。如果该过程不收敛，将使用该值停止该过程。
    • 最后，iteration_counter变量将作为迭代计数器。

  6. Now, we can apply the Newton-Raphson method, as follows:

    while previous_step_size  > precision_value  and iteration_counter  < max_iteration :
    PreviousX = actual_X 
    actual_X  = actual_X  - FirstDerivative(PreviousX)/ SecondDerivative(PreviousX)
    previous_step_size  = abs(actual_X  - PreviousX) 
    iteration_counter  = iteration_counter +1 
    print('Number of iterations = ',iteration_counter ,'\nActual value of x  is = ',actual_X )

    这个过程类似于我们在 Python 部分中实现梯度下降时所采用的解决问题的方法。这里使用了一个while循环，它将重复自身，直到两个条件都得到验证（TRUE。当两个参数中至少有一个假设为FALSE值时，循环将停止。这两个条件提供了对精度和迭代次数的检查。

    Newton-Raphson 方法更新 x 的当前值，如下所示：

    在循环的每一步中，存储前一步的x值，以计算前一步的精度，作为两个x值之差的绝对值。此外，每一步都会增加步进计数器。在每一步结束时，打印迭代次数和x的当前值，如下所示：

    Number of iterations =  1 
Actual value of x  is =  2.0
Number of iterations =  2 
Actual value of x  is =  1.625
Number of iterations =  3 
Actual value of x  is =  1.5516304347826086
Number of iterations =  4 
Actual value of x  is =  1.5485890147300967
Number of iterations =  5 
Actual value of x  is =  1.5485837703704566
Number of iterations =  6 
Actual value of x  is =  1.5485837703548635

    正如我们在数值优化的接近牛顿-拉斐逊一节中提到的，达到该解所需的迭代次数严重扭曲。事实上，我们从基于梯度下降的方法收敛所需的 526 次迭代，到牛顿-拉斐逊方法的 6 次迭代。

  7. Finally, we will print the result:

    print('X value of f(x) minimum = ', actual_X )

    返回以下结果：

    X value of f(x) minimum =  1.5485837703548635

    我们可以验证，返回值与精确值非常接近，等于1.5454545454545454。它与我们作为迭代过程的术语施加的精度值完全不同。

加深对随机梯度下降的认识

正如我们在探索梯度下降技术一节中提到的，梯度下降方法的实施包括从维度空间中随机选择的配置开始，对函数及其梯度进行初始评估。

从这里开始，我们尝试沿着梯度指示的方向移动。这确定了下降方向，其中函数趋向于最小值，并检查函数是否实际采用低于先前配置中计算的值。如果是这样，该过程将继续迭代，重新计算新的渐变。这可能与前一个完全不同。在此之后，它再次开始搜索新的最小值。

此迭代过程要求在每个步骤中更新整个系统状态。这意味着必须重新计算系统的所有参数。从计算角度来看，这相当于极其昂贵的运营成本，并大大减慢了估算过程。对于标准梯度下降法，在计算整个数据集的梯度后更新权重，在随机方法中，在一定数量的示例后更新系统参数。这些是随机选择的，以加快过程，并尝试避免任何局部最小情况。

考虑包含 N 个观测现象的数据集。这里，让f作为一个目标函数，我们想要最小化一系列参数x。在这里，我们可以写出以下等式：

通过对上一个等式的分析，我们可以推断目标函数f的评估需要对函数f进行 n 次评估，数据集中包含的每个值一次评估。

在经典梯度下降法中，在每一步中，通过以下等式计算与数据集所有值对应的函数梯度：

在某些情况下，对上一个等式中的和的评估可能特别昂贵，例如当数据集特别大且目标函数没有基本表达式时。梯度的随机下降通过引入梯度函数的近似来解决这个问题。在每一步中，梯度的计算仅在数据集的随机子集中使用，而不是与数据集中包含的数据相对应地计算梯度的总和。

因此，前面的等式替换了以下等式：

在前面的等式中，是数据集中随机选择的一个观测值的梯度。

这种技术的优点如下：

  • 该算法仅基于一部分观测值，允许对参数空间进行更广泛的探索，有更大的可能性找到新的和可能更好的最小点。
  • 该算法的一个步骤在计算上要快得多，从而确保更快地收敛到最小点。
  • 参数估计值也可以通过一次只将数据集的一部分加载到内存中来计算，从而允许将此方法应用于大型数据集。

在 Python 中发现多元优化方法

在本节中，我们将分析 Python SciPy 库中包含的一些数值优化方法。SciPy 是基于 NumPy 的数学算法和函数的集合。它包含一系列可用于操作和显示数据的命令和高级类。有了 SciPy，Python 中增加了功能，使其成为一个数据处理和系统原型环境，类似于商业系统，如 MATLAB。

使用 SciPy 的科学应用程序受益于世界各地开发人员在数值计算的众多领域开发的附加模块。数值优化问题也包含在可用模块中。

SciPyoptimize模块包含许多用于最小化/最大化目标函数的函数，包括约束函数和非约束函数。它处理非线性问题，同时支持局部和全局优化算法。此外，还讨论了线性规划、约束和非线性最小二乘、寻根和曲线自适应等问题。在下面的部分中，我们将分析其中的一些。

内尔德-米德法

大多数著名的优化算法都是基于导数的概念和可以从梯度中推断出的信息。然而，从实际应用中衍生出的许多优化问题的特点是，目标函数的解析表达式未知，这使得无法计算其导数，或者由于其特别复杂，因此对导数进行编码可能需要太长时间。为了解决这类问题，已经开发了几种算法，这些算法不尝试近似梯度，而是使用一组采样点中的函数值通过其他方式确定新的迭代。

Nelder-Mead方法通过评估构成称为单纯形的几何形状的测试点，尝试最小化非线性函数。

重要提示

单纯形被定义为欧几里德空间中的一组闭凸点，使我们能够找到线性规划的典型优化问题的解。

选择单纯形的几何图形主要是由于两个原因：单纯形能够使其形状适应目标函数自身变形的空间趋势，以及它只需要记忆n+1点。基于单纯形的直接搜索方法的每次迭代都以单纯形开始，由其n+1顶点和相关函数的值指定。计算一个或多个测试点和函数的相应值，迭代以一个新的单纯形结束，以便函数在其顶点处的值满足与先前单纯形相关的某种形式的下降条件。

Nelder-Mead 算法在每次迭代时对函数的求值方面特别节省，因为在实践中，它通常只需要对函数进行一次或两次求值即可构建新的单纯形。但是，由于它不使用任何梯度评估，因此可能需要更长的时间才能找到最小值。

使用 SciPyoptimize模块的minimize例程，该方法在 Python 中很容易实现。让我们看一个使用此方法的简单示例：

  1. Let's start by loading the necessary libraries:

    import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
from mpl_toolkits.mplot3d import Axes3D

    导入生成 3D 图形所需的库（Axes3D。

  2. Now, let's define the function:

    def matyas(x):
   return 0.26*(x[0]**2+x[1]**2)-0.48*x[0]*x[1]

    Matyas 函数是连续的、凸的、单峰的、可微的和不可分离的，定义在二维空间上。matyas功能定义如下：

    这个函数是在 x，ye[-10，10]上定义的。此函数在 f（0，0）=0 中有一个全局最小值。

  3. Let's visualize the matyas function:

    x = np.linspace(-10,10,100)
y = np.linspace(-10,10,100)
x, y = np.meshgrid(x, y)
z = matyas([x,y])
fig = plt.figure()
ax = fig.gca(projection='3d')
surf = ax.plot_surface(x, y, z, rstride=1, cstride=1, 
                      cmap=cm.RdBu,linewidth=0, antialiased=False)
ax.zaxis.set_major_locator(LinearLocator(10))
ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))
fig.colorbar(surf, shrink=0.5, aspect=10)
plt.show()

    首先，我们定义了自变量 x 和 y，在我们已经指定的[-10.10]范围内。因此，我们使用numpy meshgrid()函数创建了一个网格。此函数创建一个数组，其中行和列对应于x和y的值。我们将使用该矩阵绘制z变量的对应点，该变量对应于 Matyas 函数。在定义了x、y和z变量之后，我们绘制了一个三维图形来表示函数。绘制了以下图表：

    Figure 7.7 – Meshgrid plot to represent the function

    图 7.7–表示功能的网格图

  4. As we already mentioned, the Nelder-Mead method does not require us to calculate a derivative as it is limited to evaluating the function. This means that we can directly apply the method:

    x0 = np.array([-10, 10])
NelderMeadOptimizeResults = minimize(matyas, x0,
            method='nelder-mead',
            options={'xatol': 1e-8, 'disp': True})
print(NelderMeadOptimizeResults.x)

    为了做到这一点，我们首先定义了一个初始点，以开始搜索函数的最小值。因此，我们使用了 SciPy 优化模块的minimize()功能。此函数查找一个或多个变量的标量函数的最小值。已传递以下参数：

    • matyas：您希望最小化的功能

    • x0：初始向量

    • method = 'nelder-mead': The method used for the minimization procedure

      此外，还添加了以下两个选项：

    • 'xatol': 1e-8：定义收敛可接受的绝对误差

    • 'disp': True：设置为 True 可打印聚合消息

  5. Finally, we printed the results of the optimization method, as follows:

    Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 77
         Function evaluations: 147
[3.17941614e-09 3.64600127e-09]

    正如已经预期的那样，最小值是在值0中确定的。此外，该值与以下值一致：

    X = 3.17941614e-09 
Y = 3.64600127e-09

    正如我们所预期的，这些值非常接近于零。与该值的偏差与我们为该方法设置的误差一致。

鲍威尔共轭方向算法

共轭方向方法最初是作为迭代方法引入的，用于求解具有对称正定系数矩阵的线性系统，以及最小化严格凸二次函数。

用于最小化二次函数的共轭方向方法的主要特征是，以简单的方式生成一组方向，这些方向除了线性独立外，还具有相互共轭的更重要的特性。

鲍威尔方法的思想是，如果沿着 p（p< n) directions in a stage of the research, then when taking a step along each direction, the final displacement from the beginning up to the p的每一步找到二次函数的最小值，则第四步相对于研究的所有p子方向共轭。

例如，如果点 1 和 2 是从同一方向但不同起点的一维搜索中获得的，则由 1 和 2 形成的线将指向最大值。这些线表示的方向称为共轭方向。

让我们分析一个应用鲍威尔方法的实际案例。我们将使用我们在内尔德-米德方法一节中定义的matyas函数：

  1. 让我们从加载必要的库开始：

    import numpy as np
from scipy.optimize import minimize

  2. 现在，让我们定义函数：

    def matyas (x):
    return 0.26 * (x [0] ** 2 + x [1] ** 2) -0.48 * x [0] * x [1]

  3. Now, let's apply the method:

    x0 = np.array([-10, 10])
PowellOptimizeResults = minimize(matyas, x0, 
           method='Powell',
           options={'xtol': 1e-8, 'disp': True})
print(PowellOptimizeResults.x)

    此处使用了 SciPy 优化模块的minimize()功能。此函数查找一个或多个变量的标量函数的最小值。已传递以下参数：

    • matyas：我们想要最小化的功能

    • x0：初始向量

    • method = 'Powell': The method used for the minimization procedure

      此外，还添加了以下两个选项：

    • 'xtol': 1e-8：定义收敛可接受的绝对误差

    • 'disp': True:设置为真打印汇聚消息

  4. Finally, we printed the results of the optimization method. The following results are returned:

    Optimization terminated successfully.
        Current function value: 0.000000
        Iterations: 3
        Function evaluations: 66
[-6.66133815e-14 -1.32338585e-13]

    根据尼尔-米德方法一节的规定，最小值为 0。此外，该值与以下值一致：

    X = -6.66133815e-14
Y = -1.32338585e-13

    正如我们所预期的，这些值非常接近于零。现在，我们可以对应用于同一函数的两种方法进行比较。我们可以注意到，对于 Powell 方法，达到收敛所需的迭代次数等于3，而对于 Nelder-Mead 方法，它等于77。还注意到对该功能的评估数量急剧减少；66反对147。最后，通过 Powell 方法减小计算值与期望值之间的差异。

总结其他优化方法

SciPy optimize 软件包中的minimize()例程包含许多用于无约束和约束最小化的方法。我们在前面的章节中详细分析了其中的一些。在以下列表中，我们总结了该软件包提供的最常用方法：

  • Newton Broyden-Fletcher Goldfarb-Shanno（BFGS）：这是一种迭代无约束优化方法，用于解决非线性问题。此方法查找一阶导数为零的点。
  • 共轭梯度****CG：该方法属于共轭梯度算法家族，执行一个或多个变量的标量函数最小化。这种方法要求系统矩阵是对称的、正定的。
  • 狗腿信赖域（狗腿）：方法首先围绕当前最佳解定义一个区域，在该区域可以近似原目标函数。因此，该算法在该区域内向前迈出了一步。
  • 牛顿 CG：这种方法也称为截断牛顿法。这是一种方法，通过采用基于共轭梯度的程序确定研究方向，以粗略最小化二次函数。
  • 有限记忆 BFGS****L-BFGS：这是拟牛顿方法家族的一部分。它使用 BFGS 方法系统地保存计算机内存。
  • 线性近似约束优化（COBYLA）：运行机制为迭代，利用线性规划的原理对上一步发现的解进行细化。通过逐步降低速度来实现收敛。

总结

在本章中，我们学习了如何使用不同的数值优化技术来改进仿真模型提供的解决方案。我们首先介绍了数值优化的基本概念，定义了一个极小化问题，并学习了如何区分局部极小值和全局极小值。然后，我们继续研究基于梯度下降的优化技术。我们定义了该技术的数学公式，并给出了几何表示。此外，我们还加深了对学习率和试错相关概念的了解。通过这样做，我们解决了一个实际案例，以加强我们通过解决搜索二次函数最小值的问题所学到的概念。

随后，我们学习了如何使用 Newton-Raphson 方法搜索函数的根，以及如何利用相同的方法进行数值优化。我们还分析了这项技术的一个实际案例，以立即将我们学到的概念付诸实践。我们通过寻找凸函数的局部极小值来实现这一点。

然后，我们继续研究随机梯度下降算法，这使我们能够大大降低数值优化问题的计算成本。这一结果是通过在每一步使用单一的梯度估计来获得的，该梯度是以随机方式从可用梯度中选择的。

最后，我们探讨了 pythonscipy 包中包含的多元数值优化算法。对于其中一些，我们定义了数学公式，并给出了使用该方法的一个实例。对于其他人，起草了一份摘要，列出他们的特点。

在下一章中，我们将学习如何使用仿真模型来处理财务问题。我们将探索几何布朗运动模型是如何工作的，我们将发现如何使用蒙特卡罗方法进行股票价格预测。最后，我们将学习如何使用马尔可夫链建模信用风险。

八、将仿真模型用于金融工程

基于人工智能和机器学习的系统的爆炸性进入为金融部门开辟了新的前景。这些方法可以带来诸如用户权利保护等好处，以及宏观经济效益。

蒙特卡罗方法在金融领域找到了一个自然的应用，可以用数值方法解决有保障看涨期权的定价和问题。从本质上说，这些方法包括使用给定的数学定律和足够大的数据集仿真给定的过程或现象，这些数据是从充分代表真实变量的分布中随机创建的。其想法是，如果不可能进行分析研究，或者不可能或不方便进行充分的实验取样，则使用现象的数值仿真。在本章中，我们将介绍在金融环境中使用仿真方法的实际案例。您将学习如何使用蒙特卡罗方法预测股票价格，以及如何评估与股票组合相关的风险。

在本章中，我们将介绍以下主题：

  • 理解几何布朗运动模型
  • 用蒙特卡罗方法预测股票价格
  • 投资组合管理的风险模型研究

技术要求

在本章中，我们将学习如何在金融工程中使用仿真模型。为了理解这些主题，需要代数和数学建模的基础知识。

要使用本章中的 Python 代码，您需要以下文件（可在 GitHub 的上获得）https://github.com/PacktPublishing/Hands-On-Simulation-Modeling-with-Python ：

  • StandardBrownianMotion.py
  • AmazonStockMontecarloSimulation.py
  • ValueAtRisk.py

理解几何布朗运动模型

苏格兰植物学家罗伯特·布朗于 1827 年在显微镜下观察了悬浮在水中的花粉颗粒如何以随机和不可预测的方式连续运动，并由此命名为布朗是。1905 年，爱因斯坦对布朗观察到的运动现象给出了分子解释。他认为粒子的运动是可以用数学描述的，假设各种跳跃是由于花粉粒子与水分子的随机碰撞造成的。

今天，布朗运动首先是概率论背景下的一种数学工具。这一数学理论被用来描述一系列不断扩大的现象，这些现象是由与物理学截然不同的学科研究的。例如，金融证券的价格、热的传播、动物种群、细菌、疾病、声音和光线都是使用同一工具建模的。

重要提示

布朗运动是一种由小颗粒或胶粒大小的颗粒所产生的不间断和不规则运动组成的现象，也就是说，这些颗粒太小，肉眼无法观察到，但在浸入流体中时比原子大得多。

定义标准布朗运动

构造布朗运动模型有多种方法和布朗运动的各种等价定义。让我们从标准布朗运动（维纳过程）的定义开始。标准布朗运动的基本性质包括：

  • 标准布朗运动从零开始。
  • 标准布朗运动采用连续路径。
  • 布朗过程所遭受的增加是独立的。
  • 布朗过程在时间间隔dt中所受的增加表示高斯分布，平均值等于零，方差等于时间间隔dt。

基于这些性质，我们可以把这个过程看作是大量极小增量的总和。在选择两个瞬间t和s后，随机变量Y（s）-Y（t）遵循正态分布，mean 为（s-t），方差为（s-t），我们可以用以下等式表示：

正态性假设在线性变换中非常重要。事实上，标准布朗运动的名称来源于标准正态分布的分布类型，参数为=0 和=1。

因此，可以说，具有单位均值和方差的布朗运动*Y（t）*可以表示为标准布朗运动的线性变换，根据以下等式：

在前面的等式中，我们可以观察到以下情况：

  • 是标准的布朗运动。

该方程的弱点在于*Y（t）假设负值的概率为正值；事实上，由于Z（t）具有独立增量的特征，可以采用负号，因此Y（t）*的负性风险不是零。

现在，考虑足够小的时间间隔的 Brownian 运动（Wiener 过程）。该过程的无穷小增量如下表所示：

前面的公式可以改写如下：

这一过程不受变化的限制，因此在经典分析的背景下无法区分。事实上，对于区间dt，前一个趋于无穷大。

将维纳过程称为随机游动

维纳过程可被视为随机游动的边界情形。我们在第 5 章**基于仿真的马尔可夫决策过程中处理了一个随机游走。我们已经看到，一个粒子在瞬间n的位置将由以下等式表示：

*

在前面的公式中，我们可以观察到以下情况：

  • Yn 是行走中的下一个值。
  • Yn-1 为前一时间段的观测值。
  • Zn 是该步骤中的随机波动。

如果n随机数Zn 的均值等于零，方差等于 1，那么对于n的每个值，我们可以使用以下等式定义一个随机过程：

上述公式可用于迭代过程。对于非常大的n值，我们可以编写以下内容：

前面的公式是由于我们在第 4 章、蒙特卡罗仿真中介绍的中心极限定理。

执行标准布朗运动

那么，让我们演示如何在 Python 环境中生成简单的布朗运动。让我们从最简单的情况开始，在这种情况下，我们定义时间间隔、要执行的步骤数和标准偏差：

  1. We start by importing the following libraries:

    import numpy as np
import matplotlib.pyplot as plt

    numpy库是一个 Python 库，包含许多函数，可以帮助我们管理多维矩阵。此外，它还包含大量高级数学函数，我们可以使用这些函数对这些矩阵进行运算。

    matplotlib库是用于打印高质量图形的 Python 库。有了matplotlib，只需几个命令就可以生成图形、直方图、条形图、功率谱、误差图、散点图等。这包括一组类似于 MATLAB 软件提供的命令行函数。

  2. Now, let's proceed with some initial settings:

    np.random.seed(4)
n = 1000
SQN = 1/np.math.sqrt(n)
ZValues = np.random.randn(n)
Yk = 0
SBMotion=list()

    在第一行代码中，我们使用random.seed()函数初始化随机数生成器的种子。这样，使用随机数的仿真将是可复制的。实验的再现性将通过以下事实得到保证：生成的随机数始终相同。我们设置迭代次数（n，并计算以下等式的第一项：

    然后，我们使用random.randn()函数生成n随机数。此函数返回平均值为 0、方差为 1 的n样本的标准正态分布。最后，我们根据属性*（Y*（0）*=0）*的要求设置布朗运动的第一个值，并初始化包含布朗运动位置坐标的列表。

  3. At this point, we will use a for loop to calculate all of the n positions:

    for k in range(n):
    Yk = Yk + SQN*ZValues[k]
    SBMotion.append(Yk)

    我们只需将当前随机数乘以 SQN，添加到包含累积和的变量中。然后将当前值附加到 SBMotion 列表中。

  4. Finally, we draw a graph of the Brownian motion created:

    plt.plot(SBMotion)
plt.show()

    打印以下图表：

Figure 8.1: Brownian motion graph

图 8.1–布朗运动图

因此，我们创建了我们的第一个布朗运动仿真。它的使用特别适合于金融仿真。在下一节中，我们将演示如何做到这一点。

利用蒙特卡罗方法进行股价预测

正如我们在第 4 章、**蒙特卡罗仿真中所探讨的，蒙特卡罗方法使用事件在特定条件下可能发生的不同概率来仿真受检过程的不同演变。这些仿真探索了现象的整个参数空间，并返回了一个具有代表性的样本。对于获得的每一个样本，都要对感兴趣的数量进行测量，以评估其性能。正确的仿真意味着过程结果的平均值收敛到预期值。

探索亚马逊股价走势

股票市场提供了一个快速赚大钱的机会，也就是说，至少在一个没有经验的用户眼中是这样。股票市场上的交易所可以引起价格的大幅波动，吸引世界各地投机者的注意。为了从股票市场的投资中获得收入，有必要从多年对这一现象的深入研究中获得扎实的知识。在这种情况下，拥有一种预测股票市场证券的工具的可能性代表了所有人都感到的需要。

让我们演示如何开发世界上最著名公司之一的股票仿真模型。亚马逊由杰夫·贝佐斯（Jeff Bezos）在 20 世纪 90 年代创立，是世界上最早通过互联网销售产品的公司之一。自 1997 年以来，亚马逊股票在证券交易所上市，代号为 AMZN。AMZN 股票的历史价值可以从过去 10 年来一直在处理股票市场的各种互联网网站上获得。我们将参考 AMZN 股票在2010-04-08至2020-04-07纳斯达克 GS 股票报价中的表现。为了从2020-04-07获取数据，我们需要选择雅虎网站上的2020-04-08作为结束日期。

可从雅虎财经网站下载.csv格式的数据 https://finance.yahoo.com/quote/AMZN/history/ 。

在以下屏幕截图中，您可以看到 AMZN 股票的 Yahoo Finance 部分，其中突出显示的按钮用于下载数据：

Figure 8.2: Amazon data on Yahoo Finance

图 8.2–雅虎财经的亚马逊数据

下载的AMZN.csv文件包含很多功能，但我们只使用其中两个，如下所示：

  • 日期：报价日期
  • 收盘：收盘价

我们将逐行分析代码，以充分了解整个过程，这将引导我们仿真亚马逊股价表现的一系列预测场景：

  1. As always, we start by importing the libraries:

    import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import norm
from pandas.plotting import register_matplotlib_converters
register_matplotlib_converters()

    已导入以下库：

    pandas库是一个开放源码的 BSD 许可库，它包含用于为 Python 编程语言处理高性能数值的数据结构和操作。

    SciPy 是基于 NumPy 的数学算法和函数的集合。它有一系列命令和高级类来操作和显示数据。有了 SciPy，Python 就增加了功能，使其成为一个类似于商业系统（如 MATLAB）的数据处理和系统原型环境。

  2. Now, let's import the data contained in the AMZN.csv file:

    AmznData = pd.read_csv('AMZN.csv',header=0, 
           usecols = ['Date',Close'],parse_dates=True,
           index_col='Date')

    我们使用了pandas库的read_csv模块，该模块将数据加载到名为DataFrame的pandas对象中。传递以下参数：

    'AMZN.csv'：文件名。

    header=0：包含列名和数据开头的行号。默认情况下，如果传递了非标题行（header=0），则从文件的第一行推断列名。

    usecols=['Date',Close']：此参数通过指定列名提取数据集的子集。

    parse_dates=True：布尔值；如果True，请尝试解析索引。

    index_col='Date'：这允许我们指定将用作数据帧索引的列的名称。

  3. Now we will explore the imported dataset to extract preliminary information. To do this, we will use the info() function, as follows:

    print(AmznData.info())

    打印以下信息：

    <class 'pandas.core.frame.DataFrame'>
DatetimeIndex: 2518 entries, 2010-04-08 to 2020-04-07
Data columns (total 1 columns):
Close    2518 non-null float64
dtypes: float64(1)
memory usage: 39.3 KB
None

    在这里，返回了很多有用的信息：对象类、存在的记录数（2518）、索引的起始值和结束值（2010-04-08到2020-04-07）、列数及其包含的数据类型，以及其他信息。

    我们还可以打印数据集的前五行，如下所示：

    print(AmznData.head())

    打印以下数据：

                     Close
Date                  
2010-04-08  140.960007
2010-04-09  140.059998
2010-04-12  141.199997
2010-04-13  140.160004
2010-04-14  144.279999

    如果我们想要打印不同数量的记录，那么通过指示要打印的行数就足够了。同样，我们可以打印数据集的最后 10 条记录：

    print(AmznData.tail())

    打印以下记录：

                      Close
Date                   
2020-04-01  1907.699951
2020-04-02  1918.829956
2020-04-03  1906.589966
2020-04-06  1997.589966
2020-04-07  2011.599976

    通过对头部和尾部的初步快速比较，我们可以验证过去 10 年亚马逊股票的价值从 140 美元左右上升到 2011 美元左右。这对亚马逊股东来说是一笔极好的交易。

    使用describe()功能，我们将使用基本统计数据提取数据预览：

    print(AmznData.describe())

    返回以下结果：

                 Close
count  2518.000000
mean    723.943074
std     607.588565
min     108.610001
25%     244.189995
50%     398.995011
75%    1006.467514
max    2170.219971

    我们可以确认在过去 10 年中价值的显著增长，但我们也可以看到，鉴于标准差的值非常高，股票如何经历了显著的波动。这告诉我们，那些忠于股票并长期持有股票的股东从股票上涨中获益最大。

  4. After analyzing the preliminary data statistics, we can take a look at the performance of the Amazon shares in the last 10 years by drawing a simple graph:

    plt.figure(figsize=(10,5))
plt.plot(AmznData)
plt.show()

    使用了以下matplotlib功能：

    figure()：此函数创建一个新的图形，该图形暂时为空。我们使用figsize参数设置框架的大小，该参数以英寸为单位设置宽度和高度。

    plot()：此函数用于绘制AmznData数据集。

    show()：此功能在 IPython 和 PyLab 模式下运行时，显示所有图形并返回 IPython 提示符。

    打印以下图：

Figure 8.3: Amazon share graph

图 8.3–亚马逊份额图

亚马逊股票在过去 10 年中的显著增长是显而易见的。此外，应该注意的是，自 2015 年以来，增长幅度最大，但我们将尝试从数据中提取更多信息。

将股价走势作为时间序列处理

亚马逊股票价格随时间的趋势，如上图所示，配置为有序数据序列。这类数据可以方便地作为时间序列处理。让我们考虑一个简单的定义：时间序列包含一个变量的实验观测的时间序列。此变量可以与不同来源的数据相关。通常，它涉及失业率、利差、股市指数和股价趋势等金融数据。

将该问题作为时间序列处理将使我们能够从数据中提取有用的信息，以便开发用于管理未来情景的预测模型。比较不同年份同一时期的股票价格趋势，或者更简单地说，比较连续时期之间的股票价格趋势，可能很有用。

假设Y1、Yt、Yn 是时间序列的元素。让我们首先比较两个不同时间的数据，分别用t和t+1表示。因此，它是两个连续的时期；我们感兴趣的是评估被观察现象所经历的变化，其可通过以下比率定义：

这个百分比称为百分比变化。可定义为时间t+1的Y相对于上一次t的百分比变化率。此描述符返回有关数据在一段时间内如何发生更改的信息。百分比变化允许您监控股票价格和市场指数，而不仅仅是比较不同国家的货币：

  1. To evaluate this useful descriptor, we will use the pct_change() function contained in the pandas library:

    AmznDataPctChange = AmznData.pct_change()

    此函数返回当前元素和上一个元素之间的百分比变化。默认情况下，该函数计算前一行的百分比变化。

    时间序列百分比变化的概念与股票价格回报的概念相联系。基于回报的方法允许数据的标准化，这是在评估以不同度量为特征的变量之间的关系时具有根本重要性的操作。

    我们将在对数尺度上处理收益，因为这种选择将给我们带来几个优势：正态分布结果；返回的值（返回的对数）非常接近初始值（返回），至少对于非常小的值；随着时间的推移，会产生累加的结果。

  2. To pass the return on a logarithmic scale, we will use the log() function of the numpy library, as follows:

    AmznLogReturns = np.log(1 + AmznDataPctChange) 
print(AmznLogReturns.tail(10))

    打印以下结果：

                   Close
Date                
2020-03-25 -0.028366
2020-03-26  0.036267
2020-03-27 -0.028734
2020-03-30  0.033051
2020-03-31 -0.007272
2020-04-01 -0.021787
2020-04-02  0.005817
2020-04-03 -0.006399
2020-04-06  0.046625
2020-04-07  0.006989

  3. 为了更好地理解收益是如何随时间分布的，让我们画一个图表：

    plt.figure(figsize=(10,5))
plt.plot(AmznLogReturns)
plt.show()

打印以下图表：

Figure 8.4: Logarithmic values of the returns

图 8.4–收益的对数值

前面的图表显示，对数收益率在整个周期内呈正态分布，平均值稳定。

介绍 Black-Scholes 模型

布莱克-斯科尔斯（BS）模型无疑代表了量化金融史上最重要、最具革命性的工作。在传统金融文献中，假设几乎所有金融资产价格（股票、货币和利率）都是由布朗漂移运动驱动的。

该模型假设资产的预期收益等于无风险利率r。这种方法能够仿真资产对数规模的回报。假设我们在瞬间观察到一个资产：t（0），t（1）…，t（n）。我们注意到，使用s（i）=s（ti），资产在*t（i）*的价值。基于这些假设，我们可以使用以下等式计算回报：

然后，我们将对数标度上的收益率进行转换，如下所示：

通过将 BS 方法应用于布朗几何运动，股票价格将满足以下随机微分方程：

在前面的方程中，dB（t）是标准布朗运动，μ和σ是实常数。前面的方程在假设s（i）-s（i-1）很小的情况下有效，当股票价格发生轻微变化时，就会出现这种情况。这是因为如果z很小，ln*（1+z）大致等于z*。上一个方程的解析解如下：

通过在对数尺度上传递前面的方程，我们得到以下方程：

在前面的等式中，我们可以观察到以下情况：

  • α是漂移。
  • B（t）是标准布朗运动。
  • 为标准偏差。

我们引入了漂移的概念，它代表了股票市场中长期资产的趋势。为了理解漂移，我们将使用河流水流的概念。如果我们把液体颜色倒入河流中，它会沿着河流水流的方向扩散。类似地，漂移表示股票跟随长期资产趋势的趋势。

应用蒙特卡罗仿真

使用上一节讨论的 BS 模型，我们可以评估从前一天开始的资产日价格乘以基于系数r的指数贡献。这个系数是一个周期性的回报率。它转化为以下等式：

前面的式中的第二项，呃，叫做日收益率，根据 BS 模型，由下式给出：

无法预测资产的回报率。表示它的唯一方法是把它看作一个随机数。因此，为了预测资产的价格趋势，我们可以使用基于随机运动的模型，如 BS 方程所表示的模型。

BS 模型假设股票价格的变化取决于一段时间内的预期回报。日收益率有两项：固定漂移率和随机变量。这两个术语规定了波动的确定性和不确定性。

为了计算漂移，我们将使用预期收益率，即最可能发生的收益率，使用对数收益和方差的历史平均值，如下所示：

根据前面的等式，资产的日变化率是收益的平均值，小于一段时间内方差的一半。我们继续我们的工作，计算在处理股票价格趋势作为时间序列部分中计算的亚马逊证券回归的漂移：

  1. To evaluate the drift, we need the mean and variance of the returns. Since we also calculate the standard deviation, we will need the calculation of the daily return:

    MeanLogReturns = np.array(AmznLogReturns.mean())
VarLogReturns = np.array(AmznLogReturns.var()) 
StdevLogReturns = np.array(AmznLogReturns.std())

    使用了三个numpy功能：

    mean()：计算沿指定轴的算术平均值，并返回数组元素的平均值。

    var()：计算沿指定轴的方差。它返回数组元素的方差，该方差是分布范围的度量。

    std()：计算沿指定轴的标准偏差。

    现在我们可以计算漂移如下：

    Drift = MeanLogReturns - (0.5 * VarLogReturns)
print("Drift = ",Drift)

    返回以下结果：

    Drift =  [0.00086132]

    这是布朗运动的固定部分。漂移返回预期值的年化变化，并补偿与直接布朗运动相比结果的不对称性。

  2. 为了计算布朗运动的第二个分量，我们将使用随机变量。这对应于平均值和事件之间的距离，表示为标准偏差的数量。在执行此操作之前，我们需要设置间隔和迭代次数。间隔的数量将等于观测的数量（2518），而迭代的数量（代表我们打算开发的仿真模型的数量）为 20:

    NumIntervals = 2518
Iterations = 20

  3. Before generating random values, it is recommended that you set the seed to make the experiment reproducible:

    np.random.seed(7)

    现在，我们可以生成随机分布：

    SBMotion = norm.ppf(np.random.rand(NumIntervals, Iterations))

    返回一个 2518×20 矩阵，包含我们要执行的 20 个仿真的随机贡献和我们要考虑的 2518 个时间间隔。回想一下，这些区间对应于过去 10 年的每日价格。

    使用了两个功能：

    norm.ppf()：此 SciPy 函数给出了累积概率具有给定值的变量的值。

    np.random.rand()：此 NumPy 函数计算给定形状中的随机值。它创建给定形状的数组，并使用[0，1]上均匀分布的随机样本填充该数组。

    我们将按以下方式计算每日回报：

    DailyReturns = np.exp(Drift + StdevLogReturns * SBMotion)

    每日收益率是衡量股票价格变化的一个指标。它以前一天收盘价的百分比表示。正回报意味着股票增值，而负回报意味着它失去了价值。np.exp()函数用于计算输入数组中所有元素的指数。

  4. After long preparation, we have arrived at a crucial moment. We will be able to carry out predictions based on the Monte Carlo method. The first thing to do is to recover the starting point of our simulation. Since we want to predict the trend of Amazon stock prices, we recover the first value present in the AMZN.csv file:

    StartStockPrices = AmznData.iloc[0]

    pandasiloc()函数用于使用基于位置的索引返回纯整数进行选择。然后，我们将初始化包含预测的数组：

    StockPrice = np.zeros_like(DailyReturns)

    numpy``zeros_like()函数用于返回与给定数组形状和类型相同的零数组。现在，我们将设置StockPrice数组的起始值，如下所示：

    StockPrice[0] = StartStockPrices

  5. To update the predictions of the Amazon stock prices, we will use a for loop that iterates for a number that is equal to the time intervals we are considering:

    for t in range(1, NumIntervals):
    StockPrice[t] = StockPrice[t - 1] * DailyReturns[t]

    对于更新，我们将根据以下等式使用 BS 模型：

最后，我们可以查看结果：

plt.figure(figsize=(10,5))
plt.plot(StockPrice)   
AMZNTrend = np.array(AmznData.iloc[:, 0:1])
plt.plot(AMZNTrend,'k*')   
plt.show()

打印以下图表：

Figure 8.5: Amazon trend graph

图 8.5–亚马逊趋势图

在上一张图中，黑色突出显示的曲线代表了过去 10 年亚马逊股票价格的趋势。其他曲线是我们的仿真。我们可以看到，一些曲线偏离预期曲线，而另一些曲线看起来更接近实际趋势。

投资组合管理风险模型研究

良好的风险度量在金融学中至关重要，因为它是评估金融资产的主要工具之一。这是因为它允许您监控证券，并为投资组合的构建提供了标准。多年来被广泛使用的一个衡量标准是方差。

使用差异作为风险度量

多元化投资组合在风险和预期价值方面的优势使我们能够为证券找到正确的配置。我们的目标是在相同的风险下获得最高的期望值，或者最小化获得相同期望值的风险。为了实现这一点，有必要将风险的概念追溯到一个可测量的数量，通常称为方差。因此，通过最大化每个方差水平的投资组合收益预期值，可以重建一条称为有效边界的曲线，该曲线确定了可用于构建每个风险水平的投资组合的证券所能获得的最大预期值。

最小方差投资组合表示具有最低可能方差值的投资组合，而不考虑预期值。该参数的目的是优化由投资组合方差表示的风险。仅当收益分布为正态时，仅将风险追踪到方差度量才是最优的。事实上，正态分布具有一些特性，使得方差成为足以代表风险的度量。它完全可以通过两个参数（均值和方差）来确定。因此，知道均值和方差就足以确定分布的任何其他点。

引入风险价值度量

在非正态和极限值的情况下，考虑方差作为唯一的风险度量。一个被广泛使用了二十多年的风险度量是风险价值（VaR。VaR 的诞生与金融机构管理风险并因此能够衡量风险的需求日益增长有关。这是因为金融市场的结构日益复杂。

事实上，引入这一指标并不是为了阻止方差极限作为风险指标，因为计算 VaR 值的方法正是从正态性假设开始的。然而，为了更容易理解，让我们通过对不同类型的风险采用单一指标，将证券的总体风险包含在单个数字或金融资产组合中。

在金融环境中，VaR 是给定置信区间的一种估计，即在每个时间范围内，证券或投资组合的损失可能有多高。因此，VaR 侧重于收益分布的左尾，即实现概率较低的事件所在的位置。指出损失，而不是预期价值周围收益的分散，使其更接近于风险而非方差的一般概念。

重要提示

摩根大通被认为是将 VaR 作为一个广泛衡量标准的银行。1990 年，摩根大通总裁丹尼斯·韦瑟斯通（Dennis Weatherstone）对每天收到的冗长的风险分析报告感到不满。他想要一份简单的报告，总结银行在整个交易组合中的总风险敞口。

在计算 VaR 后，我们可以说，根据置信区间给出的概率，我们在未来 N 天内的损失不会超过投资组合的 VaR。VaR 是不超过置信区间给定概率的损失水平。

例如，在 95%的置信水平下，一年 100 万欧元的 VaR 意味着在 95%的情况下，下一年投资组合的最大损失将为 100 万欧元。没有什么能告诉我们剩下的 5%的病例发生了什么。

以下图显示了投资组合收益的概率分布，并显示了 VaR 的价值：

Figure 8.6: Probability distribution of the portfolio returns

图 8.6——投资组合收益的概率分布

VaR 是以下两个参数的函数：

  • 时间范围
  • 信心水平

必须指定 VaR 的某些特征：

  • VaR 不能描述最严重的损失。

  • VaR 没有说明损失在左尾的分布情况。

  • VaR is subject to sampling errors.

    重要提示

    采样误差告诉我们采样值与实际总体值的偏差有多大。这种偏差是因为样本不能代表总体或存在扭曲。

VaR 是一种广泛使用的风险度量，用一个数字概括了金融工具组合风险的重要方面。它的计量单位与计算它所依据的投资组合的回报相同，并且很容易理解，回答了一个简单的问题：金融投资会有多糟糕？

现在让我们来看一个计算 VaR 的实际案例。

估计一些纳斯达克资产的 VaR

纳斯达克是世界上最著名的股票市场指数之一。其名称是全国证券交易商协会报价单的首字母缩写。这是代表美国科技行业股票的指数。考虑到投资者心目中的纳斯达克（NASDAQ），美国主要科技和社交机构的品牌很容易出现。想想谷歌、亚马逊、Facebook 和其他许多公司；他们都在纳斯达克上市。

在这里，我们将学习如何恢复纳斯达克上市的六家公司的报价数据，然后我们将演示如何评估与购买这些证券的股票组合相关的风险：

  1. As always, we start by importing the libraries:

    import datetime as dt
import numpy as np
import pandas_datareader.data as wb
import matplotlib.pyplot as plt
from scipy.stats import norm

    已导入以下库：

    datetime库包含用于处理日期和时间的类。它包含的函数使我们能够轻松提取用于格式化和操作日期的属性。

    熊猫图书馆的pandas_datareader.data模块包含的功能允许我们提取财务信息，而不仅仅是从一系列提供此类数据的网站中提取。收集的数据以数据帧格式返回。pandas库是一个开放源代码的 BSD 许可库，它包含用于为 Python 编程语言处理高性能数值的数据结构和操作。

  2. We will set the stocks we want to analyze by defining them with tickers. We also decide the time horizon:

    StockList = ['ADBE','CSCO','IBM','NVDA','MSFT','HPQ'] 
StartDay = dt.datetime(2019, 1, 1)
EndDay = dt.datetime(2019, 12, 31)

    一个数据帧中包含了六个代码。股票代码是一个缩写，用于唯一标识特定股票市场上特定证券在证券交易所上市的股票。它由字母、数字或两者的组合组成。这些股票用于指代全球科技行业的六家领先公司：

    ADBE：Adobe Systems Inc.——世界上最大、最具差异化的软件公司之一。

    CSCO：思科系统公司——基于互联网协议（IP）的网络和其他通信和信息技术产品的生产。

    IBM：国际商用机器——信息技术相关产品的生产和咨询。

    NVDA：Nvidia 公司-视觉计算技术。这就是发明 GPU 的公司。

    MSFT：微软公司——这是该行业最重要的公司之一，也是全球营业额最大的软件生产商之一。

    HPQ：惠普公司——全球领先的产品、技术、软件、解决方案和服务提供商，面向个人消费者和大型企业。

    在确定股票行情之后，我们设定了分析的时间范围。通过定义 2019 年全年，我们只需设置分析的开始日期和结束日期。

  3. Now we can recover the data:

    StockData =  wb.DataReader(StockList, 'yahoo',
                           StartDay,EndDay)
StockClose = StockData["Adj Close"]
print(StockClose.describe())

    为了检索数据，我们使用了pandas_datareader.data模块的DataReader()功能。此函数用于将来自各种 internet 源的数据提取到数据帧中。已通过以下主题：

    StockList：待回收存货清单

    'yahoo'：收集数据的网站

    StartDay：监测开始日期

    EndDay：监测结束日期

    恢复的数据输入到熊猫数据框中，该数据框将包含 36 列，对应于 6 种股票中每种股票的 6 条信息。每个记录将包含每天的以下信息：高值、低值、打开值、关闭值、音量和调整后的关闭。

    对于投资组合的风险评估，只有一个值足够：调整后的收盘价。此列从起始数据帧中提取并存储在StockData变量中。然后，我们使用describe()函数为每只股票开发了基本统计数据。已返回以下统计信息：

    Figure 8.7- Statistics of the portfolios

    图 8.7-投资组合统计

    分析前面的表，我们可以注意到有 252 条记录。这是 2019 年证券交易所开张的日子。让我们注意到这一点，因为这些数据以后会很有用。我们还注意到，由于股票的价值不同，列中的值具有非常不同的范围。更容易理解股票的走势，最好画图表。让我们接下来做这个：

    fig, axs = plt.subplots(3, 2)
axs[0, 0].plot(StockClose['ADBE'])
axs[0, 0].set_title('ADBE')
axs[0, 1].plot(StockClose['CSCO'])
axs[0, 1].set_title('CSCO')
axs[1, 0].plot(StockClose['IBM'])
axs[1, 0].set_title('IBM')
axs[1, 1].plot(StockClose['NVDA'])
axs[1, 1].set_title('NVDA')
axs[2, 0].plot(StockClose['MSFT'])
axs[2, 0].set_title('MSFT')
axs[2, 1].plot(StockClose['HPQ'])
axs[2, 1].set_title('HPQ')

    为了在 6 只股票的趋势之间进行简单的比较，我们追踪了 6 个子批次，这些子批次按 3 行 2 列顺序排列。我们使用了matplotlib库的subplots()函数。此函数返回包含地物对象和轴的元组。所以，当你使用fig、axs = plt.subplots()时，你将这个元组解压成fig和axs变量。如果您想在地物级别更改属性或稍后将地物保存为图像文件，则拥有fig非常有用。变量axs允许我们设置每个子批次的轴的属性。事实上，我们调用这个变量是为了通过调用图表矩阵的行-列索引来定义在每个子图中绘制什么。此外，对于每个图表，我们还打印了标题，这使我们能够了解它所指的股票代码。

    完成此操作后，我们绘制图形：

    plt.figure(figsize=(10,5))
plt.plot(StockClose)
plt.show()

    使用了以下matplotlib功能：

    figure()：此函数创建一个新的图形，目前为空，我们使用figsize参数设置框架的大小，该参数以英寸为单位设置宽度和高度。

    plot()：此函数用于绘制AmznData数据集。

    show()：此功能在 PyLab 模式下以 IPython 运行时，显示所有图形并返回 IPython 提示符。

    打印以下图表：

    Figure 8.8: Graphs of the statistics

    图 8.8–统计图

    分析前面的数字，一切都更清楚了。股票走势明显。撇开绝对值不谈，每个股票的绝对值差异很大，我们可以注意到，在整个监测期间，只有微软股票记录了几乎增加的趋势。相反，其他股票表现出波动趋势。我们还注意到，HPQ 股票记录了三次突然下跌。

  4. After taking a quick look at the trend of stocks, the time has come to evaluate the returns:

    StockReturns = StockClose.pct_change()
print(StockReturns.tail())

    pct.change()函数返回当前收盘价与上一个值之间的百分比变化。默认情况下，该函数计算前一行的百分比变化。

    时间序列百分比变化的概念与股票价格回报的概念相联系。基于回报的方法提供了数据的标准化，这对于评估以不同度量为特征的变量之间的关系具有根本重要性。这些概念已在本章使用蒙特卡罗方法进行股票价格预测一节中进行了探讨。请注意，我们只提到了其中的一些。

    然后我们打印返回的数据帧的队列，以分析其内容。返回以下结果：

    Figure 8.9: The stock returns DataFrame

    图 8.9–库存返回数据帧

    在上表中，减号表示负回报或亏损。

  5. Now we are ready to assess the investment risk of a substantial portfolio of stocks of these prestigious companies. To do this, we need to set some variables and calculate others:

    PortfolioValue = 1000000000.00
ConfidenceValue = 0.95
Mu = np.mean(StockReturns)
Sigma = np.std(StockReturns)

    首先，我们设定投资组合的价值；这是十亿美元。这些数字不应该吓到你。对于一家管理着众多投资者的银行来说，实现这一投资价值并不困难。因此，我们设置了置信区间。之前，我们说过 VaR 是基于这个值的。随后，我们开始计算 VaR 计算的一些基本量。我指的是回报的平均值和标准差。为此，我们使用了相关的numpy函数：np.mean ()和np.std。

    我们继续设置计算 VaR 所需的参数：

    WorkingDays2019 = 252.
AnnualizedMeanStockRet = MeanStockRet/WorkingDays2019
AnnualizedStdStockRet = 
                StdStockRet/np.sqrt(WorkingDays2019)

    之前，我们看到从雅虎网站的财务部分提取的数据包含 252 条记录。这是 2019 年证券交易所的工作日数，所以我们设置了这个值。那么，让我们继续对刚刚计算的平均值和标准偏差进行年化。这是因为我们要计算股票的年度风险指数。对于平均值的年化，除以工作天数就足够了，而对于标准差，我们必须除以工作天数的平方根。

  6. Now we have all the data we need to calculate the VaR:

    INPD = norm.ppf(1-ConfidenceValue,AnnualizedMeanStockRet,
                   AnnualizedStdStockRet)
VaR = PortfolioValue*INPD

    首先，我们计算置信度、均值和标准差的风险水平为 1 的逆正态概率分布。这项技术涉及从我们提到的三个参数开始构建概率分布。在本例中，我们从一些分布统计数据开始向后工作，并尝试重建起始分布。为此，我们使用 SciPy 库的norm.ppf()函数。

    norm()函数返回一个正常的连续随机变量。首字母缩略词ppf代表百分点函数，这是分位数函数的另一个名称。与随机变量概率分布相关的分位数函数规定了随机变量的值，以便变量小于或等于该值的概率等于给定概率。

    此时，VaR 的计算方法是将获得的逆正态概率分布乘以投资组合的价值。为了使获得的值更具可读性，将其四舍五入到小数点后的前两位：

    RoundVaR=np.round_(VaR,2)

    最后，打印获得的结果，每行一个，以简化比较：

    for i in range(len(StockList)):
    print("Value-at-Risk for", StockList[i], 
                   "is equal to ",RoundVaR[i])

    返回以下结果：

    Value-at-Risk for ADBE is equal to  -1547.29
Value-at-Risk for CSCO is equal to  -1590.31
Value-at-Risk for IBM is equal to  -2047.22
Value-at-Risk for NVDA is equal to  -1333.65
Value-at-Risk for MSFT is equal to  -1286.01
Value-at-Risk for HPQ is equal to  -2637.71

回报风险最高的股票是惠普和 IBM，而回报风险最低的股票是微软股票。

总结

在本章中，我们应用了基于蒙特卡罗方法的仿真概念，更一般地，将随机数的生成应用到与金融工程领域相关的实际案例中。我们首先定义了基于布朗运动的模型，该模型描述了小颗粒在流体中不间断和不规则的运动。我们学习了如何描述数学模型，然后导出了一个实际应用程序，将随机行走仿真为维纳过程。

随后，我们讨论了另一个相当有趣的实际案例，即如何使用蒙特卡罗方法预测著名亚马逊公司的股票价格。在过去 10 年中，我们开始探索亚马逊共享的趋势，并进行简单的统计，以提取我们通过视觉分析确认的任何趋势的初步信息。随后，我们学会了将股票价格趋势视为时间序列，计算每日回报。然后，我们用 BS 模型解决了这个问题，定义了漂移和标准布朗运动的概念。最后，我们应用蒙特卡罗方法预测与股票价格趋势相关的可能情景。

作为最后一个实际应用，我们评估了一些在纳斯达克市场上市的最著名科技公司的股票组合的风险。我们首先定义了与金融资产相关的转介概念，然后引入了 VaR 的概念。随后，我们实现了一种算法，在给定置信区间和时间范围的情况下，根据股票价格历史数据返回的每日回报计算 VaR。

在下一章中，我们将学习人工神经网络的基本概念，如何将前馈神经网络方法应用于数据，以及神经网络算法如何工作。然后，我们将了解深层神经网络的基本概念，以及如何使用神经网络仿真物理现象。*

avatar
文章
阅读
粉丝