大家好,我叫戴尔·休曼 (Dale Human),是通用电气医疗集团 (GE Healthcare) 的一名工程师,工作地点在威斯康星州密尔沃基市。今天,我将和大家简单聊聊如何构建模型,如何使用自助法 (bootstrap) 来评估预期结果,以及模型预测结果的偏差。如果想更深入地了解模型,还可以使用蒙特卡罗模拟,根据输入变量的变化来进一步评估系统预期偏差。内容很多,我们开始吧。首先,我先简单介绍一下我的兴趣所在。
这一切都归结于新产品开发或研发等情况下,当你的样本量很小时,你只能进行少量测试,并试图根据这些测试或测量结果来决定下一步的计划或系统可能的表现。在这种情况下,像引导法和蒙特卡罗模拟这样的工具就非常适用。让我们来看看如何在 JMP 中实现这些操作。
我们先从一个相当简单的系统入手。假设你有一个夹紧系统,这里有一把螺丝刀,用来拧紧某个组件。我们在组件中安装了一个测力传感器,用来测量拧紧螺丝时组件所受的力。这是一个非常典型的机械应用。我们想知道,如果我测量螺丝上的扭矩,并由此得到测得的力,我们能否建立起扭矩与力之间的关系,并尝试模拟该系统未来的运行情况?
我这里有一个示例数据表,我会在这里打开它。您可以看到,我们做了一些测试,在这个例子中是九次测试,我们对螺杆驱动器施加了一定的扭矩,并测量了力。我们希望了解力和施加的扭矩之间的关系。让我们建立一个模型,因为我们只有一个响应变量和一个输入变量(扭矩),所以我们可以使用“拟合 Y 对 X”功能。
我将进入“分析”菜单,然后选择“Y 对 X 拟合”。我们可以将测量力作为响应变量,将施加的扭矩作为输入变量输入到对话框中。JMP 提示说,由于只有两个变量,它将进行双变量分析。我点击“确定”。和 JMP 一样,我们首先得到一个图表,一个数据散点图。我们可以直观地观察施加的扭矩和测量力之间的关系。正如您在这里看到的,它们之间似乎存在某种关系。随着扭矩的增加,测量力也趋于增加。这很合理。
看起来随着扭矩越来越大,这里似乎开始出现一个平台期。在这个系统中,你可以想象一下,当扭矩达到一定程度时,组件会被完全夹紧,因此继续施加更大的扭矩从力学角度来看并没有什么帮助,但我们确实看到了这种关系。我们来尝试建立一个模型。我将选择这里的红色三角形。你可以使用“拟合直线”来拟合线性关系,但如果我认为存在这种平台期效应,我将使用“特殊拟合”对话框,它允许我选择 X 轴和 Y 轴的一些变换。
在这种情况下,我将选择 X 的自然对数,这有助于捕捉该系统可能存在平台效应的趋势,所以我选择自然对数。我还关闭了中心多项式选项。您可以保留它,它不会改变任何结果,只是看待方程的另一种方式。目前,我已将其关闭。我只需点击“确定”。我们可以看到,这是我们的最佳拟合线。这本质上是使用分析中施加扭矩的自然对数进行最小二乘回归得到的,这是我们的最佳拟合曲线。它似乎与我们的数据吻合得很好。
很好。我们已经有了根据分析结果建立的方程,所以我们可以看到这一点。如果您之前做过回归分析,您可能已经看过 R 平方值。我们这里的 R 平方值看起来非常好,都远高于 90%。调整后的 R 平方值和原始 R 平方值非常接近。在这种情况下,只有一个 X 变量,这也在意料之中。
如果你看一下均方根误差(RMSE),它实际上类似于模型的标准偏差。它可能有点大,可能不太容易看出来。我接下来会查看我们应用的拟合结果,并通过点击“置信区间”选项来启用该拟合的置信区间。我们可以看到模型的不确定性。例如,假设我施加了 100 牛米的扭矩,那么我测量的平均力的预期变化范围大约在 1000 牛米左右。这就是置信区间的宽度。
这个数值有点大。您可能在想,只有九个数据点,所以这个模型存在很大的不确定性,这很可能是事实。我们下面还有参数估计分析,其中显示了模型参数的值,包括方程中的斜率和截距值。此外,您还可以看到这些值的标准误差,它量化了这些值的不确定性。
您可以看到,截距的标准误差约为 2000 牛顿,斜率项的不确定度约为每牛顿 500 牛顿米。我们正在尝试理解,未来是否真的会预期出现如此大的波动。如果我们收集越来越多的数据,不确定度真的会如此之大吗?
鉴于我们构建这个模型时仅使用了很小的样本(本例中只有九个数据点),我们该如何评估预期结果的真实性?或者说,由于样本量小,我们是否高估了模型的误差?一种方法是使用一种叫做自助法(Bootstrap)的技术。
我们来谈谈自助法(Bootstrap)。自助法其实是一种“简单”的技术,它通过对原始数据进行多次重采样,来尝试重新估计待计算参数的值。在这个例子中,我们会多次重新估计模型的斜率和截距,以便更好地了解该系统未来可能存在的真正不确定性。这就是自助法的基本思想。如图所示,我们用三个数据点作为样本。
自助法(bootstrapping)的作用是从原始数据集中随机抽取一个样本。这里,我们在这个方框里放了一个随机样本,它包含三个数据点。不过,这次我是随机选择的,碰巧抽取了第3行两次,第1行一次。第2行没有出现。这个过程需要重复很多次,成百上千次。每次抽取新样本后,都需要重新计算系统中感兴趣的参数。在本例中,我们最终需要查看模型的斜率和截距值。
无论你感兴趣的参数是什么,你都要对每个不同的样本计算该参数,并且重复这个过程很多很多次,以帮助理解预期的变化范围,即使我的数据集很小。这就是自助法的基本思想。它是一种常用的技术。它的优点在于它是非参数的。它不依赖于正态性或其他类似的假设,因此适用于许多不同的情况。
当你的样本量很小时,就可以使用这种方法。之所以会想到它,是因为你试图估计模型或分析计算中的误差。通常,你会使用方差之类的指标来观察系统的变异性。记住,方差的计算方法是除以 n 或 n 减 1。如果样本量很小,那么除数通常也很小,因此你往往会高估预期的变异性,但如果继续收集数据,这种预期可能并不现实。自助法(Bootstrap)就是一种尝试在这种情况下估计更实际变异性的方法。
这就是自举法的基本概念。举个简单的例子,在我们回到之前建立的模型之前,我们可以先看一下这里的数值分布。我们同样有九个测量力值。让我们先对它们进行分布分析。我将进入“分析分布”菜单,在这个例子中,我们先只查看力的分布。让我们看一下数值分布。
我这里有数据。你可以看到一个汇总统计表,我暂时隐藏了分位数。这并不是我关注的重点,但如果你看一下,你会发现这个例子中平均值约为 7100 牛顿,标准差约为 2000 牛顿。标准差相当大,因为我们的数据样本量很小。这是由于样本量 n 较小造成的。
我们能否确定这是否真的是我们预期未来总体变异的趋势?考虑到我们只测量了总体中的九个样本,我们可以尝试使用自助法来判断这是否是总体变异的合理预期。具体该如何操作呢?JMP 软件让这一切变得非常简单。我们已经有了包含目标参数的表格,我们感兴趣的是均值和标准差。这两个值都在表格中,所以我只需右键单击该表格。你会看到之前可能见过的所有菜单项,但在最底部是“自助法”。
这样可以让你多次对数据进行重采样。对于每个单独的样本,它都会重新计算所有统计数据。然后你可以查看这些值的分布,以便更好地了解你的系统在样本量越来越大的情况下可能会如何运行。默认的样本数为 2500,因此 JMP 会查看我测量的力值,并从这九个初始值中随机抽取九个。同样,有时它会多次抽取同一个值,就像我们在图中看到的那样。
你不可能每次都得到完全相同的九个数据点。这就是自助法的优势所在。默认情况下,它会运行 2500 次。我点击“确定”,让它运行。这只需要几秒钟。这是一个非常简单的计算。但请注意,如果你的模型非常复杂,那么进行大量计算显然会花费一些时间。如果你要计算的模型计算量非常大,那么……在这个例子中,我们只是计算均值和标准差。这很简单。现在,JMP 创建了一个新的数据表。
你会看到它有 2501 行。因为它对我们的原始数据进行了 2500 次重采样。第一行就是原始数据本身的值,也就是第 0 行。然后所有其他行都是针对此表中每个统计量的重采样计算值。如果我把这些值单独列出来,由于所有这些值都在同一个表中,JMP 会尝试对每个统计量进行引导抽样。目前我们并不关心其他所有值。我们真正想关注的是均值,因为我真正想了解的是系统平均性能(在本例中为平均力)的预期变化范围。
我将重点关注均值。JMP 会自动为表格中的每个变量构建分布分析。在这里,我们可以运行脚本,它会打开所有这些不同变量的分析。现在,我们真正想关注的是均值。我将按住 Ctrl 键并隐藏所有这些分布分析。然后只显示均值,这样屏幕会更清晰,方便大家查看。我将把它拖到分布分析旁边。
我打算再次隐藏分位数。我并不关注分位数。也许你会对分位数感兴趣,但我们来看看这里计算出的值。这里有2500个不同的均值,它们是通过对这个样本进行重采样得到的。我看均值大约是7100牛顿,非常接近我们最初的值。然而,标准差却小得多。我们最初的标准差大约是2100牛顿。
现在数值只有大约 670 牛顿,这是因为我们从这个庞大的总体中进行了多次重采样。我们对标准差的估计可能越来越接近预期值,也就是这个总体中真实的变异范围。这里我们夸大了变异范围,是因为样本量较小。我们是从同一个总体中重采样的。所有这些计算中使用的所有数据点都来自与这九个数据点相同的总体,因为它们都取自相同的数据。
这或许更能真实地反映我们预期在这个群体中看到的变异情况。这就是自助法的威力。现在你可能会想,如果我抽取一个更大的自助样本,我的标准差会继续缩小吗?我们来试试。我回到原始数据。在这个表格中右键单击,选择“自助法”。这次我们抽取 10,000 个自助样本。我点击“确定”。
运行还需要几秒钟,所以我们让它运行。可能不止几秒钟。现在我们有一个包含 10,000 行的数据表,同样,每一行都是同一总体的不同样本。我们可以查看这个表中的分布情况。我把它打开。我仍然只关注均值。我关闭所有分布窗口,只打开均值窗口。我们把它放在其他窗口旁边。抱歉,我稍微调整一下窗口位置,这样我们就可以把它们放在一起看了。
现在我们有了一个包含 10,000 个值的自助抽样样本。我们可以看到均值仍然在 7100 左右,但标准差仍然在 670 左右。这其实并不奇怪,因为所有这些值都是从同一个总体中计算出来的。这个总体有一定的均值、一定的变异性和一定的标准差。我们不应该期望这些值会发生变化,因为这两项分析的样本量都非常大。
这与之前的研究不同,因为这次的样本量非常小。但这确实能让你看到未来的发展趋势。你可能会预期变异范围更接近 670 牛顿,而不是 2000 牛顿。当然,图表底部也计算了一些置信区间,不同的置信水平可以让你了解总体预期变异范围。
你会再次注意到,自助法样本的值彼此非常接近,而且比我们小样本均值的原始置信区间略窄。这是因为我们使用了可能更符合实际的标准差估计值。这就是自助法的强大之处。它能让你了解系统中的变异情况,即使你最初只有一小份数据样本。好了,现在我们已经基本了解了自助法的作用,我将关闭这些表格,以免它们妨碍我们继续讨论。
回到我们刚才提到的日志,我已经讲解过如何进行自助法,并且已经演示过了。你的标准差和均值实际上应该不会改变,因为它们代表了最初提取数据的总体。这就是自助法的原理。让我们回到模型。我打开那部分,我们可以用同样的方法处理这个模型。在这个模型中,我们计算了最小二乘回归模型的斜率和截距参数的估计值。问题是,我真的认为这些估计值的误差应该这么大吗?
我们可以像之前计算这些力的均值那样,对这些估计值进行自助法抽样。同样,你只需要右键单击表格中的“估计值”列,选择“自助法”,然后选择抽样次数。这次我选择 2500,这样速度会快一些。计算需要几秒钟时间。现在我们又得到了一个包含 2500 个值的表格。接下来,我将查看 JMP 生成的分布分析。由于我们只有两个估计值(一个用于截距,一个用于斜率),因此该分析中只有两个值。
我再次按住 Ctrl 键,同时关闭这两个窗口,即可隐藏分位数。这里截距的标准误差约为 2000 牛顿,斜率的标准误差约为 450 牛顿。那么,我们从自助法中得到了什么结果呢?如果我们查看自助样本的标准差,它仍然约为 2000,在本例中为 2200。类似地,斜率自助法的标准差约为 475,略有不同,但接近最小二乘回归分析得到的原始值 452。你可能会想,等等,自助法并没有起到作用,你说得对。在这种情况下,自助法并没有起到太大作用。
这其实是意料之中的,因为如果你还记得最小二乘回归,它的目的是选择能够最小化模型方差的拟合优度。斜率和截距的估计值之所以被选中,正是因为它们能够最小化模型误差。对于最小二乘回归来说,自助法的作用并不大。事实证明确实如此。我只是想用这个例子来说明如何将自助法应用于模型。
现在,如果您使用其他无法提供误差估计的建模技术,例如神经网络分析或主成分分析,而分析结果无法直接告诉您预期变异是多少,也就是标准误差是多少,那么您可以使用自助法来帮助了解模型中的误差。但在这里,我要向您展示的是,自助法的作用并不大。它通常与最小二乘分析的结果一致,因为最小二乘分析旨在最小化这种变异。
这个概念适用于任何建模技术。再次强调,因为自助法并非参数方法,它不对数据的分布做任何假设。因此,它可以应用于你可能使用的所有不同类型的建模技术。例如,你可以使用建模(或者更准确地说是自助法)来帮助理解最小二乘回归的误差。但这并非特别有用。如果你只是想了解均值的变化,就像我们在这里处理力的情况一样,那么你可以使用自助法来更真实地了解未来的变化趋势,甚至在进行任何更详细的建模之前就可以这样做。你也可以这样使用它。
这就是自助法的思想。我已经用最简单的例子展示了它,即观察均值的预期变异。此外,你还可以用它来理解模型估计的误差。它不仅适用于最小二乘法,也适用于其他一些方法。在这方面,它非常有用。
另一种评估系统变化影响的方法是蒙特卡罗模拟,JMP 也提供了非常简便的蒙特卡罗模拟功能,它是建模平台的一部分。我们来看一下蒙特卡罗模拟。从根本上讲,蒙特卡罗模拟的作用是什么?假设你有一个模型,或者更确切地说,你有一个方程。这里,我们根据分析结果建立了一个模型方程,这里使用的是同一个模型。蒙特卡罗模拟的作用是,将变化注入到该模型的输入中,并评估这种变化对输出因子(即响应变量)的影响。
你的做法是假设输入存在某种变化。例如,我假设施加的扭矩服从正态分布,中心值是一个感兴趣的值,标准差为某个值。现在,我从这个分布中随机抽取大量值,然后将这些值代入方程并计算输出,以便了解输入的变化情况。以下是我预期输出的结果。这就是蒙特卡罗模拟的基本原理。
同样,通常需要对大量样本进行此操作,以帮助您了解响应中预期会出现哪些变化。JMP 还允许您在输出中注入其他类型的变化。这里我就不赘述了,您只需知道这是一个可用的选项。我们如何使用 JMP 进行蒙特卡罗模拟呢?非常简单。首先,您需要建立一个模型。在本例中,我们构建了一个最小二乘模型。接下来,我将进入模型拟合界面并启用分析器。这就是您的起点。如果您之前见过分析器,它只是根据您之前的分析结果确定的模型方程的图形化表示,您可以调整输入值 X,并观察 Y 的变化。
您还可以通过性能分析器进行蒙特卡罗模拟。点击性能分析器旁边的红色三角形,即可启用模拟器。这就是蒙特卡罗模拟器。它允许您指定 X 值的变化方式,并显示响应变量在您进行此操作时将发生的变化。在这里,我可以从下拉菜单中选择 X 的变化方式,在本例中,X 代表施加的扭矩。我选择“随机”。选择“随机”后,您可以指定要用于该输入变化的分布。
例如,假设我们希望将该系统的扭矩设置为 125 牛米,这是我们的典型扭矩值。您已经使用即将使用的扭矩扳手进行了一些快速测量,并确定如果我将其编程为 125 牛米,则扭矩通常为 125 牛米,但标准偏差约为 3 牛米。您可以输入您预期值的平均值和标准偏差。
然后你只需要点击这里的“模拟”按钮,JMP 就会计算一个很大的数值,这个数值就是你在这个框里输入的数值。在这个例子中,它计算了 10,000 个值。JMP 从均值为 125、标准差为 3 的正态分布中随机抽取了 10,000 个不同的值。它使用之前确定的参数估计值,将这 10,000 个值代入你的模型方程,并计算所有这些变量的响应值。这个分布图展示了结果的样子。
在这个视图中,图表看起来有点折叠,只是因为图表的比例尺比较窄,您可以根据需要放大。您也可以让 JMP 将这些值显示在下面的表格中,这里是“模拟到表格”选项。我将打开它并单击“生成表格”按钮。这里显示的是从正态分布中随机选择的 10,000 个不同的施加扭矩值。这些值被代入公式,然后用于计算力。
JMP 已经为您生成了一份分布分析报告,并保存在这里。现在我打开它。这里显示的是测量力的分布情况。这里展示的是更高分辨率的分布图,以便您更清楚地了解系统中的运行情况。这份测量力分布图展示了根据预设分布,输入扭矩变化对系统响应的影响。
现在,您可以进行各种不同的分析,例如过程能力分析、公差区间分析等等。所有这些分析都可以在此时进行,以帮助您了解在输入变量存在这种变化的情况下,系统响应会如何变化。具体选择哪种分析方法取决于您自己。希望您可以使用一些历史数据。您已经查看了用于该系统的扭矩扳手,并且已经注意到,扭矩确实会根据这种分布而变化,因此您可以选择合适的分布和参数。
你可能拥有领域专业知识,能够思考系统中的运行机制,并得出结论:这些值应该服从均匀分布;当然,你也可以进行猜测,因为这是一个模型。如果你不确定,可以先假设服从某种正态分布,看看会发生什么。然后,你可以进行许多你可能感兴趣的假设分析。这就是蒙特卡罗模拟。在 JMP 中操作非常简单。蒙特卡罗模拟和自举法都允许你评估系统中的变异性:自举法基于小样本,而蒙特卡罗模拟则基于已知的输入变异性来评估输出的变化。
我们今天讨论了几个不同的方面。再次强调,自助法是一种可以用来了解系统变异情况的技术,即使样本量很小。而蒙特卡罗分析则允许你将变异注入模型的输入,从而了解这会对模型的输出产生什么影响。这些都是非常常用、非常流行且功能强大的技术。一旦你建立了一个模型,这些工具就非常实用。这就是我今天要讲的内容。非常感谢。