原文：Testing and Tuning Market Trading Systems Algorithms in C++
一、介绍
在我们深入研究这本书的内容(或者豆腐，如果你喜欢的话)之前，我们应该清楚在这里你会发现什么，不会发现什么，以及读者应该做什么样的准备。
目标受众和内容概述
这本书的目标读者是那些统计学背景一般的人(统计学 101 已经足够了)，对任何语言都有一些编程技巧的人(这里的例子中使用了非常倾向于传统 C 语言的 C++)，以及对金融市场交易感兴趣的人，他们的数学严谨程度远远超过大多数交易者。在这里，你会发现一个有用的算法集合，包括样本代码，这将帮助你调整你的想法到交易系统，有高于平均水平的盈利能力。但是有很多东西你在这本书里找不到。我们从概述这本书的内容开始。
这本书里有什么
本书涵盖了以下主题:
-
如果你的系统涉及到参数的优化，大多数都是这样，你将学会如何确定你的优化系统是否捕捉到了真实的市场模式，或者它只是学习了随机的噪声模式，而这些模式不会再出现。
-
您将学习如何修改线性回归，使其更不容易过度拟合，并且作为一个奖励，将预测器分为有价值的和无价值的。您还将学习如何修改线性回归，以便在适度非线性的情况下使用它。
-
你会发现一个非常通用和强大的非线性优化算法，适用于基于预测模型的交易系统和传统的算法系统。
-
所有的交易系统都假定交易市场有一定程度的一致性；如果你的系统所基于的模式在最近的历史中经常发生，我们必须假设同样的模式至少会持续到不久的将来。一些交易系统对市场模式的适度变化很稳健，而另一些系统会因为市场模式的微小变化而变得一文不值。您将学习如何评估您的系统对这种变化的稳健程度。
-
如果你已经设计了自己的专有指标，你将学会如何确认它们是合理稳定的(这是任何有效指标的关键属性),或者如果它们不稳定，如何将它们调整到平稳状态。你还将学习如何计算它们，以最大化它们的信息含量，最小化它们的噪音，并以有效的方式提供给你的交易系统，以最大化它们的效用。
-
大多数交易系统开发人员都熟悉 walkforward 测试。但没有那么多人意识到，普通的前向算法通常不足以正确验证交易系统候选，而且可能因为微妙的原因产生危险的乐观结果。您将了解如何在第二层 walkforward 分析中嵌入一个 walkforward 算法，或者在 walkforward 分析中嵌入一个交叉验证层。这种“验证中验证”的场景不仅是测试交易系统的最佳方式，也是唯一真正正确的方式。
-
你将学习如何估计你的系统可能产生的未来利润的范围。如果你发现你的系统有几乎确定的未来盈利能力，但是这种盈利能力相对于产生的风险来说很有可能很小，你就知道你的系统还没有准备好交易。
-
你将学习如何估计灾难性下降的概率，即使你的系统运行“正确”
-
您将了解严格的统计测试算法，这些算法能够抵抗偶尔的大赢大输，这些算法会使许多“传统”验证算法无效。
-
许多交易系统开发人员喜欢使用“墙上的意大利面条”的方法来开发交易系统。尽管经常遭到蔑视，但这实际上是一种合法的方法，只要明智地去做。您将学习如何确定众多竞争系统中的“最佳”是否真正值得。
这本书里没有什么
本书不包括以下主题:
-
这本书不是“市场交易者统计学入门”类型的书。假设读者已经熟悉了均值和标准差、正态分布、假设检验的 p 值等概念。不需要比这些概念更高级的东西；这里介绍的高级统计技术是建立在基本概念的基础上的，任何通过了统计学 101 甚至是心理学统计学课程的人都可以掌握。但是如果你不知道什么是标准差，你会发现这本书很难懂。
-
这也不是一本“金融市场交易入门”的书。假设你知道一些术语的含义，比如开始交易的和结束交易的，多头和空头，表示每笔交易的收益。如果你完全是金融市场交易的新手，在阅读这本书之前，你需要学习背景材料。**
** 在这里，你会发现很少或没有实际的、经过验证的交易系统。那些一毛钱一打，通常值这个价。但是如果你对交易系统有自己的想法，你将学习如何实现、测试和调整它，以最大化它的利润潜力。
* 你在这本书里找不到绝密的超级骗子保险指标。提出的几个指标要么是常识，要么在公共领域广为流传。但是如果你对指标有自己的想法，你就会学会如何最大化它们的效用。*
*## 关于交易系统
由于本文中介绍了不同的测试程序，它们必然会在不同的交易系统中进行演示。请注意以下感兴趣的项目:
-
我并不赞同这些系统赚钱。相反，我尽可能地保持系统的简单，这样就可以把重点放在它们的测试上，而不是它们的实际效用上。这本书假设读者对交易系统有自己的想法；这里的目标是为调整和严格测试现有系统提供先进的统计方法。
-
所有用于演示的交易系统都假设我们在使用日线，但这不是必须的。条形可以是任何长度，从几分之一秒到几个月。事实上，大多数演示只使用每根棒线的开盘价或收盘价，所以将这些算法应用于交易分笔成交点数据也是可行的。日柱是最方便的，测试数据也是最容易得到的。
-
大多数演示系统在酒吧收盘时开始和结束交易。自然，在现实生活中这是困难或不可能的；更公平和保守的方法是在一根棒线收盘时做交易决定，在下一根棒线开盘时开始或结束交易。但是这会给这里显示的算法增加不必要的混乱。请记住，我们的目标是以最直接的方式呈现统计算法，将焦点放在统计测试上。在大多数情况下，对实现的小修改不会从本质上改变严格统计测试的结果。
-
在这些测试中，交易成本(滑点和佣金)被故意忽略了，这也是为了保持对统计测试的关注，而不增加混乱。提供的代码和附带的描述清楚地说明了如果需要，如何将交易成本纳入计算中。
市场价格和回报
大多数股票市场的价格范围很广，或许开始时每股交易几美元，经过分割调整后，今天每股交易数百或数千美元。当我们计算交易的回报时，我们不敢只减去开盘价和收盘价。从 1 美元到 2 美元的 1 美元波动是巨大的，而从 150 美元到 151 美元的波动几乎是微不足道的。因此，许多人计算百分比变动，将价格变动除以起始价格，然后乘以 100。这就解决了尺度问题，而且很直观。不幸的是，它有一个问题，这使得它在许多统计分析中是一个糟糕的方法。
百分比移动的问题在于它们是不对称的。如果我们在一笔交易中赚了 10 %,然后在下一笔交易中亏损了 10 %,我们就没有回到起点。如果我们从 100 点涨到 110 点，但是损失了 110 点的 10 %,我们就是 99 点。这似乎并不严重，但是如果我们从不同的角度来看，我们就会明白为什么它会成为一个大问题。假设我们有一个长期交易，市场从 100 移动到 110，我们的下一个交易从 110 回到 100。我们的净权益变化为零。然而，我们记录了 10%的收益，然后是 9.1%的损失，净收益接近 1%。如果我们记录一系列的交易回报进行统计分析，这些错误会很快增加，结果是一个完全没有价值的交易系统可以显示出惊人的净收益！这将使几乎所有的性能测试无效。
有一个简单的解决方案被专业开发人员使用，我将在本书中通篇使用:将所有价格转换为价格的对数，并将交易回报计算为这些对数的差值。这解决了所有的问题。例如，捕捉到从 10 到 11 的市场波动的交易是 2.39789–2.30258 = 0.09531，捕捉到从 100 到 110 的波动的交易是 4.70048–4.60517 = 0.09531。如果一笔交易让我们从 110 回到 100，我们损失 0.09531，净收益为零。太好了。
这种方法的一个好处是，较小的对数价格变化乘以 100，几乎等于百分比变化。例如，从 100 到 101，1%的变化，相当于 100 *(4.61512–4.605)= 0.995。甚至早些时候提到的 10%的变动也相当于 9.531%。出于这个原因，我们将把从日志中计算出的回报视为近似的百分比回报。
两种类型的自动交易系统
最初，所有形式的自动市场交易都是所谓的算法或基于规则的 *。*系统开发人员提出了一套严格定义的规则来指导开仓和平仓。规则可能会规定，如果某个条件组合变为真，那么您将打开一个多头头寸并持有该头寸，直到某个其他条件组合变为真。算法交易的一个经典例子是均线交叉系统。计算短期和长期均线，如果短期均线高于长期均线，就做多，否则就做空。训练这个原始的交易系统是通过寻找在历史数据集上提供最佳性能的短期和长期回顾来执行的。算法系统，许多涉及几十个条件，今天仍在广泛使用。
最近，许多开发人员(包括我自己)已经形成了这样的观点，即基于模型的系统更强大，尽管它们有一个共同的缺点，即它们经常盲目信任内部工作方式深不可测的黑盒。在基于模型的自动交易中，我们计算一个或多个(通常更多)?? 指标，这些指标是回顾过去并衡量市场特征的变量。这些可能包括趋势、波动性、短期循环行为等等。我们还计算了一个展望未来并描述近期市场行为的目标变量。目标可以是下一根或几根棒线的大小和市场运动方向。目标也可能是一个二元标志，它告诉我们市场在触及保护性止损点之前是否先触及了一个预定义的利润目标。然后，在给定指标变量的值的情况下，我们训练一个预测模型来估计目标变量的值。为了交易该系统，我们用指标的当前值呈现训练的模型，并考虑模型的预测。如果预测足够强(表明有信心)，我们根据预测的走势建立市场头寸。
基于模型的交易相对于基于规则的算法交易的优势在于，我们可以利用人工智能领域的许多最新发展，让运行在强大计算机上的复杂程序发现交易系统，这些系统可能非常复杂或模糊，没有人可能希望发现并以明确的规则编程。当然，这需要付出很高的代价:我们通常不知道模型到底发现了什么“规则”，我们必须盲目地接受模型的决定。
因为这两种类型的交易系统开发在今天都被广泛使用，所以本文将迎合这两种思想流派。不可避免的是，这里给出的一些统计检验只适用于其中一个。但是总是会尝试设计测试过程，让从业者可以使用这两种风格。
相信电脑的痛苦
对许多人来说，尤其是经验丰富的凭直觉交易者，走向自动化交易最困难的部分是当计算机的交易决策与他们的直觉冲突时接受它们，更不用说他们多年成功的交易了。我将从我自己的亲身经历中给出一个具体的例子。我根据合同开发了一个短期日内交易系统。我对该系统进行了极其彻底、严格的统计测试，结果明确显示，当该系统在冒着偶尔出现巨额亏损的风险(一种非常宽松的保护性止损)的情况下获取大量小额利润时，其利润最大化。这惹恼了负责将信号交易传到场内的交易员。他不断地用他的口头禅“减少你的损失，让你的胜利奔跑”来轰炸我。这是一些交易风格的真理，但不是这个特殊的系统。他控制不住自己；他一直否决电脑的交易决定。系统会要求平仓赢利的交易，但他会保持平仓，希望有更大的收益。或者，市场会朝着一个未平仓的位置移动，他会在系统止损点到达之前很久就把它平仓。他一直告诉我，如果他让它继续下滑，而不是早点止损，会损失多少钱。并行运行的计算机模拟赚了很多钱，而他的修改版本赚得少得多，这一事实对他的观点没有影响。多年来，他一直是一个成功的全权交易者，他知道如何交易，没有#$%^的电脑会告诉他别的方法。我们的关系从未成功过。这个故事的寓意是:如果你没有勇气相信自动化交易，那就忘掉它吧。
未来的泄露比你想象的更危险
未来泄露是未来知识非法泄露到测试程序中。在交易系统的开发和测试中，当未来市场行为的某个方面模拟了交易系统在现实生活中的表现时，就会发生这种情况。因为当我们交易我们的系统时，我们显然不知道未来，这种泄漏导致乐观的性能估计。
我不止一次地惊讶于原本严肃的系统开发人员如此随意地采取这种形式的欺骗。我有过聪明的、受过教育的开发人员耐心地向我解释说，是的，他们确实理解一些小程度的未来知识参与了他们的性能模拟。但是他们又煞费苦心地解释这种“不可避免的”泄漏是如此之小，以至于微不足道，不可能对他们的结果产生任何实质性的影响。他们不知道。这就是为什么本文反复强调的重点是避免未来泄漏的方法，即使是最微小的接触。在我早期的系统开发中，我经常惊讶于这种泄漏是如此的微妙。
为了强调这一点，图 1-1 显示了一个几乎随机的赢 1/输 1 交易系统的股票曲线，只有 1%的优势。这条曲线，如果它真的是随机的(无价值的)，平均来说会是平坦的，仅仅从这个微小的边缘来看是相当可观的。未来的泄露远比你想象的更致命。认真对待。
图 1-1
有 1%边的随机系统的权益曲线
百分比获胜谬论
有一个简单的数学公式，对交易系统的开发和评估至关重要，但对许多人来说似乎很难接受，即使他们在智力上理解它。见方程式 1-1 。
)
(1-1)
这个公式表示，交易的预期回报(如果这种情况重复多次，我们将获得的平均回报)等于我们将赢的金额乘以赢的概率减去我们将输的金额乘以我们将输的概率。
很容易接受的是，如果我们掷一枚公平的硬币，正面赢一美元，反面输一美元，我们的预期收益是零；如果我们重复抛硬币很多次，从长期来看，每次抛硬币的平均回报是零。同样容易接受的是，如果硬币是公平的，我们赢了两美元，但只输了一美元，我们就处于令人羡慕的地位。
现在考虑交易一个真正随机游走的市场；在其他属性中，从一个条形到下一个条形的变化都是相互独立的，并且平均值为零。开发一个除了零期望(当然忽略交易成本)之外的交易系统是不可能的。但是我们可以很容易地改变预期的输赢，以及输赢的频率。
例如，假设我们建立了一个多头头寸，并在进场价格之上设置了 1 个点的利润目标，在进场价格之下设置了 9 个点的止损出场。每一次我们经历的损失，都是巨大的，是我们赢的 9 倍。但是如果我们在假设的随机市场上执行大量这样的交易，我们会发现我们赢的次数是输的 9 倍。我们赢了 90%的机会。根据等式 1-1 ，我们每笔交易的预期回报仍然是零。这里的要点是，输赢的大小和概率是密不可分的。如果有人吹嘘他们的交易系统多长时间赢一次，问问他们的盈亏情况。如果他们吹嘘他们的赢和他们的输相比是多么巨大，问问他们多久赢一次。二者都不是孤立存在的。*
二、预优化问题
评估和改善平稳性
从本质上来说，时间序列的平稳性(如市场价格变化、指标或单个交易收益)是指其统计属性随时间保持不变的程度。统计学家可能会对如此宽松的定义感到畏缩，但这抓住了这个术语的实际意义。当我们使用市场历史来创建一个(最好是)有利可图的交易系统时，我们隐含地指望产生回溯测试盈利能力的历史模式至少在不久的将来仍然有效。如果我们不愿意做出这样的假设，我们还不如放弃交易系统的设计。
这个概念有许多方面与金融市场的自动交易特别相关。
-
市场，以及从市场历史中得出的指标和交易回报，本质上是不稳定的。它们的属性不断变化。唯一的问题是:情况有多糟？我们能处理吗？我们能解决问题让它变得更好吗？
-
对非平稳性进行任何严格的传统统计检验都是没有意义的。事实上，我们进行的任何测试都会显示出非常显著的统计非平稳性，所以我们不必费心；我们已经知道答案了。
-
非平稳性可以有无数种形式。也许方差在一段时间内是相当恒定的，而平均值则是漂移的。反之亦然。或者偏斜度可能改变。或者…
-
一些类型的不稳定可能对我们无害，而另一些可能对我们的交易系统是毁灭性的。一个交易系统可能有一种不稳定的弱点，而另一个交易系统可能被不同的东西拖累。在评估平稳性时，我们必须尽可能地考虑上下文。
-
评估一个完成的交易系统的耐用性的最好方法是使用第 142 页给出的渐进向前算法。
但是我们将忽略最后一点。这一章专门讨论在交易系统的发展过程中，我们应该在之前考虑的问题。渐进向前行走出现在开发的最后，是几个最终验证过程中的一个。
排除了非平稳性的传统统计测试，那么你应该怎么做？你绝对必须仔细研究你的指标。你可能会对你所看到的感到惊讶。它们的中心趋势可能会慢慢上下波动，使得预测模型在一个或两个极端都没有用。日复一日的徘徊是正常的，但缓慢的徘徊，或方差的缓慢变化，是一个严重的问题。如果一个指标在回到更“正常”的行为之前花了几个月甚至几年的时间，模型可能会在这些延长的时间内关闭或做出错误的预测。我们必须警惕这种灾难性的情况，如果我们不小心的话，这种情况很容易发生。
有时我们可能没有指标来绘制。下一节中显示的 STATN 程序是一个有价值的替代方案。但是理解非平稳性的潜在问题是很重要的。设计一个年复一年运行良好的自动交易系统，不需要调整，甚至完全重新设计，是非常困难的。市场总是变化的。我们很容易陷入的陷阱是设计一个在回溯测试中表现良好的系统，但其令人鼓舞的表现仅仅是因为在我们回溯测试历史中有利的一段时间里表现突出。因此，我们必须研究我们系统的权益曲线。如果它只在一小部分时间里表现出色，而在其他时间表现平平，我们就应该仔细考虑这种情况。当然，如果前一段时间表现出色，而最近的表现有所恶化，那就更是如此了！
关键的一点是，当我们在某种市场条件下开发一个交易系统时，只有在这种市场条件下，我们才能期待持续的良好表现。因此，*我们希望在我们的开发和测试期间，市场条件变化足够频繁，以便所有可能的市场条件都能得到体现。*而且即使所有条件都表现出来了，缓慢的徘徊也可能引起周期性的延长的不利表现。长时间的出色表现，接着是长时间的糟糕表现，会令人沮丧。
STATN 计划
对于我们这些渴望硬数字的人来说，有一个很好的测试，它比基于目测情节的武断决定更可靠。我已经在 STATN.CPP 程序中提供了这个算法的示例。这个版本读取市场历史文件，并检查市场随时间的趋势和波动性。您可以通过添加其他市场属性(如 ADX 或您使用的任何自定义指标)来轻松修改它。
这个程序的原理很简单，但却令人惊讶地揭示了市场异常。它基于这样一种想法，即在特定市场条件下开发的交易系统(如上升或下降趋势，高或低波动性)在其他市场条件下可能会失去盈利能力。在大多数情况下，我们希望看到反映在我们的指标中的这些条件在规则和合理随机的基础上变化，以便我们开发的系统将尽可能多地经历投入使用时会遇到的各种条件。缓慢的徘徊是危险的非平稳性的本质；市场属性可能会在一段时间内保持在一个状态，然后在另一段时间内变为另一个状态，同样会影响我们的指标。这使得开发健壮的模型变得困难。粗略地说，平稳性等于行为的一致性。
使用以下命令调用该程序:
STATN Lookback Fractile Version Filename
让我们来分解这个命令:
-
Lookback:历史棒线的数量，包括当前棒线，用于计算市场的趋势和波动性。 -
Fractile:趋势和波动性的分位数(0–1)，用作差距分析的上/下阈值。 -
Version: 0 表示原始指标，1 表示差异原始指标，> 1 表示指定原始减去扩展原始。详见第 14 页。 -
Filename:格式为YYYYMMDD Open High Low Close的市场历史文件。
使用真实市场数据的例子将出现在第 17 页。首先，我们探索一些代码片段。参见 STATN。完整上下文的 CPP。
该程序遍历市场历史，计算趋势(最小二乘线的斜率)和波动性(平均真实范围)。它查找对应于指定分位数的分位数；0.5 将是中间值。对于每根棒线，它决定趋势和波动的当前值(或它们的修改值，如后所述)是小于分位数还是大于或等于分位数。每次状态改变时(从上面到下面或从下面到上面),它会记录已经过了多少个小节并记录下来。例如，如果下一个条形的状态发生变化，则计数为 1。如果状态在下一个条之后改变了一个条，则计数为 2，依此类推。对于 1、2、4、8、16、32、64、128、256、512 和大于 512 的条计数，定义了 11 个仓。当程序结束时，它打印仓位计数，一个趋势表和一个波动表。
Version参数需要更多一点的解释，其理由将推迟到下一节。现在，请理解，如果用户将其指定为 0，趋势和波动指标将完全按照计算结果使用。如果它是 1，每个指示器的当前值通过减去它在lookback条之前的值来调整，使它成为一个经典振荡器。如果它大于 1，则通过使用Version * Lookback的回看减去该值来调整当前值，使其成为另一种振荡器。后两个版本要求实际回看大于用户指定的回看，如以下代码所示:
if (version == 0)
full_lookback = lookback ;
else if (version == 1)
full_lookback = 2 * lookback ;
else if (version > 1)
full_lookback = version * lookback ;
nind = nprices - full_lookback + 1 ; // This many indicators
如果nprices是价格条的数量，我们丢失其中的full_lookback–1，得到指标的nind值，如前面代码的最后一行所示。
以下代码块显示了趋势指标(可能已修改)的计算。波动性也是如此。对于每一遍，k是指示器当前值的索引。我们必须从足够远的指标历史开始，以涵盖完整的回顾。
for (i=0 ; i<nind ; i++) {
k = full_lookback - 1 + i ;
if (version == 0)
trend[i] = find_slope ( lookback , close + k ) ;
else if (version == 1)
trend[i] = find_slope ( lookback , close + k ) –
find_slope ( lookback , close + k - lookback ) ;
else
trend[i] = find_slope ( lookback , close + k ) –
find_slope ( full_lookback , close + k ) ;
trend_sorted[i] = trend[i] ;
}
对值进行排序以找到用户指定的分位数，然后计算每个箱中的计数。
qsortd ( 0 , nind-1 , trend_sorted ) ;
k = (int) (fractile * (nind+1)) - 1 ;
if (k < 0)
k = 0 ;
trend_quantile = trend_sorted[k] ;
gap_analyze ( nind , trend , trend_quantile , ngaps , gap_size , gap_count ) ;
在调用gap_analyze()之前，我们必须做一些准备，为其提供间隙尺寸的边界。如果你想的话，可以随意修改。分析代码出现在下一页。
#define NGAPS 11 /* Number of gaps in analysis */
ngaps = NGAPS ;
k = 1 ;
for (i=0 ; i<ngaps-1 ; i++) {
gap_size[i] = k ;
k *= 2 ;
}
该例程仅保留一个标志above_below，如果当前值等于或高于阈值，则该标志为真 (1)，如果低于阈值，则为假 (0)。对于每次循环，如果指示器仍然在阈值的同一侧，则计数器递增。如果它转换了方向，则相应的 bin 递增，并且计数器复位。到达数组末尾相当于翻转两边，所以最后一个系列才算。
void gap_analyze (
int n ,
double *x ,
double thresh ,
int ngaps ,
int *gap_size ,
int *gap_count
)
{
int i, j, above_below, new_above_below, count ;
for (i=0 ; i<ngaps ; i++)
gap_count[i] = 0 ;
count = 1 ;
above_below = (x[0] >= thresh) ? 1 : 0 ;
for (i=1 ; i<=n ; i++) {
if (i == n) // Passing end of array counts as a change
new_above_below = 1 - above_below ;
else
new_above_below = (x[i] >= thresh) ? 1 : 0 ;
if (new_above_below == above_below)
++count ;
else {
for (j=0 ; j<ngaps-1 ; j++) {
if (count <= gap_size[j])
break ;
}
++gap_count[j] ;
count = 1 ;
above_below = new_above_below
;
}
}
}
通过振荡改善位置平稳性
提高指标稳定性的一个简单但通常有效的方法，至少就其中心趋势而言，是计算其相对于某个相关“基础”值的值*。最常见也是最有效的方法是减去滞后值，滞后值通常(但不一定)是指标的回望值。例如，我们可以计算最近 20 个价格的趋势，并从中减去 20 根棒线之前该指标的值。*
一个相似但不完全相同的方法是计算当前时间的指标，但有两个不同的回顾，一个短，一个长。从短期指标中减去长期指标，得到一个更稳定的修正指标。
这两种方法都涉及到重要的权衡。可能是指标的实际值承载了重要的信息。刚刚描述的两个修改放弃了实际值，而选择了相对值。根据我的经验，后一个值通常比实际值携带更多的预测信息，并且在几乎所有情况下都有更好的稳定性。但这并不具有普遍性，这种权衡必须牢记在心。
如果这种权衡是一个问题，请记住，第一种方法，即找出指标的当前值和滞后值之间的差异，是最“强大”的，因为它通常导致最大的稳定性，同时也丢弃了关于真实当前值的大多数信息。第二种方法更像是一种妥协。此外，通过调整长期回顾，人们可以对这种权衡施加很大的控制。增加长期回顾可以更好地保存当前值的信息，但代价是平稳性的改善较少。
在下一页，我们可以看到 STATN 程序制作的两个表格，回顾值为 100，S&P 100 指数 OEX 的分位数为 0.5(中位数)。上表是趋势，下表是波动。第一列是原始指标；第二个是Version =1，滞后差；而第三个是Version =3，给出 300 棒的长期回看。
Trend with Lookback=100, Fractile=0.5
Gap Version=0 Version=1 Version=3
1 3 1 0
2 3 1 0
4 2 2 2
8 5 2 1
16 4 3 4
32 14 2 12
64 22 14 25
128 29 54 33
256 18 15 21
512 3 1 1
>512 0 0 0
Volatility with Lookback=100, Fractile=0.5
Gap Version=0 Version=1 Version=3
1 13 41 19
2 6 13 6
4 2 9 13
8 2 8 6
16 4 9 4
32 2 10 10
64 3 12 8
128 5 25 10
256 9 23 18
512 2 5 9
>512 6 0 1
在这个趋势表中，我们看到原始指标有三个长时间段，在这三个时间段中，该指标保持在其中值的同一侧。这些周期大于 256 根连续棒线，也许长达 512 根棒线，超过两年！两个修改版本只有一个这样的周期。
波动性的情况甚至更严重，原始指标有六个时间段大于 512 棒线，波动性位于中值的同一侧。修改极大地改善了这种情况，尽管在下一个较低的级别会有显著的恶化。波动一般具有极端的非平稳性。
极端平稳归纳
刚刚描述的两种方法只在指标的中心趋势中引入平稳性。这很重要，可以说是平稳性最重要的品质。如果一个指标缓慢波动，长时间停留在高值，然后长时间移动到低值，这个指标在交易系统中的效用可能会受损。基于此类指标的系统很容易长期亏损，甚至停止交易。当然，在某些情况下，停止交易是有用的；如果你有几个互补的系统，如果每个系统都在盈利交易和不交易之间交替，那就太好了。不幸的是，在现实生活中，这样的系统只是例外，而不是规则。
但是一个指标有无数种不稳定的方式。集中趋势(平均值)通常是最重要的，其次是方差。如果一个指标在很长一段时间内变化很小，然后在随后的很长一段时间内变化很大，则该指标将受损。
有一种简单的方法可以使均值、方差或两者达到一个极端但可控的程度。只需回顾指标最近值的移动窗口，并计算该窗口内的平均值(如果指标表现良好)或中值(如果偶尔出现极值)。从当前值中减去该值，以在中心趋势中引入稳定性。如果窗口很短，效果会很明显，足以克服几乎任何程度的不稳定性。同样，您可以计算移动窗口的标准差(如果指标表现良好)或四分位间距(如果出现野值)。将(可能居中的)当前值除以该量，以引起方差的平稳性。
没有提供这种方法的例子，因为它是一种简单的计算。请记住，长窗口将保留有关指标实际值的大量信息，同时提供很少的非平稳性减少。相反，短窗口会破坏几乎所有关于实际值的信息，使一切都与最近的历史相关，从而导致巨大的平稳性。
用熵度量指标信息
几十年前，贝尔实验室的 Claude Shannon 开发了一种严格且极其强大的方法来量化一条消息可以传达的信息量。这与交易系统的开发有关，因为从最近的市场历史计算的指标可以被认为是来自市场的信息，传达了关于市场的当前和可能的未来状态的信息。如果我们可以量化指标中的平均信息，我们就可以了解该指标的潜在价值。更好的是，我们可以修改指标，增加其信息含量。并非巧合的是，这些增加信息的修改与众所周知的提高预测模型性能的修改完全相同。这是一个值得研究的领域。
我们将采取一种肤浅的、直观的方法来量化指标中的平均信息。关于这个主题的更详细的探索，请参阅我的两本书c++中的数据挖掘算法或者评估和改进预测和分类。
假设需要传达一条信息，这条信息就是一道选择题的答案。也许这是一个简单的二元选择，例如“市场处于上升趋势状态”与“市场处于下降趋势状态”或许再详细一点，比如四种可能的情况:“市场处于强涨/弱涨/弱跌/强跌”状态。现在添加一个限制，即消息必须是二进制的，即一个或多个 1 和 0 的字符串。显然，第一个问题的答案可以用一个二进制位给出，而第二个问题的答案将需要两个位来涵盖四种可能的市场状态(00，01，10，11)。一般来说，如果有 K 个可能的答案，那么我们将需要消息中的 log 2 ( K )位来传达正确的答案。
量化信息价值的一个好方法是它所传达的信息的位数。为消息赋值的一种不太清楚但更有用的方法是通过接收消息来消除不确定性的位数。假设你参加了一个总共有 1024 张彩票的抽奖，其中一张是你的。获胜者的身份可以用 log 2 (1024)=10 位编码。在你收到任何信息之前，你对获胜者的身份有 10 点不确定性。相当于，每个条目都有 1/1024 的机会成为赢家。
收到一条信息，回答一个简单的问题:你是否赢得了彩票。让我们计算这两个可能答案的价值。如果答案是你中了彩票，则解决了概率为 1/1024 的事件，给该特定消息的值为 log2(1024)=–log2(1/1024)= 10 位。如果答案是您没有赢，则概率为 1023/1024 的事件已经解决，给该特定消息的值为–log2(1023/1024)= 0.0014 位。
大多数人(和计算机)不使用 base 2 中的日志。相反，他们使用自然对数。当这样做时，信息的单位是 nat，而不是比特。所以，在讨论的例子中，你赢了的一个答案的值是–log(1/1024)= 6.93 NATs，失望答案的值是–log(1023/1024)= 0.00098 NATs。
我们只是计算了每个答案的价值。但是我们也对消息的预期价值感兴趣。回想一下，离散随机变量的期望值是每个值乘以该值的概率的乘积之和。因此，消息的期望值是你赢得答案的概率乘以它的值，加上你没有赢得答案的概率乘以它的值*。这是 1/1024 –log(1/1024)+1023/1024 –log(1023/1024)= 0.0077 纳特。这个期望值被称为消息的熵，并且被符号化为 H 。*
我们可以更严谨一些。设 χ 是一个集合，它枚举了一个消息流 X 中每一个可能的答案。因此， χ 可能是{ 大幅度上升趋势、小幅度上升趋势、小幅度下降趋势、大幅度下降趋势 }。当我们观察到一个 X 的值时，我们称之为 x ，根据定义它总是 χ 的成员。这就写成了 x ∈ χ 。设 p ( x )为 x 被观测到的概率。那么 X 的熵由方程 2-1 给出。在这个等式中，0*log(0)被定义为零。
)
(2-1)
我们在没有证明的情况下陈述，当每个可能的答案( x 的值)具有相等的概率时，消息流 X 的熵(平均信息量)最大，并且这个最大熵是 log( K ，其中 K 是 x 的可能值的数量。因此，我们将对H(X)/log(K)的值最感兴趣，因为这个数字的范围从零(消息流不传递任何信息)到一(消息流传递最大可能的信息量)。这个比值称为相对熵或比例熵。
最后，我们可以将所有这些(高度简化的)理论与自动市场交易联系起来。我们要做的是筛选我们交易系统中使用的所有指标的相对熵。如果相对熵很小，我们应该考虑用不同的方法计算指标，也许可以采用一种简单的方法，比如应用非线性变换来增加相对熵。在我自己的工作中，我希望相对熵至少为 0.5，甚至更高，尽管这个阈值是非常随意的。
有几个注意事项要记住。首先，理解熵是信息内容的度量，但是我们不知道这些信息是否与手头的任务相关。一个指标可能在预测市场波动是否会在接下来的一周爆发方面做得非常出色。但是，如果我们的目标是确定我们是做多还是做空，这个信息丰富的指标对我们的项目来说可能毫无价值。然而，熵可以被认为是信息内容的上限，所以如果熵很小，我们的指标可能没有什么价值。
第二，可能发生的情况是，无论我们做什么修改来增加指标的熵，实际上都会妨碍指标的性能。也许我们最初的想法作为一个指标做得很好，但是当我们应用一个看似无害的改变来大大增加它的熵时，它在我们交易系统中的效用下降了。这是可能发生的。但是请理解，这两种情况，尤其是第二种情况，是不寻常的例外。在绝大多数情况下，增加指标的熵可以显著提高其性能。
计算指标的相对熵
从指标的历史值计算指标的相对熵的最简单也可能是最好的方法是，将指标的整个范围划分为多个区间，这些区间以相等的间距划分该范围，计算落入每个区间的事例的比例，并使用等式 2-1 来计算熵。将这个量除以箱数的对数，得到相对熵。请注意，将范围划分为包含相同数量事例的箱是没有意义的，因为这将始终给出一个相对熵。相反，仓必须由总范围的相等数字部分来定义。下面是一个简单的子例程:
double entropy (
int n , // Number of data values
double *x , // They are here
int nbins , // Number of bins, at least 2
int *count // Work area nbins long
)
{
int i, k ;
double minval, maxval, factor, p, sum ;
minval = maxval = x[0] ;
for (i=1 ; i<n ; i++) {
if (x[i] < minval)
minval = x[i] ;
if (x[i] > maxval)
maxval = x[i] ;
}
factor = (nbins - 1.e-10) / (maxval - minval + 1.e-60) ;
for (i=0 ; i<nbins ; i++)
count[i] = 0 ;
for (i=0 ; i<n ; i++) { // Count the number of cases in each bin
k = (int) (factor * (x[i] - minval)) ;
++count[k] ;
}
sum = 0.0 ;
for (i=0 ; i<nbins ; i++) { // Sum Equation 2-1
if (count[i]) {
p = (double) count[i] / n ;
sum += p * log ( p ) ;
}
}
return -sum / log ( (double) nbins ) ;
}
在前面的代码中，我们必须对将数据值映射到 bin 的因子的计算进行两次微小的旋转。分子被略微减小，以确保在最后一个 bin 之后没有映射到不存在的“bin”。修改分母，以确保在所有数据值相等的病态情况下，我们不会被零除。最后一个循环只是对等式 2-1 求和，我们通过将熵除以其最大可能值来得出相对熵。
熵影响预测模型的质量
熵作为衡量指标信息含量的有用性不仅仅是理论上的空谈。不管是不是巧合，熵与我们训练有效预测模型的能力高度相关。这是因为高熵与指标范围内数据值的大致相等分布相关，并且在大多数情况下，当模型的指标具有这样的分布时，模型训练最有效。
最常见的有问题的低熵情况是当有一个或多个极端异常值时。许多模型训练算法会将一个异常值视为说了一些重要的事情，并对该异常值投入大量注意力。这减少了对大量“正常”病例的关注。图 2-1 展示了这种情况的一个有点简单但通常很现实的例子。这是一个处理两个类的线性分类器，应该是一个玩具般简单的问题。虚线示出了实现完美分类的线性边界。但是左下角的那个例子，是 X2 指标的异常值，把边界线往它的方向拖，严重损害了分类质量。尽管这个特定的例子以线性分类器为特征，但是即使是具有弯曲边界的非线性分类器，也会经常遭受同样的退化。
图 2-1
离群值会降低性能
由于低熵，不需要离群值来产生性能下降。假设我们有这样一种情况，在这种情况下，与预测/分类完全无关的一些外部条件，在大约一半的情况下，在另外一个优秀的预测器中引起一个大的固定偏移。也许这个变量在一半的情况下有 1.0 左右的值，在这些情况下有很好的表现。同样假设在另一半的情况下，它的值在 100.0 左右，并且它在这个组中也有极好的功效。没有多少模型能够处理这种极低熵的情况。他们会将聚类在 1.0 左右的病例和聚类在 100.0 左右的病例之间的分离视为主导因素，并专注于使用该聚类成员来尝试预测/分类。结果不会很好。
提高指标的熵
如果你测试一个指标，发现它有危险的小熵(小于 0.5 是可疑的；低于 0.1 是严重的，应该进行调查并可能解决)，那么你的第一步应该是重新考虑你的指标想法是如何实现的。也许对你的计算算法做一个简单的修改就能解决这个问题，而不会影响你的想法。以下是一些需要考虑的其他想法:
-
如果你的指标计算除以一个可能变小的值，你就如履薄冰了。
-
你的修改应该与你最初的想法单调相关。换句话说，如果修订前的情况 A 小于修订前的情况 B，那么在修订后应该保持相同的顺序。在其他期望的属性中，这确保了如果某个阈值在修订前的基础上分离病例，则存在将在修订后执行相同分离的阈值。这是保存信息的一个重要品质。
-
截断(将极值重新映射到单个极限值)是解决异常值问题的一种糟糕方法。除此之外，它违反了刚刚列出的优先原则！
-
如果您只有一些罕见的异常值，单调仅尾修改是一个很好的解决方案，它极大地提高了熵，但对指标值的影响相对较小。选择一个中等的百分位数，低异常值可能在 1-10%之间，高异常值可能在 90-99%之间。在这一阈值的“好”侧的情况保持不变。此阈值“异常值”一侧的情况会受到极端单调压缩，如对数压缩。这将在第 29 页详细讨论。
-
如果只有右尾较重或为正偏斜(仅在异常大的情况下)，平方根或立方根变换将处理中等偏斜或异常值，而对数变换应处理严重情况。
-
如果两个尾部都很重，考虑立方根变换。
-
如果两个尾部都非常重或有严重的异常值，双曲正切函数(方程 2-2 和图 2-2 )或逻辑函数(方程 2-3 和图 2-3 )可以提供很好的结果，前提是在应用该函数之前对指标值进行了适当的预标定。如果使用逻辑函数，最好在变换后减去 0.5，使其位于零的中心，这是许多训练算法所欣赏的。
)
(2-2)
)
(2-3)
图 2-3
逻辑函数
图 2-2
TANH 函数
-
如果您的指标应该具有类似于常见统计分布的分布有理论上的原因，那么通过应用该分布的累积分布函数进行转换可能是有效的。例如，许多指标(即两个移动平均线之差的振荡器)有一个漂亮的钟形曲线形状，除了适度沉重的尾部，这几乎是正常的，尾部不严重，但坏到足以引起麻烦。在 STATS 中应用普通 cdf (
normal_cdf())。CPP)会做得很出色。其他指标可能是两个类方差量的比值，在这种情况下，在 STATS 中为 F CDF (F_CDF)。CPP)是理想的。 -
有时，您的指标的分布可能会以一种不直接的方式出现问题。例如，考虑由第 25 页描述的外部条件引起的聚集，其中指标具有良好的紧凑分布，完全没有异常值，但是数据聚集成几个小的聚集。或者它可能有这个问题加上一个重尾，或者两个重尾。当这种情况发生时，有一种粗暴的方法，这种方法很笨拙，但是非常有效和通用，特别是如果您有大量具有代表性的指标值样本。按升序对您的样本进行排序，并选择性地保存以备将来使用。然后，要转换一个值，请使用二分搜索法来限制排序数组中的值。转换后的值是小于或等于转换前值的已排序元素的数量。这产生了具有非常接近完美的相对熵的变换指标。当样本很大，具有充分的代表性，并且很少或没有联系时，这种方法效果最好。作为最后一步，将这个数除以元素总数，然后减去 0.5。这给出了一个范围从-0.5 到 0.5 的值，这个范围对许多训练程序特别友好。
-
刚刚介绍的许多技术力求产生一种尽可能在其范围内均匀的指示剂分布。但有时这并不理想，尽管它有最大熵。当指标的极端值确实具有特殊意义，但这种极端值阻碍甚至阻止预测模型的正确训练时，就会发生这种情况。在这种情况下，你想做的只是驯服尾巴，而不是消灭它们。如果您已经采用了一种变换来生成几乎均匀的分布，但是您希望原始极值映射到足够突出但又不会极端到有问题的值，有一个简单的解决方法:变换到正态分布。这种分布具有钟形曲线形状，其中大多数情况集中在内部，但在两个尾部都有适度的极值。为此，首先应用将指示器映射到几乎均匀分布的任何变换。然后使用逆正态累积分布函数进行第二次变换。这可以通过调用 STATS.CPP 中的函数
inverse_normal_cdf()来完成。得到的指标仍会有极值，但不足以降低模型训练。
单调仅尾部清洗
有时候，你对你的指标的分布总体上是满意的，除了它偶尔会有一个野值，需要驯服以产生适当的熵。或者，预测模型交易系统中的目标偶尔会出现极端情况，妨碍你的训练算法，但你不想过多干预目标值，以免过分扭曲业绩数据。这种情况要求只影响最极端值的转换，而不影响大多数情况。这里有一个很好的处理方法。
这种转变分两步进行。第一步识别标记低尾和高尾的数据值，第二步修改尾(仅)。这里显示的版本以相同的方式处理两个尾部。读者可以很容易地修改它，只处理一个尾部，或者不同地处理上下尾部。
如下所示调用子例程。用户提供了一个工作向量，因为我们必须对原始数据进行排序来定位尾部。调用者指定将被单调压缩的每个尾部的分数(通常很小，可能是 0.01 到 0.1)，保留顺序关系，同时强烈引入离群值。变量cover是保持不变的病例比例。我们将原始数据复制到工作区，并对其进行分类。
void clean_tails (
int n , // Number of cases
double *raw , // They are here
double *work , // Work area n long
double tail_frac // Fraction of each tail to be cleaned (0-0.5)
)
{
int i, istart, istop, best_start, best_stop ;
double cover, range, best, limit, scale, minval, maxval ;
cover = 1.0 - 2.0 * tail_frac ; // Internal fraction preserved
for (i=0 ; i<n ; i++)
work[i] = raw[i] ;
qsortd ( 0 , n-1 , work ) ; // Sort data ascending
识别尾部的一个好方法是检查排序数组中具有指定覆盖范围(包含所需数量的“内部”事例)的每个可能的连续事例集。找出范围(最大值减去最小值)最小的集合。那么，将位于这个最小范围内部集合之外的那些情况标为尾部是合理的。该内部设置将被识别为从istart到istop运行。
istart = 0 ; // Start search at the beginning
istop = (int) (cover * (n+1)) - 1 ; // This gives desired coverage
if (istop >= n) // Happens if careless user has tail=0
istop = n - 1 ;
我们从最左边到最右边，在每一个可能的位置运行这个内部集合。对于每组试验终点，找出范围并跟踪哪个位置的范围最窄。
best = 1.e60 ; // Will be minimum span
best_start = best_stop = 0 ; // Not needed; shuts up LINT
while (istop < n) { // Test every possible position
range = work[istop] - work[istart] ; // This is what we minimize
if (range < best) {
best = range ;
best_start = istart ;
best_stop = istop ;
}
++istart ;
++istop ;
}
此时，我们已经找到了最窄的内部集合。获取其下限值和上限值，并防止粗心的呼叫者。
minval = work[best_start] ; // Value at start of interior interval
maxval = work[best_stop] ; // And end
if (maxval <= minval) { // Rare pathological situation
maxval *= 1.0 + 1.e-10 ;
minval *= 1.0 - 1.e-10 ;
}
最后一步是修改尾部(只修改尾部)。通过使用maxval–minval作为缩放常数，我们保持该过程不受数据缩放变化的影响。limit变量控制变换后的尾部值位于内部范围之外的程度。使用(1.0–cover)这个因素是我自己的启发，在我看来是合理的。不同意也可以，想改就改。
读者应该检查这段代码并确认limit确实定义了转换值的偏离极限，即最小值和最大值(以及它们之间的值！)保持不变，并且变换是单调的(保序的)。
limit = (maxval - minval) * (1.0 - cover) ;
scale = -1.0 / (maxval - minval) ;
for (i=0 ; i<n ; i++) {
if (raw[i] < minval) // Left tail
raw[i] = minval - limit * (1.0 - exp ( scale * (minval - raw[i]) ) ) ;
else if (raw[i] > maxval) // Right tail
raw[i] = maxval + limit * (1.0 - exp ( scale * (raw[i] - maxval) ) ) ;
}
}
熵程序
文件熵。CPP 包含一个完整的程序，演示了从市场价格历史文件计算的各种指标的熵的计算。其中两个指标是第 13 页描述的 STATN 程序中的趋势和波动指标。此外，这里的版本参数与 STATN 程序中的相同，尽管它在熵的上下文中不如在平稳性中有趣。
使用以下命令调用该程序:
ENTROPY Lookback Nbins Version Filename
让我们来分解这个命令:
-
Lookback:历史棒线的数量，包括当前棒线，用于从市场价格历史计算指标。 -
Nbins:用于计算熵的箱数。对于一个至少有几千条记录的市场历史来说，大约 20 个左右的仓位是好的，尽管实际上这个数字并不太重要。如果少量改变容器的数量会导致计算出的熵发生很大的变化，那么数据或读者设计的任何自定义指标都有问题。绘制直方图！ -
Version: 0 表示原始指标，1 表示差异原始指标，> 1 表示指定原始减去扩展原始。详见第 14 页。 -
Filename:格式为YYYYMMDD Open High Low Close的市场历史文件。
计算以下指标，并打印其最小值、最大值、中值和相对熵。
-
趋势是由最小二乘拟合定义的每根棒线的(对数)价格变化。
-
波动率是根据标准定义计算的平均真实范围。
-
扩展是一个故意设计得很差的指标，它展示了而不是如何定义一个指标，以及低熵如何揭示问题。收盘价的范围(最大收盘价减去最小收盘价)是为覆盖指定回看距离一半的最近价格计算的。然后计算相同的量，但是滞后于回看的一半。扩展指标是最近的范围除以旧的范围，分母略微增加以防止被零除。该指标揭示了波动性(价格范围)的粗略度量是增加、减少还是保持不变。
-
RawJump 衡量最近收盘价与最近指数平滑收盘价的对比情况。这个量揭示了市场是突然上涨还是下跌，还是保持不变。它的两端偶尔会有异常值，因此熵值很低。
-
CleanedJump 是第 29 页描述的单调尾部平滑应用于每个尾部的外部 5%后的 RawJump 。
当熵程序在标准普尔 500 市场历史上运行时，使用 20 根棒线和 20 个仓来计算熵，计算第一列中的相对熵值。当回看下降到七个小节时，我们得到的结果显示在第二列中。
| *趋势* | Zero point five eight | Zero point four eight three | | *波动性* | Zero point six three nine | Zero point five five nine | | *膨胀* | Zero point four six one | Zero | | *拉弗 Jump* | Zero point four eight four | Zero point three nine five | | 清洁跳跃 | Zero point nine five eight | Zero point nine five two |特别是对于较短的回望，趋势和波动指标的相对熵勉强可以接受。他们可以稍微调整一下。一些温和的东西可能会做。扩张指标，通过使用不稳定的比率故意设计得很差，在回顾七根棒线时变得一文不值。请特别注意这样一个事实，即 RawJump 指标的相对熵从差到优，只不过是清理了外部 5%的尾部，其他什么也没碰。
三、优化问题
正则化线性模型
如果我只能让这本书的读者有一个想法，那就是:你的指标的强度比用它们来指示交易的预测模型的强度重要得多。这些年来我见过的一些最好的、最稳定的、最赚钱的交易系统，都是使用简单的线性或接近线性的模型，以高质量的指标作为输入。我也见过太多的人向一些现代的、高度复杂的非线性模型提供边际指标，徒劳地希望该模型将奇迹般地克服垃圾进、垃圾出规则。它不会发生。当我在开发一个新的交易系统时，我首先转向线性模型，只有当我看到明显的优势时，才转向非线性模型。
与复杂的非线性模型相比，在基于预测模型的交易系统中使用线性模型有很多优点。
-
线性模型不太可能过度拟合训练数据。因此，训练偏差被最小化。这个问题在第 121 页开始的部分有更详细的论述。
-
线性模型比许多或大多数非线性模型更容易理解。理解指标值与交易决策的关系是预测模型的一个非常有价值的属性。
-
线性模型的训练速度通常比非线性模型快。在第七章中，我们将探索需要频繁重新训练的强大测试算法，因此快速训练是一个主要优势。
-
很容易将线性模型转换成非线性模型，而不会严重破坏刚刚列出的特性。这将在第 67 页讨论。
-
随着算法复杂性的适度增加，很容易惩罚过度复杂/强大的线性模型，这是一个重要但经常被忽视的正确训练的一部分，称为正则化。
正则化模型概述
正如从第 121 页开始的部分将详细讨论的那样，当模型错误地将随机噪声与真实的、可重复的模式混为一谈时，预测模型的设计和训练中就会出现一个常见而严重的问题。这叫做过拟合。因为根据定义，噪声不会重复，一个过度拟合的模型在投入使用时会表现不佳。
普通的线性模型比大多数非线性模型更不容易过度拟合，尤其是那些极其复杂和强大的模型。但是，即使普通的线性模型也可能过度拟合训练数据，这通常是因为使用了太多的预测器。至少有两种常见且适度有效的方法来处理由于预测器数量过多而导致的过拟合问题。
-
减少预测器的数量:最常用的方法是向前逐步选择。选择一个最有效的预测值。然后，在第一个预测因子存在的情况下，选择一个增加最大预测能力的预测因子。然后选择第三个，依此类推。这种方法的一个严重问题是很难找到第一个单独的预测器。在大多数应用中，是几个预测器的相互作用提供了动力；没有一个单一的预测器做得很好。事实上，可能是 A 和 B 一起做了一个神话般的工作，而其中任何一个都是毫无价值的。但如果 C 恰好表现一般， C 可能会先被选中，结果是 A 和 B 都没有进入比赛，这个优秀的组合也就失去了。还有其他更高级的变异，如逆向选择或子集保留。这些在我的书c++ 中的数据挖掘算法中有详细讨论。然而，每种方法都有自己的问题。
-
将线性回归方程中的系数向零收缩，远离它们的“最佳”最小二乘值。这可能非常有用，因为它让我们保留所有相关的预测器和它们的联合关系信息，同时降低它们学习随机噪声以及更显著的真实模式的能力。但这不是一项微不足道的任务。
有效的模型设计和训练(线性或非线性)的目标是执行这两个修正中的一个或(通常)两个。有许多方法可以做到这一点，其中大多数都涉及到对模型的复杂性施加惩罚，这一过程被称为正则化。区分各种方法的是复杂性的定义和相关惩罚的性质。当应用于线性模型时，此处显示的方法特别强大，因为它可以根据用户的判断进行一种或两种修复，并且它以一种易于理解和快速训练的方式进行。此外，有一个简单的交叉验证方案，让我们优化复杂性降低超参数。它真的很漂亮。
首先我们必须设计一些符号。不失一般性，所有后续开发将假设所有预测值 x 已被标准化为具有零均值和单位方差。这极大地简化了相关的方程以及相关的计算机代码。如果需要，简单的代数运算可以恢复原始变量的系数。
-
N -案例的数量。
-
K -预测变量的数量。
-
xij——案例 i 的预测值 j 。
-
xI-预测向量( K long)为例 i 。这是一个列向量，是表示所有 j 的集合xij的一种方便的表示法。
-
yI——事例 i 的目标变量。
-
β -方程 3-1 表示的线性模型中的 K 系数。这是一个列向量。
-
β0——方程 3-1 表示的线性模型中的标量常数。
-
α -控制正则化的类型的范围从 0 到 1 的常数。
-
λ -控制正则化的度的非负常数。
基本的线性模型认为目标变量的期望值等于预测值的加权组合，加上一个常数。这在方程 3-1 中以矢量形式显示。
)
(3-1)
我们已经假设预测值已经标准化为零均值和单位方差。在这种情况下，β 0 等于目标的平均值。(这一简单结果的推导在许多标准统计学文本中都有说明。)如果我们假设目标变量也已经标准化，我们可以在开发和编程中获得更多的简单性。我们由此知道β 0 = 0，因此在所有后续工作中可以忽略它。在这种情况下，我们预测目标的标准化值。要获得原始目标的预测值，只需进行非标准化:乘以标准偏差，然后加上平均值。即使有了这个额外的假设，所有原始值的系数和β 0 也很容易用基本代数得到。
找到β权重的最佳值的传统方法是计算最小化均方误差的那些值，即每个预测值yˇI和真实目标值 y i 之间的均方差。但是在我们的正则化版本中，我们给误差增加了一个惩罚项，惩罚是β权重集的函数。这显示在方程式 3-2 中。在这个等式中，乘数 2 可以被吸收到λ或 P α ()中，但以这种方式显示是为了澄清一些中间推导，我们在这里不做介绍，因为我们的重点将放在对模型编程至关重要的等式上。要了解全部细节，请参阅 Friedman、Hastie 和 Tibshirani 的优秀论文“通过坐标下降的广义线性模型的正则化路径”(统计软件杂志，2010 年 1 月)。注意，λ控制惩罚项的影响，如果λ = 0，我们有普通的最小二乘解。将 α 的下标应用于罚函数 P 来阐明它控制罚函数性质的事实。还要注意，方程 3-2 是那篇论文中误差的两倍，没有实际后果。
)
(3-2)
罚函数 P α 是权重向量的二范数(平方和)和一范数(绝对值和)的加权和，相对权重由 α 参数确定。这显示在方程式 3-3 中。
)
(3-3)
α 的值对罚函数的性质有深远的影响，该值的范围可以从 0 到 1。这两个极端值有一个共同的名字，许多开发人员都很熟悉。当 α = 0 时，我们有岭回归，当 α = 1 时，我们有套索。这两个极端模型之间的差异最好通过考虑当存在高度相关的预测因子时会发生什么来说明。岭回归倾向于给相关集合中的所有预测因子分配大致相等的权重，从所有预测因子中提取大致相等的贡献。事实上，如果有一组 m 完全相关(归一化后相同)的预测值，岭回归将为每个预测值分配一个β权重，该权重等于 1/ m 乘以在没有其他预测值的情况下分配给其中一个预测值的权重。
套索( α = 1)对一组高度相关的变量做出相反的反应。它倾向于挑选对模型最有用的一个，给它分配一个相对较大的权重，给相关集合的其他成员分配零权重，实际上是将它们从模型中移除，而不是给所有成员分配较小的、相似的权重。
α = 1 的一个潜在问题是，如果碰巧有两个或更多的预测器完全相关，套索就会失去判断哪个是最好的想法，因为它们都同样有用。训练算法在数值上变得不稳定。由于这个原因，除非你确信数据中不存在这种退化，否则如果你想使用套索模型，你应该将 α 设置为非常接近 1 但不完全是 1 的值。这个模型将几乎等同于一个真正的套索，但不会因为完美或接近完美的相关性而不稳定。
在大多数金融交易开发中，有大量的预测器以“墙上的意大利面条”的方式扔向模型，通常最好将 α 设置为 0 到 1 之间的一个值，以达到最佳效果。对于任何固定的λ，零系数(从模型中排除的变量)的数量随着 α 从零到一单调增加。当 α = 0 时，所有变量都包括在内，然后对于更大的 α 值，它们倾向于一个接一个地退出(它们的β权重变为零)。以这种方式，开发者可以设置 α 的值来支持期望的稀疏度。
在将这里描述的模型与普通的线性回归进行比较时，读者应该记住三件事:
-
当我们这样惩罚模型时，我们得到的解不再是最小二乘解。计算的β权重将产生超过普通线性回归的均方误差。在大多数实际应用中，这是一件好事，因为它产生了更好的泛化能力。这就是这种方法的全部意义！然而，从表面上看，这似乎是违反直觉的，好像我们在故意削弱这个模型。但这正是我们正在做的，以减少它错误地学习随机噪声的能力。
-
我们应该对这个模型如何处理强相关的预测值感到特别高兴。普通的线性回归通常对这种情况有可怕的反应，将一些系数放大到巨大的正值，然后通过将其他系数放大到巨大的负值来进行补偿，使一个相关变量与另一个相关变量保持微妙的平衡关系。
-
这种正则化模型通常找到候选预测值的一个子集，通常是一个小的子集，就像普通逐步包含的情况一样。但它的方法非常不同，远远优于逐步包含。后者采取有序的全有或全无的方法；一旦一个变量被包含，它将永远存在。但是正则化的线性模型是逐渐运行的，缓慢地收敛于预测因子的理想子集。在其他变量存在的情况下，变量的值会时好时坏。这使得最终子集更有可能是真正最优的。
保证收敛的β调整
有一个简单的公式，根据该公式，给定训练数据，连同超参数λ和 α ，我们可以有效地计算出调整后的β权重，从而降低等式 3-2 中所示的误差标准。在所有实际条件下，该误差准则具有单个局部最小值，该局部最小值也是全局最小值。因此，即使对于大问题，简单的权重旋转也能保证收敛，通常非常快。在这一节中，我们将介绍这一调整公式，省略了在前面引用的论文中可以找到的许多推导细节。我们将很快看到如何使用这个公式来实现一个高效稳定的训练算法。
首先将模型的残差定义为其预测误差，如方程 3-4 所示。
)
(3-4)
为每个预测值 j 定义一个我称之为自变量 j 的项，如等式 3-5 所示。这是计算的缓慢部分，因为它需要对所有情况下的乘积求和。定义软阈值算子 S ()，如公式 3-6 所示。然后，减少误差标准的β j 的新值由等式 3-7 给出。
)
(3-5)
)
(3-6)
)
(3-7)
差异案例加权
在某些应用程序中(虽然在市场交易应用程序中并不常见)，将一些案例评定为比其他案例更重要，从而指导训练算法更加关注减少重要案例的错误，这可能是有用的。上一节中展示的测试版更新公式很容易修改来实现这一功能。
让 N 个案例权重被表示为wI，其中这些权重总和为 1。等式 3-8 给出了软阈值算子的自变量，等式 3-9 给出了更新后的β权重。
)
(3-8)
)
(3-9)
感兴趣的读者最好做一个简单的练习。假设所有权重都等于 1/ N ，那么不存在差分加权。解决方程 3-8 简化为方程 3-5 ，方程 3-9 简化为方程 3-7 的问题。如果没有马上看到(剧透预警！)，记住预测因子已经标准化为单位方差了。因此，对于每个 j ，所有情况下的总和xij的平方等于 N 。
协方差更新的快速计算
如果情况( N )比预测值( K )多得多，这是市场交易中的常见情况，那么有一个替代公式来计算β权重更新，它比“原始”公式 3-5 和 3-8 快得多。基本公式由公式 3-10 给出。
)
(3-10)
如果不使用差分案例加权，则所有 j 、茵纳、 j 、欣纳、JK、 3-12 由公式 3-11 给出 Xss 、j、 = 1。这些表达式的推导在前面引用的论文中给出。
)
(3-11)
)
(3-12)
如果我们使用差分加权，我们需要方程 3-13 、 3-14 和 3-15 。这些推导在引用的论文中没有给出，但是它们很容易从方程 10 开始并遵循非加权步骤得到，记住预测值是标准化的。
)
(3-13)
)
(3-14)
)
(3-15)
注意，等式 3-13 到 3-15 仅取决于训练数据和权重，因此它们可以在训练开始时仅计算一次。而且方程 3-10 ，每次迭代都必须求值，只涉及对 K 项求和，而不是对 N 项求和，当 K < < N 时，节省的时间是巨大的。
预备代码
我们以一些片段开始展示代码，这些片段说明了我们如何准备训练模型的关键部分。封装该模型及其所有训练算法的整个CoordinateDescent类的完整源代码在 CDMODEL.CPP 文件中。
程序员将首先调用构造函数，如下所示。这里我们将跳过它的代码，因为它只涉及内存分配和其他简单的内务处理。暂时忽略nl参数；这个后面会讨论。其他参数不言自明。
CoordinateDescent::CoordinateDescent (
int nv , // Number of predictor variables
int nc , // Number of cases we will be training
int wtd , // Will we be using case weights? 1=Yes, 0=No
int cu , // Use fast covariance updates rather than slow naive method
int nl // Number of lambdas we will be using in training
)
在我们构造了一个CoordinateDescent对象之后，必须调用一个成员函数来输入训练数据，计算一些准备中的东西。
void CoordinateDescent::get_data (
int istart , // Starting index in full database for getting nc cases of training set
int n , // Number of cases in full database (we wrap back to the start if needed)
double *xx , // Full database (n rows, nvars columns)
double *yy , // Predicted variable vector, n long
double *ww // Case weights (n long) or NULL if no weighting
)
在这个调用中，我们可以在数据集中指定一个起始索引(用于预测值、目标值和可选权重)。构造函数调用(nc)中指定的事例数将从索引istart开始从xx、yy和ww(如果使用)中获取。如果在获得nc病例之前到达数据的结尾，它会绕到数据集的开头。我们稍后会看到这种包装是如何有用的。
get_data()例程首先将预测值和目标值保存在私有数组中，然后通过减去平均值并除以标准偏差来标准化它们。这里没有显示这些简单的操作。如果要使用差分加权，权重被缩放为总和为 1(因此用户无需担心这一点)，并且使用等式 3-13 计算XSSvec。与重量相关的代码如下:
if (w != NULL) {
sum = 0.0 ;
for (icase=0 ; icase<ncases ; icase++) {
k = (icase + istart) % n ; // Wrap to start if needed
w[icase] = ww[k] ;
sum += w[icase] ;
}
for (icase=0 ; icase<ncases ; icase++)
w[icase] /= sum ;
for (ivar=0 ; ivar<nvars ; ivar++) {
xptr = x + ivar ;
sum = 0.0 ;
for (icase=0 ; icase<ncases ; icase++) // Equation 3-13
sum += w[icase] * xptr[icase*nvars] * xptr[icase*nvars] ;
XSSvec[ivar] = sum ;
}
}
如果我们使用快速协方差更新方法，这是当案例多于预测值时的明智做法，我们必须按照上一节所述计算Yinner和Xinner。注意Xinner是一个对称矩阵，但是我们还是保存了整个矩阵。这浪费了非常便宜的存储器，但是更简单的寻址节省了非常昂贵的时间。
在下面的代码中，我们一次处理一个变量。在第一种情况下，通过使用指针xptr获得当前变量的偏移量，可以简化寻址。此后，我们只需向下跳一格就可以得到这个变量。
for (ivar=0 ; ivar<nvars ; ivar++) {
xptr = x + ivar ;
sum = 0.0 ; // Do Yinner
if (w != NULL) { // Weighted cases
for (icase=0 ; icase<ncases ; icase++)
sum += w[icase] * xptr[icase*nvars] * y[icase] ; // Equation 3-14
Yinner[ivar] = sum ;
}
else {
for (icase=0 ; icase<ncases ; icase++)
sum += xptr[icase*nvars] * y[icase] ; // Equation 3-11
Yinner[ivar] = sum / ncases ;
}
// Do Xinner
if (w != NULL) { // Weighted
for (jvar=0 ; jvar<nvars ; jvar++) {
if (jvar == ivar)
Xinner[ivar*nvars+jvar] = XSSvec[ivar] ; // Already computed, so use it
else if (jvar < ivar) // Matrix is symmetric, so just copy
Xinner[ivar*nvars+jvar] = Xinner[jvar*nvars+ivar] ;
else {
sum = 0.0 ;
for (icase=0 ; icase<ncases ; icase++)
sum += w[icase] * xptr[icase*nvars] * x[icase*nvars+jvar] ; // Eq (3-15)
Xinner[ivar*nvars+jvar] = sum ;
}
}
} // If w
else { // Unweighted
for (jvar=0 ; jvar<nvars ; jvar++) {
if (jvar == ivar)
Xinner[ivar*nvars+jvar] = 1.0 ; // Recall that X is standardized
else if (jvar < ivar) // Matrix is symmetric, so just copy
Xinner[ivar*nvars+jvar] = Xinner[jvar*nvars+ivar] ;
else {
sum = 0.0 ;
for (icase=0 ; icase<ncases ; icase++)
sum += xptr[icase*nvars] * x[icase*nvars+jvar] ; // Equation 3-12
Xinner[ivar*nvars+jvar] = sum / ncases ;
}
}
} // // Else not weighted
} // For ivar
Beta 优化流程概述
在前面的几节中，我们看到了对于任何选择的β权重，我们如何计算一个修正值，该值将误差标准降低到唯一的全局最小值。因此，在最简单的层面上，我们可以轮换权重，依次调整每个权重，直到获得令人满意的收敛。但是我们可以更智能地做这件事，利用这样一个事实:一旦β权重变为零，它在后续迭代中倾向于保持为零。这里显示了训练算法的概要，解释如下。稍后会出现更详细的代码。
do_active_only = 0 ; // Begin with a complete pass
for (iter=0 ; iter<maxits ; iter++) { // Main iteration loop; maxits is for safety only
active_set_changed = 0 ; // Did any betas go to/from 0.0?
for (ivar=0 ; ivar<nvars ; ivar++) { // Descend on this beta
if (do_active_only && beta[ivar] == 0.0)
continue ;
[ Compute correction ]
if (correction != 0.0) { // Did this beta change?
if ((beta[ivar]==0.0 && new_beta != 0.0) || (beta[ivar] != 0.0 && new_beta==0.0))
active_set_changed = 1 ;
}
} // For all variables; a complete pass
converged = [ Convergence test ] ;
if (do_active_only) { // Are we iterating on the active set only?
if (converged) // If we converged
do_active_only = 0 ; // We now do a complete pass
}
else { // We just did a complete pass (all variables)
if (converged && ! active_set_changed)
break ;
do_active_only = 1 ; // We now do an active-only pass
}
} // Outer loop iterations
这种训练算法的基本思想是，我们可以通过将大部分精力集中在那些非零的β权重(称为活动集)上来节省大量计算工作。粗略地说，我们通过所有的预测，调整每个β权重。过了这一关之后，经常会出现这样的情况，一些测试，也许是许多测试，是零。因此，我们进行额外的传递，只调整那些非零的(活动集)，直到获得收敛。当我们收敛时，我们会遍历所有预测值，以防修正的β权重导致一个或多个β变为零或从零开始。如果没有这样的变化发生，并且我们通过了收敛测试，我们就完成了。否则，我们将返回到仅在活动集中循环。
我们从do_active_only 假开始，以便调整所有预测值。为了安全起见，主迭代循环受到maxits的限制，尽管实际上这个限制永远不会达到。我们使用active_set_changed来标记是否有任何 beta 权重变为零或从零开始。
ivar循环对所有预测器进行一次遍历。如果我们只做活动集，而这个β是零，跳过它。否则，我们计算一个修正的β。如果 beta 发生了变化，我们可以看到变化是从零开始还是到零，如果是这样，我们可以通过设置active_set_changed标志来记录。
在我们通过预测器之后，我们执行一个收敛测试。如果我们一直只检查活动集，并且如果我们已经收敛，我们重置do_active_only以便下一次我们检查所有预测器。
另一方面，如果我们的最后一遍是对所有预测值的完整检查，获得了收敛，并且活动集没有改变，那么我们就都完成了。否则，我们设置do_active_only标志，以便我们回到只关注活动集。
这种专注于活动集的奇特算法只有在有大量零 beta 权重时才有优势。但是，在使用这种模型的应用程序中，情况往往就是这样。此外，在贝塔系数很少或没有为零的情况下，惩罚很少或没有，所以我们不妨使用花哨的版本。
Beta 优化代码
前一节给出了 beta 优化算法的概要，省略了细节，以便过程的基本逻辑将是清楚的。在本节中，我们将详细介绍整个优化代码。它的名称如下:
void CoordinateDescent::core_train (
double alpha , // User-specified alpha (0-1) (0 problem for descending lambda)
double lambda , // Can be user-specified, but usually from lambda_train()
int maxits , // Maximum iterations, for safety only
double eps , // Convergence criterion, typically 1.e-5 or so
int fast_test , // Convergence via max beta change vs explained variance?
int warm_start // Start from existing beta, rather than zero?
)
alpha ( α )和lambda (λ)参数已经见过很多次了。我们使用maxits只是为了限制迭代次数，以防止意外挂起。实际上，它会被设置得非常大。eps参数控制收敛信号发出前结果的精确度。fast_test参数控制使用两个收敛测试(稍后描述)中的哪一个。最后，warm_start允许从 beta 权重的当前值开始训练，而不是从零开始(默认)。该例程从一些初始化开始。
S_threshold = alpha * lambda ; // Threshold for the soft-thresholding S() of Eq (3-6)
do_active_only = 0 ; // Begin with a complete pass
prior_crit = 1.0e60 ; // For convergence test
if (warm_start) { // Pick up with current betas?
if (! covar_updates) { // If not using covar updates, must recompute residuals
for (icase=0 ; icase<ncases ; icase++) {
xptr = x + icase * nvars ;
sum = 0.0 ;
for (ivar=0 ; ivar<nvars ; ivar++)
sum += beta[ivar] * xptr[ivar] ;
resid[icase] = y[icase] - sum ;
}
}
}
else { // Not warm start, so initial betas are all zero
for (i=0 ; i<nvars ; i++)
beta[i] = 0.0 ;
for (i=0 ; i<ncases ; i++) // Initial residuals are just the Y variable
resid[i] = y[i] ;
}
先前初始化代码最值得注意的方面是，如果我们正在进行热启动，并且我们没有使用快速协方差更新方法，那么我们必须重新计算残差。回想一下，方程 3-7 和 3-9 的简单更新方法需要残差。当然，如果我们从零开始所有的β权重，那么所有的预测也是零，残差只是目标。
随着迭代的进展，我们将计算解释的目标方差的分数，以启迪用户。为此，我们需要目标的均方差，如果用户选择了按重要性对案例进行加权，则需要对均方差进行适当加权。下面的代码计算这个数量:
if (w != NULL) { // We need weighted squares to evaluate explained variance
YmeanSquare = 0.0 ;
for (i=0 ; i<ncases ; i++)
YmeanSquare += w[i] * y[i] * y[i] ;
}
else
YmeanSquare = 1.0 ; // The target has been normalized to unit variance
我们现在开始主要的外部循环，迭代直到获得收敛。迭代限制maxits应该设置得非常大(几千或更多)以便它不会导致过早退出；只是“挂保险”。我们重置标志，该标志将指示活动集是否改变，并且我们将使用max_change来跟踪收敛测试的最大 beta 变化。
for (iter=0 ; iter<maxits ; iter++) {
active_set_changed = 0 ; // Did any betas go to/from 0.0?
max_change = 0.0 ; // For fast convergence test
对所有预测值进行一次遍历的循环现在开始。如果我们只处理活动集(非零 beta ),并且这个 beta 为零，跳过它。加权情况下的方程 3-9 和未加权情况下的方程 3-7 将需要分母中的update_factor，因此现在计算它。回想一下XSSvec[]是通过等式 3-13 计算出来的。
for (ivar=0 ; ivar<nvars ; ivar++) { // Descend on this beta
if (do_active_only && beta[ivar] == 0.0)
continue ;
// Denominator in update
if (w != NULL) // Weighted?
Xss = XSSvec[ivar] ;
else
Xss = 1 ; // X was standardized
update_factor = Xss + lambda * (1.0 - alpha) ;
我们计算软阈值函数的自变量。有三种可能。要么我们使用快速协方差更新方法，要么我们使用具有不同情况加权的朴素方法，要么我们使用具有相等加权的朴素方法。我们不必在这里将协方差更新方法分为有权重和无权重，因为在计算Xss、Xinner和Yinner时已经考虑了任何权重，如第 42 页所示。
if (covar_updates) { // Any sensible user will specify this unless ncases < nvars
sum = 0.0 ;
for (kvar=0 ; kvar<nvars ; kvar++)
sum += Xinner[ivar*nvars+kvar] * beta[kvar] ;
residual_sum = Yinner[ivar] - sum ;
argument = residual_sum + Xss * beta[ivar] ; // Equation 3-10
}
else if (w != NULL) { // Use slow naive formula (okay if ncases < nvars)
argument = 0.0 ;
xptr = x + ivar ; // Point to column of this variable
for (icase=0 ; icase<ncases ; icase++) // Equation 3-8
argument += w[icase] *
xptr[icase*nvars] * (resid[icase] + beta[ivar] * xptr[icase*nvars]) ;
}
else { // Use slow naive formula (okay if ncases < nvars)
residual_sum = 0.0 ;
xptr = x + ivar ; // Point to column of this variable
for (icase=0 ; icase<ncases ; icase++)
residual_sum += xptr[icase*nvars] * resid[icase] ; // X_ij * RESID_i
residual_sum /= ncases ;
argument = residual_sum + beta[ivar] ; // Equation 3-5
}
我们刚刚计算了软阈值函数的自变量，方程 3-6 。应用该函数，并使用公式 3-7 或公式 3-9 计算该 beta 的新值。不久前，我们计算出update_factor是这些方程中的分母。
if (argument > 0.0 && S_threshold < argument)
new_beta = (argument - S_threshold) / update_factor ;
else if (argument < 0.0 && S_threshold < -argument)
new_beta = (argument + S_threshold) / update_factor ;
else
new_beta = 0.0 ;
修正量是新的 beta 值和旧值之间的差值。记录这一过程中的最大变化，因为我们可能会用它来进行收敛测试。如果我们使用缓慢的朴素更新方法，我们也将使用这种校正来快速重新计算残差，这是朴素方法所需要的。
correction = new_beta - beta[ivar] ;
if (fabs(correction) > max_change)
max_change = fabs(correction) ; // Used for fast convergence test
if (correction != 0.0) { // Did this beta change?
if (! covar_updates) { // Must we update the residual vector?
xptr = x + ivar ; // Point to column of this variable
for (icase=0 ; icase<ncases ; icase++) // Update residual per this new beta
resid[icase] -= correction * xptr[icase*nvars] ;
}
if ((beta[ivar]==0.0 && new_beta!=0.0) || (beta[ivar]!=0.0 && new_beta==0.0))
active_set_changed = 1 ;
beta[ivar] = new_beta ;
}
} // For all variables; a complete pass
根据do_active_only的说法，我们已经完成了一次测试，要么是所有的测试，要么只是活动集。我们现在做收敛测试，无论是快速，简单的版本或慢得多的版本。快速测试仅基于β的最大(所有预测值)变化。但是慢速测试更复杂。
如果我们使用快速协方差更新方法，我们不需要 beta 更新的残差，所以我们不需要(巨大的！)是时候计算它们了。但是我们需要残差来进行缓慢收敛测试，所以如果到目前为止我们还没有计算它们，我们必须计算它们。使用残差计算(可能加权的)均方误差。
if (fast_test) { // Quick and simple test
if (max_change < eps)
converged = 1 ;
else
converged = 0 ;
}
else { // Slow test (change in explained variance) which requires residual
if (covar_updates) { // We have until now avoided computing residuals
for (icase=0 ; icase<ncases ; icase++) {
xptr = x + icase * nvars ;
sum = 0.0 ;
for (ivar=0 ; ivar<nvars ; ivar++)
sum += beta[ivar] * xptr[ivar] ; // Cumulate predicted value
resid[icase] = y[icase] - sum ; // Residual = true - predicted
}
}
sum = 0.0 ; // Will cumulate squared error for convergence test
if (w != NULL) { // Are the errors of each case weighted differently?
for (icase=0 ; icase<ncases ; icase++)
sum += w[icase] * resid[icase] * resid[icase] ;
crit = sum ;
}
else {
for (i=0 ; i<ncases ; i++)
sum += resid[i] * resid[i] ;
crit = sum / ncases ; // MSE component of optimization criterion
}
模型的基本质量度量是模型解释的目标方差的分数。这是通过从目标的均方(方差)中减去刚刚计算的均方误差来计算的，以获得所解释的方差的量。将其除以目标均方差，得到模型所解释的目标方差的分数。这是严格用于可选的用户启发；它在优化算法中不起作用。
使用第 39 页的等式 3-3 计算正则化罚分，然后将该罚分与均方误差相加，得到我们最小化的标准，如第 38 页的等式 3-2 所示。
这个“慢”收敛标准是基于优化标准中从一次迭代到下一次迭代的变化。如果变化很小(其中“小”是由用户指定的eps定义的)，那么我们被认为已经收敛。
explained_variance = (YmeanSquare - crit) / YmeanSquare ;
penalty = 0.0 ;
for (i=0 ; i<nvars ; i++)
penalty += 0.5 * (1.0 - alpha) * beta[i] * beta[i] + alpha * fabs (beta[i]) ;
penalty *= 2.0 * lambda ; // Regularization component of optimization criterion
crit += penalty ; // This is what we are minimizing
if (prior_crit - crit < eps)
converged = 1 ;
else
converged = 0 ;
prior_crit = crit ;
}
现在，我们可以使用上一节中描述的控制逻辑完成外循环，在仅活动集和完整预测器通道之间交替。
if (do_active_only) { // Are we iterating on the active set only?
if (converged) // If we converged
do_active_only = 0 ; // We now do a complete pass
}
else { // We just did a complete pass (all variables)
if (converged && ! active_set_changed)
break ;
do_active_only = 1 ; // We now do an active-only pass
}
} // Outer loop iterations
我们基本上完成了。为了启发用户，我们计算并保存由模型解释的目标方差的分数。如果我们做了快速收敛测试和协方差更新，我们必须计算残差来得到解释的方差。这两个选项不需要常规的残差计算，所以我们目前没有残差。
if (fast_test && covar_updates) { // Residuals have not been maintained?
for (icase=0 ; icase<ncases ; icase++) {
xptr = x + icase * nvars ;
sum = 0.0 ;
for (ivar=0 ; ivar<nvars ; ivar++)
sum += beta[ivar] * xptr[ivar] ;
resid[icase] = y[icase] - sum ;
}
}
sum = 0.0 ;
if (w != NULL) { // Error term of each case weighted differentially?
for (i=0 ; i<ncases ; i++)
sum += w[i] * resid[i] * resid[i] ;
crit = sum ;
}
else {
for (i=0 ; i<ncases ; i++)
sum += resid[i] * resid[i] ;
crit = sum / ncases ; // MSE component of optimization criterion
}
explained = (YmeanSquare - crit) / YmeanSquare ;
沿λ路径下行
如同通常具有超参数的模型的情况一样，为正则化强度λ(λ)选择有效值可能并不简单。在下一节中，我们将探索一种自动选择好值的好方法。在本节中，我们将介绍一个工具，该工具将由自动化例程调用，并且还可以用于帮助手动选择一个好的 lambda。
考虑一下，如果 lambda 很大，那么任何非零 beta 的惩罚都将很大，以至于所有 beta 权重都将被强制为零。(如果α正好为零，情况可能不是这样，所以从现在开始我们将假设α为 0。)这个模型显然是零解释方差。相反，如果λ = 0，那么我们有普通的线性回归，它具有最小可能的均方误差或最大可能的解释方差。因此，我们可以从一个较大的λ开始，训练模型，稍微降低λ，然后再次训练，以此类推，直到λ很小，几乎为零。我们通常会看到非零贝塔的数量稳步增加，以及解释方差。即使对于相同数量的非零贝塔，解释的方差也会随着λ的减小而增加。如果我们打印一张图表，显示非零贝塔值的数量和解释的方差作为λ的函数，我们也许能够对λ作出明智的选择。
这种方法有一个有趣的额外好处，即使我们事先知道我们想要使用的λ。这种方法增加了训练算法已经相当好的稳定性，而在速度方面没有太大的代价。事实上，我们可以这样训练更快。我们所做的是从一个大的 lambda 开始，这个 lambda 只给我们一个或很少几个有效的预测器。那种简单的模型会训练得很快。然后，当我们稍微降低 lambda 时，我们不是从头开始，而是热启动，从现有的 betas 开始迭代。所以，每次我们用稍小的λ重新开始训练，我们都是从已经非常接近正确的 betas 开始。因此，收敛将很快获得。
很容易为下降找到一个好的起始λ，最小的λ使得所有的贝塔为零。整个过程从所有 betas 为零开始。查看方程 3-7 ，以及自变量和软阈值算子的两个先前方程。对于差分加权的情况，在下一节中查看它们的类似物。回想一下，当所有贝塔系数都为零时，残差等于目标值， y 。根据软阈值函数的定义，很明显，如果未加权情况下的等式 3-16 或差分加权情况下的等式 3-17 为真，则 β 将保持为零。
)
(3-16)
)
(3-17)
将这些方程的两边除以α得到任何预测因子的阈值λ，如果我们找到所有预测因子中的最大λ，我们就有了起始λ。下面是执行此操作的代码:
double CoordinateDescent::get_lambda_thresh ( double alpha )
{
int ivar, icase ;
double thresh, sum, *xptr ;
thresh = 0.0 ;
for (ivar=0 ; ivar<nvars ; ivar++) {
xptr = x + ivar ;
sum = 0.0 ;
if (w != NULL) {
for (icase=0 ; icase<ncases ; icase++) // Left side of Equation 3-17
sum += w[icase] * xptr[icase*nvars] * y[icase] ;
}
else {
for (icase=0 ; icase<ncases ; icase++) // Left side of Equation 3-16
sum += xptr[icase*nvars] * y[icase] ;
sum /= ncases ;
}
sum = fabs(sum) ;
if (sum > thresh) // We must cover all predictors
thresh = sum ;
}
return thresh / (alpha + 1.e-60) ; // Solve for lambda; protect from division by zero
}
在 lambda 上下降很简单。需要注意的一点是，我们保存了每个试用 lambda 的 beta 权重，因为我们以后可能需要访问它们。此外，如果调用者设置了print_steps标志，这个例程将打开一个文本文件并附加结果，以便于用户检查。
我们使用get_lambda_thresh()找到最小的λ，确保所有的贝塔值保持为零，并稍微减小它以得到我们的起始λ。我们任意地将最小λ设置为该量的 0.001 倍。在构造函数调用中指定了试验次数。代码如下:
void CoordinateDescent::lambda_train (
double alpha , // User-specified alpha, (0,1) (Greater than 0)
int maxits , // Maximum iterations, for safety only
double eps , // Convergence criterion, typically 1.e-5 or so
int fast_test , // Convergence via max beta change vs explained variance?
double max_lambda , // Starting lambda, or negative for automatic computation
int print_steps // Print lambda/explained table?
)
{
int ivar, ilambda, n_active ;
double lambda, min_lambda, lambda_factor ;
FILE *fp_results ;
if (print_steps) {
fopen_s ( &fp_results , "CDtest.LOG" , "at" ) ;
fprintf ( fp_results , "\n\nDescending lambda training..." ) ;
fclose ( fp_results ) ;
}
if (n_lambda <= 1) // Nonsensical parameter from caller
ireturn ;
/*
Compute the minimum lambda for which all beta weights remain at zero
This (slightly decreased) will be the lambda from which we start our descent.
*/
if (max_lambda <= 0.0)
max_lambda = 0.999 * get_lambda_thresh ( alpha ) ;
min_lambda = 0.001 * max_lambda ;
lambda_factor = exp ( log ( min_lambda / max_lambda ) / (n_lambda-1) ) ;
/*
Repeatedly train with decreasing lambdas
*/
if (print_steps) {
fopen_s ( &fp_results , "CDtest.LOG" , "at" ) ;
fprintf ( fp_results , "\nLambda n_active Explained" ) ;
}
lambda = max_lambda ;
for (ilambda=0 ; ilambda<n_lambda ; ilambda++) {
lambdas[ilambda] = lambda ; // Save in case we want to use later
core_train ( alpha , lambda , maxits , eps , fast_test , ilambda ) ;
for (ivar=0 ; ivar<nvars ; ivar++) // Save these in case we want them later
lambda_beta[ilambda*nvars+ivar] = beta[ivar] ;
if (print_steps) {
n_active = 0 ; // Count active predictors for user’s edification
for (ivar=0 ; ivar<nvars ; ivar++) {
if (fabs(beta[ivar]) > 0.0)
++n_active ;
}
fprintf ( fp_results , "\n%8.4lf %4d %12.4lf", lambda, n_active, explained ) ;
}
lambda *= lambda_factor ;
}
if (print_steps)
fclose ( fp_results ) ;
}
通过交叉验证优化 Lambda
如果不是最流行的，也是最流行的优化模型超参数的方法之一是交叉验证，所以这就是我们在这里要做的。原理很简单。对于每个折叠，我们调用lambda_train()来测试一组递减的λ，保存每个试验λ的β系数。然后，我们计算每个试验λ的样本外解释方差，并累积该量。当所有折叠完成后，我们检查合并的 OOS 性能，并选择哪个λ给出最佳 OOS 性能。不过，有一些事情需要注意，所以我们将把这段代码分成几个单独的部分，分别进行解释。下面是调用参数列表:
double cv_train (
int n , // Number of cases in full database
int nvars , // Number of variables (columns in database)
int nfolds , // Number of folds
double *xx , // Full database (n rows, nvars columns)
double *yy , // Predicted variable vector, n long
double *ww , // Optional weights, n long, or NULL if no weighting
double *lambdas , // Returns lambdas tested by lambda_train()
double *lambda_OOS , // Returns OOS explained for each of above lambdas
int covar_updates , // Does user want (usually faster) covariance update method?
int n_lambda , // This many lambdas tested by lambda_train() (at least 2)
double alpha , // User-specified alpha, (0,1) (greater than 0)
int maxits , // Maximum iterations, for safety only
double eps , // Convergence criterion, typically 1.e-5 or so
int fast_test // Convergence via max beta change vs explained variance?
)
注意，这不是CoordinateDescent类的成员；这是一个独立的程序。大多数参数都是不言自明的，以前也见过很多次。最后四个参数和covar_updates只是传递给核心训练例程。我们必须提供两个向量n_lambdas long: lambdas将返回测试的λ值，lambda_OOS将返回对应于每个测试的λ的 OOS 解释的方差分数。我们应该指定尽可能大的n_lambdas来进行彻底的测试；50 不是不合理的。大量的 lambda 不会明显降低训练速度，因为使用了热启动，这意味着每次 lambda 降低时，beta 优化都在先前的最优值处开始。这非常快。最后，为了获得最佳精度，折叠次数也应该尽可能多；五个是最低限度，十个是合理的，如果计算机时间允许，更多更好。
我们从一些初始化开始。自然地，我们希望对每个折叠使用相同的降序 lambdas 集合，因此我们使用整个数据集来寻找阈值。如果案例被加权，我们复制归一化的权重用于 OOS 评分。第一次培训将从第一个案例开始，我们还没有做过任何 OOS 案例。我们将累积在lambda_OOS中解释的方差的分数，因此对于每个试验λ将这个向量初始化为零。我们将累积YsumSquares中的(可能加权的)目标平方和。
cd = new CoordinateDescent ( nvars , n , (ww != NULL) , covar_updates , n_lambda ) ;
cd->get_data ( 0 , n , xx , yy , ww ) ; // Fetch the training set for this fold
max_lambda = cd->get_lambda_thresh ( alpha ) ;
if (ww != NULL) {
for (icase=0 ; icase<n ; icase++)
work[icase] = cd->w[icase] ;
}
delete cd ;
i_IS = 0 ; // Training data starts at this index in complete database
n_done = 0 ; // Number of cases treated as OOS so far
for (ilambda=0 ; ilambda<n_lambda ; ilambda++)
lambda_OOS[ilambda] = 0.0 ; // Will cumulate across folds here
YsumSquares = 0.0 ; // Will cumulate to compute explained fraction
折叠循环从这里开始。OOS 案例的数量是剩余要做的数量除以剩余折叠的数量。其余情况是样本内的，OOS 集从 IS 集之后开始。
for (ifold=0 ; ifold<nfolds ; ifold++) {
n_OOS = (n - n_done) / (nfolds - ifold) ; // Number of cases in OOS (test set)
n_IS = n - n_OOS ; // Number IS (training set)
i_OOS = (i_IS + n_IS) % n ; // OOS starts at this index
我们现在用这个样本集训练，在 lambda 上下降。这个集合从索引i_IS开始，如果到达数据集的结尾，它将循环回到开始。
cd = new CoordinateDescent ( nvars , n_IS , (ww != NULL) , covar_updates ,
n_lam bda ) ;
cd->get_data ( i_IS , n , xx , yy , ww ) ; // Fetch the training set for this fold
cd->lambda_train ( alpha , maxits , eps , fast_test , max_lambda , 0 ) ;
训练已经完成，所以我们在 OOS 集上评估性能。下面是代码；下一页是逐步解释:
for (ilambda=0 ; ilambda<n_lambda ; ilambda++) {
lambdas[ilambda] = cd->lambdas[ilambda] ; // This will be the same for all folds
coefs = cd->lambda_beta + ilambda * nvars ;
sum = 0.0 ;
for (icase=0 ; icase<n_OOS ; icase++) {
k = (icase + i_OOS) % n ;
pred = 0.0 ;
for (ivar=0 ; ivar<nvars ; ivar++)
pred += coefs[ivar] * (xx[k*nvars+ivar] - cd->Xmeans[ivar]) / cd->Xscales[ivar] ;
Ynormalized = (yy[k] - cd->Ymean) / cd->Yscale ;
diff = Ynormalized - pred ;
if (ww != NULL) {
if (ilambda == 0)
YsumSquares += work[k] * Ynormalized * Ynormalized ;
sum += work[k] * diff * diff ;
}
else {
if (ilambda == 0)
YsumSquares += Ynormalized * Ynormalized ;
sum += diff * diff ;
}
}
lambda_OOS[ilambda] += sum ; // Cumulate for this fold
} // For ilambda
delete cd ;
n_done += n_OOS ; // Cumulate OOS cases just processed
i_IS = (i_IS + n_OOS) % n ; // Next IS starts at this index
} // For ifold
上一页的代码处理单个折叠的 OOS 集。随着有效的λ下降算法的进展，训练例程保存每个试验λ的β权重。因此，我们遍历 lambda，将每个 lambda 的 betas 放入coefs。我们将遍历所有 OOS 案例，累计sum中的误差平方和。
我们循环遍历数据集，当到达末尾时循环回到开始，因此k是将要测试的 OOS 案例的索引。OOS 情况(目标和所有预测值)必须以与训练数据相同的方式进行归一化，使用相同的平均值和标准差。
这种情况下的误差diff，是真实值减去预测值。我们累积平方误差，如果使用了不同的权重，则乘以用户指定的案例权重。我们同时累积归一化目标的平方和。这必须只做一次，因为它当然对所有试验 lambdas 都是一样的。当 case 循环完成时，我们将误差总和加到被测试的 lambda 的总和上。这个向量将累积所有折叠的总和。lambda 循环完成后，我们删除这个文件夹的CoordinateDescent对象，并前进到下一个文件夹。
剩下要做的就是计算每个 lambda 的 OOS 解释方差分数，并将表现最好的 lambda 返回给调用者。目标平方和减去误差平方和得到解释的平方和。将其除以目标 SS，得到解释方差的分数。
best = -1.e60 ;
for (ilambda=0 ; ilambda<n_lambda ; ilambda++) {
lambda_OOS[ilambda] = (YsumSquares - lambda_OOS[ilambda]) / YsumSquares ;
if (lambda_OOS[ilambda] > best) {
best = lambda_OOS[ilambda] ;
ibest = ilambda ;
}
}
return lambdas[ibest] ;
}
CD_MA 程序
文件 CD_MA。CPP 包含一个程序，该程序读取市场价格文件，根据移动平均线振荡器计算大量指标，并使用CoordinateDescent正则化线性模型来找到指标的最佳子集，以预测第二天的(对数)价格变化。历史文件末尾的一年的市场数据将作为测试集使用。
使用以下命令调用该程序:
CD_MA Lookback_inc N_long N_short Alpha Filename
让我们来分解这个命令:
-
Lookback_inc:长期回看将从这个数量的棒线(包括当前棒线)开始回看。后续的长期回顾将以此数量递增。例如，如果指定为 3，则长期回看将为 3、6、9，… -
这么多的长期回顾将会被采用。最大长期回顾将是
Lookback_inc * N_long。 -
这么多的短期回顾将会被采用。它们是当前的长期回顾时间 i 然后除以
N_short+1，对于从 1 到N_short的 i ，被截断成整数。注意，当当前长期回顾小于N_short+1时，将会有多个相同的短期回顾值，从而产生完全相关的预测值。指标总数为N_long*N_short。 -
Alpha:控制正则化类型所需的 alpha。如果指定小于或等于零，lambda 将被设置为零，产生普通线性回归(无正则化)。它决不能大于或等于 1。 -
Filename:格式为YYYYMMDD Open High Low Close的市场历史文件。
将打印两个表格。第一部分显示了选择最佳λ所涉及的计算。此表中的左栏列出了试用的 lambdas。右栏显示解释方差的相应样本外分数。
第二个表列出了β系数。每一行对应一个长期回顾，回顾打印在每一行的开头。每根柱对应一个短期回看。这些回看不会被打印，因为它们随每行而变化。它们可以用前页的公式很容易地计算出来。精确为零的系数(通常但不总是因为训练算法将它们从模型中移除)用虚线表示。
图 3-1 显示了当λ= 0，无正则化时，为 OEX 产生的β系数表。这实际上等同于普通的线性回归。图 3-2 显示 alpha=0.1 时的结果，图 3-3 针对 alpha=0.9。讨论如下。
图 3-3
阿尔法=0.9
图 3-2
阿尔法=0.1
图 3-1
λ= 0(无正则化)
这次运行使用 OEX S&P 100 指数作为它的市场历史。回望增量为 2，有 30 次长期回望和 10 次短期回望。
-
回想一下本节第一页的讨论，对于小于短期回顾数加 1 的长期回顾，一些短期回顾必须重复，这意味着一些指标是其他指标的精确副本。
-
这种重复使得普通的线性回归变得不可能，因为一些权重是未定义的。将需要诸如奇异值分解的特殊技术。这里 lambda=0 的算法处理得很好，甚至有效地消除了一些重复。但是绝大多数指标都参与了这个模型。
-
因为当λ= 0 时，没有正则化，这是完全最小二乘拟合。这意味着解释方差的样本内部分应该是最大可能值，事实上我们看到了这种情况，解释了 1.63%的目标方差。
-
由于大量的指标参与(没有正规化)，我们预计会看到 OOS 表现不佳。是的，这是三次测试中得分最差的一次。
-
当我们应用 alpha=0.1 的正则化(近似岭回归)时，样本内解释方差下降，但 OOS 性能飙升至最佳。
-
在 alpha=0.1 的情况下，我们看到重复的指示符接收到相同的 beta 系数，如预期的那样。
-
当 alpha=0.9(接近套索)时，模型会最小化保留的指标数量，以尝试使模型尽可能简单，即使以牺牲性能为代价。我们看到这种情况发生，甚至选择的指标也发生了变化。OOS 业绩暴跌，意味着该模型被迫放弃一些有用的指标。
-
正则化的模型都是负系数，说明这个交易系统是均值回归系统，不是趋势跟随者！
使线性模型变得非线性
尽管线性模型往往比非线性模型更受青睐，但有时我们的两个或更多指标在与目标的关系中存在不可避免的非线性相互作用。要知道，一个指标仅仅与目标有一个非线性关系，就其本身而言，通常不是问题。我们可以通过某种方式改变指标，使其与目标之间的关系变得基本上是线性的。至少尝试一下总是件好事。当然，也可能发生这样的情况，我们只是怀疑一个单独的非线性关系，但我们无法证明它足以能够合理地转换指标。但是在绝大多数情况下，当指标在与目标的联合关系方面以非线性方式相互作用时，线性模型就会失效。在这种情况下，我们别无选择，只能放弃严格的线性模型。
但是并没有失去一切。线性模型的优势，尤其是正则化模型的优势(对其工作原理的简单理解、快速训练、较低的过拟合可能性)如此之大，以至于值得以适度非线性的方式转换指标及其相互作用，并将这些新值应用于正则化线性模型。我们几乎从来不想采用如此极端的措施，以至于交易决策的界限到处游走，为了抓住每一个错误的训练案例而扭曲。但是有一种简单的方法来应用适度的非线性变换，允许我们以温和的非线性方式使用正则化的线性模型。
自然，我们可以用一个或多个原始预测值的一个或多个非线性函数来补充模型的预测值。如果我们有选择某些特定函数的理论理由，我们当然应该这样做。但那种情况很少见。最常见和最有效的一般程序是使用低次多项式，我将很快讨论两个特殊的扭曲。总体思路是这样的:我们选择低学位，一般是二级，很少是三级。此外，选择一个子集的预测，我们希望允许非线性相互作用。这可能是所有的预测因素，尽管当我们包含更多的预测因素时，事情会变得更糟。然后在选定的程度上，用它们的每一种可能的组合来补充原始预测值。
例如，假设我们有三个预测值，我们希望考虑非线性。称他们为 A、B 和 C 。还假设我们想允许二级，最常见的选择。那么我们发给模型的预测器是 A 、 B 、 C 、 A 2 、 B 2 、 C 2 、 AB 、 AC 、 BC 。如果我们决定向上移动到第三级，附加的预测因子是 A 3 、 B 3 、 C 3 、A2B、A2C、B 显而易见，增加非线性集合中预测因子的数量，或者增加多项式的次数，会导致新预测因子数量的爆炸性增长。
当使用多项式展开时，有两件事应该做。这两者都不是数学上必需的，但如果我们要防止硬件浮点不准确性，以及提高大多数模型训练算法的速度和稳定性，这两者都很重要。首先，我们必须确保转换后的指标具有大约负一比一的自然范围。如果做到了这一点，所有的多项式变换值都具有相同的自然范围。如果我们的原始指标没有这个范围，至少近似地，我们应该从理论考虑或从大的代表性集合的检查中找到它们的真实自然范围，最小到最大。则 X 的量程调整值为 2 *(??)X–Min)/(Max–Min)–1。
我们应该采取的另一个行动只有在我们达到第三级(或者，但愿不会如此，更高级)时才需要。问题是，即使有范围调整， X 和 X 3 可以有足够的相关性，稍微妨碍一些训练算法。它很少是严重的，并且将要描述的技术可能会被一些人认为是过度的，但是它是一个具有良好回报的廉价投资。不使用X3，而是使用 0.5 (5X3–3X*)。这仍然是一个范围为负一比一的三次多项式，它将允许与X3 相同的有效非线性，但它通常与 X 的相关性要小得多，因此将由许多训练算法更有效地处理。你不会有任何损失，而且可能会有很大收获。
超越三级几乎总是毫无意义的。如果你有那么多的非线性，就用非线性模型。但是如果出于某种原因你坚持的话，查查勒让德多项式并用它们来表示高次项。
差分进化:一种通用的非线性优化器
无论你的交易系统是基于非线性预测模型还是传统的算法(基于规则)系统，你都需要一个快速稳定的方法来优化你选择的性能标准。在大多数情况下，优化算法的选择涉及到一个重要的权衡。一个多变量函数通常有几个(也许很多！)局部 optima。最快的优化者是爬山者，迅速爬到最近的山顶，不管那个特定的山是不是最好的。更有可能在众多山顶中找到最佳山顶的优化者比简单的爬山者要慢得多。那么，该怎么办呢？
幸运的是，有一种算法是两个极端之间的一个很好的妥协。在众多山顶中，它找到最好的，或者至少接近最好的山顶的可能性相对较高，但它的速度也相当快。这种算法是一种特殊的遗传或进化优化，称为差分进化。我不会在这里提供任何参考，因为互联网上充满了例子和讨论。相反，我将把重点放在这个算法的一个高度调整的变体上，我已经在我自己的工作中使用了多年，并且我发现它是一个可靠的执行者。
像所有的进化算法一样，它从一群个体开始，每个个体都是交易系统的一个完全指定的参数集。然后，它遍历整个种群，以一种很有可能产生优于双亲的个体的方式组合种群中不同成员的品质。
不同于大多数植物和动物的繁殖，差异进化需要四个个体才能产生一个孩子。其中一个被称为父代 1 ，是确定性选择的。另外三个父母 2 、差异 1 和差异 2 是随机选择的。两个差值决定了一个方向，并且 Parent2 在这个方向上被扰动。通过从 Parent1 中选择一些参数以及从受干扰的 Parent2 中选择其他参数来创建新子节点。如图 3-4 和图 3-5 所示。
图 3-5
差异子代
图 3-4
差异进化的一步
图 3-4 显示一阶和二阶微分相减，它们的差乘以一个常数，通常小于 1。这个缩小的差被加到第二个亲本上，和在随机交叉操作中与第一个亲本合并。将这个子代的表现与主要父代的表现进行比较，并为下一代保留优越的个体。
这在图 3-5 中用图形显示了两个变量。两个微分之间的差决定了一个方向，次级母体在这个方向上被扰动。在这个例子中，这个操作被称为突变，尽管这不是一个通用术语。然后水平变量取自第一亲本，垂直变量取自变异的第二亲本。
这个方案有一个重要的性质:它将扰动缩放到参数的自然尺度。假设功能函数在某个方向上有一个狭窄的脊，这是一种常见的情况。那么人口会被吸引到同样的布局。个体(完整的参数集)将在脊的方向上广泛分布，而在垂直方向上被压缩。其结果是，控制次级母体扰动程度的微分差将沿着脊变大，沿着脊变小，这正是我们想要的。
不幸的是，差分进化有一个大多数随机过程共有的弱点:它很快收敛到全局最优附近，但从来没有完全达到精确的最优。这是因为它本质上无法利用函数的本地知识。爬山法在收敛到局部最优方面做得很好，但是它们容易因为错过全局最优而失败。所以，我的方法是混合的，主要实现差异进化，但偶尔对单个个体执行爬山步骤。这极大地加速了收敛，同时对算法的全局影响可以忽略，因为该操作在任何时候都保持在单个个体的吸引域内。
此处显示了下一页的算法概述。
for (ind=0 ; ind<popsize+overinit ; ind++) { // Generate the basis population
Generate a random parameter vector
value = performance of this random vector
If this individual fails to meet a minimum requirement {
--ind ; // Skip it entirely
continue ;
}
if (ind >= popsize) { // If we finished population, now doing overinit
Find the worst individual in population (popsize individuals)
if (value > worst)
Replace the worst with this new individual
} // If doing overinit
} // For all individuals (population and overinit)
for (generation=1 ; ; generation++) {
for (ind=0 ; ind<popsize ; ind++) { // Generate all children
parent1 = individual 'ind'
Generate three different random indices for parent2 and the two differentials:
parent2, diff1, diff2
for all variables j { // This is the mutation and crossover step
with small probability
trial[j] = parent2[j] + constant * (diff1[j] - diff2[j]) ;
else
trial[j] = parent1[j] ;
}
value = performance of trial parameter set
if (value > parent1's performance)
replace parent1 with trial parameter set
Optionally pick one variable in one individual (favoring the best so far)
and find the optimal value of that variable
} // Create all children in this generation
} // For all generations
Return the best individual in the final population
第一步是生成一个初始群体，这是由这段伪代码中的第一个循环完成的。用户指定群体中的个体数量popsize。传统算法不包括过度初始化，测试overinit额外的个体。我发现，将overinit设置为大约popsize会以相对较小的额外成本产生一个具有更快收敛速度和更好全局代表性的显著更优的初始种群。
这个群体生成循环创建一个随机个体(完整的交易系统参数集)并计算其性能。如果这个人没有满足任何用户指定的要求，比如最低交易次数，他就会被拒绝，我们会再试一次。
当我们生成了popsize个个体并进入过度初始化时，对于每个新的候选人，我们在现有的群体中搜索表现最差的个体。如果新的候选人优于群体中最差的个体，则新的候选人取代最差的个体。这稳定地提高了群体的质量，也使我们更有可能在全局最优的吸引域中有一个或多个个体。
然后我们来看代码的进化部分。有两个嵌套循环。外部循环处理世代，在每一代中，我们将当前群体中的每个个体作为主要父代。次要亲本和两个差异个体是随机选择的，当然这四个个体必须是不同的。
为了提高效率，变异(通过缩小差异来干扰次级亲本)和交叉(用相应的变异变量随机替换初级亲本中的一些变量)在同一循环中完成。我们遍历所有变量来创建一个试验个体。对于每一个，我们掷骰子，通常以很小的概率将变量设置为变异值。否则，我们从主父节点复制变量。
我们计算这个试验个体的表现。如果它优于原始亲本，它将进入下一代群体。否则，主要的父代会延续到下一代。
最后，我们可选地执行传统算法中没有出现的步骤，但是我发现该步骤在加速收敛方面是有用的，同时在存在多个次优局部最优解的情况下对算法找到全局最优解的能力影响很小或者没有影响。我们在人群中挑选一个个体，对当前最好的个体有所偏爱，我们也挑选一个变量。我们使用爬山算法来找到这个变量的值，以优化这个个体的性能。这给了我们两个世界(随机优化与爬山)的最佳选择，因为它让算法准确地收敛到精确的最优值，比纯粹的随机算法快得多，同时它不干扰差分进化找到全局最优值的能力。这是因为当它完成时，它只发生在一个个体上，这使该个体保持在其局部最优的吸引域内，而不触及群体中可能具有其自己的吸引域的其他个体。因此，吸引域保持分离，使全球最佳者最终占据主导地位。
在所有代完成后，我们选择最终种群中最好的个体，并将其返回给用户。
为了清楚起见，刚刚显示的算法被简化了。我的实现要复杂得多，因为这些年来我已经用许多方法对它进行了改进，以调整它的性能，尤其是在优化交易系统的情况下。从下一页开始，我们将完成整个子程序，分别列出和注释每个部分。这个代码可以在 DIFF_EV.CPP 文件中找到。注意，这个文件还包含一些与差异进化无关但同时执行起来很有效的其他代码。我们将在这里忽略这段代码，在第 91 页详细讨论。
DIFF_EV 差异进化的 CPP 程序
使用以下参数列表调用差分进化子例程:
int diff_ev (
double (*criter) ( double * , int ) , // Crit function maximized
int nvars , // Number of variables (trading system parameters)
int nints , // Number of first variables that are integers
int popsize , // Population size
int overinit , // Overinitialization for initial population
int mintrades , // Minimum number of trades for candidate system
int max_evals , // For safety, max number of failed initial performance evals
int max_bad_gen , // Max number of contiguous gens with no improvement
double mutate_dev , // Deviation for differential mutation
double pcross , // Probability of crossover
double pclimb , // Probability of taking a hill-climbing step, can be zero
double *low_bounds , // Lower bounds for parameters
double *high_bounds , // And upper
double *params , // Returns nvars best parameters, plus criterion at end
int print_progress // Print progress to screen?
)
调用者提供的criter()函数计算交易系统的性能标准，该标准将被最大化。它采用交易系统可优化参数的向量。在我的实现中，还提供的整数是用户指定的最小交易数。读者应该会发现，在为生成的交易系统设置最低要求时，很容易添加其他变量。
参数可以是整数或实数；正如将要看到的，它们在内部被不同地处理。所有整数参数必须在参数数组中排在第一位，nints指定整数的个数。
用户可以将overinit设置为零，以使用传统版本的算法。然而，我发现将其设置为等于popsize附近的值是有利的。这有助于加速收敛并增加找到真正全局最大值的概率。但请注意，收益递减点很快就会达到。很快就会发生这样的情况，稳定改进的群体中最差的个体通常优于大多数过度初始化的个体，使得继续过度初始化是一种浪费。
用户在mintrd中指定所需交易的最小数量。从代码演示中可以看出，如果优化器很难找到满足这个要求的系统，那么指定的数量可能会自动减少。因此，用户应该检查通过最佳系统获得的交易数量，以确认它是令人满意的。如果程序员愿意，消除这种自动减少是很容易的，但我发现它很有用。
max_evals参数是一种安全措施。如果交易系统本质上很差，以至于大多数试验参数都产生了被拒绝的系统，那么产生初始群体就要花费大量的时间。为了防止这种情况，请将max_evals设置为一个较大但合理的值。这不应该被认为是一个收敛测试；实际上，这个极限应该永远不会遇到。
收敛由max_bad_gen参数定义。如果这许多连续的代在最佳个体中没有改进，则获得收敛，并且算法停止。这通常应该是相当大的，也许 50 甚至更多，因为事情可能会因为运气不好而变得糟糕，然后突然再次起飞。
当一个突变的参数替换了主亲本中的相应参数时，就会发生一个交叉，这种情况发生的概率由pcross给出。这通常应该很小，可能最多为 0.1 到 0.5。
爬山步骤的概率由pclimb给出。这可以是零，以严格避免爬山，保持传统版本的差异进化。它可以被设置为一个很小的正值，比如 0.00001，在这种情况下，当前最好的个体(没有其他人)偶尔会爬山。这最大限度地提高了末期精确收敛。最后，可以将其设置为稍大但仍然较小的值，如 0.2。这样，除了调整最佳个体，它还会偶尔随机调整其他个体。将它设置为较大的值通常不是很有利，因为爬山是一种昂贵的操作，尤其是对于实参数，偶尔多做一次的回报通常不足以补偿增加的成本。此外，如果爬山太频繁，对真正的全局最大值的检测可能会有些受阻，尽管这通常不是问题。
调用者必须分别使用low_bounds和high_bounds向量来指定参数的下限和上限。
长度必须为nvars +1 的params向量返回最佳参数。这个数组中的最后一项是这个最优参数集的标准函数值。
如果print_progress输入非零，频繁的进度报告将被打印到控制台屏幕上。
只分配了三个工作数组:一个用于保存“当前”群体，一个用于保存正在创建的群体，一个短数组用于跟踪最佳个体。我们使用failures来计算初始群体中随机产生的个体被拒绝的次数，通常是因为交易系统的交易太少。我们很快就会看到，它将被用来降低最低交易要求。为了安全起见，n_evals统计我们评估随机生成的个体以创建初始群体的总次数。这样可以紧急逃生，避免挂电脑。第一个popsize个体填充pop1数组，过度初始化进入pop2[0]。
dim = nvars + 1 ; // Each individual is nvars variables plus criterion
pop1 = (double *) malloc ( dim * popsize * sizeof(double)) ;
pop2 = (double *) malloc ( dim * popsize * sizeof(double)) ;
best = (double *) malloc ( dim * sizeof(double)) ;
failures = 0 ; // Counts consecutive failures
n_evals = 0 ; // Counts evaluations for catastrophe escape
for (ind=0 ; ind<popsize+overinit ; ind++) {
if (ind < popsize) // If we are in pop1
popptr = pop1 + ind * dim ; // Point to the slot in pop1
else // Use first slot in pop2 for work
popptr = pop2 ; // Point to first slot in pop2
我们现在生成一个随机个体，并将其放入popptr。第一个nints参数为整数，其余为实数。这两种类型都是通过在各自的指定范围内统一选择值来生成的。然而，整数和实数的处理略有不同。函数unifrand()生成一个 0-1 范围内的均匀随机数。
for (i=0 ; i<nvars ; i++) { // For all variables (parameters)
if (i < nints) { // Is this an integer?
popptr[i] = low_bounds[i]+(int)(unifrand() * (high_bounds[i]-low_bounds[i] + 1.0));
if (popptr[i] > high_bounds[i]) // Virtually impossible, but be safe
popptr[i] = high_bounds[i] ;
}
else // real
popptr[i] = low_bounds[i] + (unifrand () * (high_bounds[i] - low_bounds[i])) ;
} // For all parameters
评估此人的交易系统表现，参数设置在popptr中。将此性能保存在popptr的最后一个槽中，紧接在参数之后。回想一下，每个槽都是nvars+1长。在构建初始群体时，计算性能评估的数量，以便我们可以将其作为紧急出口，以避免陷入表面上(或实际上！)死循环。最后，初始化第一个被测试者的最佳、最差和平均表现。这个memcpy()将这个人的参数和表现复制到短数组中，在这里我们记录下有史以来最好的成绩，最终返回给用户。
value = criter ( popptr , mintrades ) ;
popptr[nvars] = value ; // Also save criterion after variables
++n_evals ; // Count evaluations for emergency escape
if (ind == 0) {
grand_best = worstf = avgf = value ;
memcpy ( best , pop1 , dim * sizeof(double) ) ; // Best so far is first!
}
下一个代码块处理被拒绝的个体。请注意，此代码使用零阈值来拒绝参数集，例如显示亏损或未能满足最小交易计数要求。如果您想要使用一个不同的性能标准，一个不适合这个阈值的标准，您应该修改这个代码，或者更好的是，转换您的性能标准。例如，如果您想要最大化利润因子，那么合适的阈值应该是 1 而不是 0，您可以将您的绩效定义为利润因子的对数。
在下面显示的拒绝处理代码中，我们首先检查我们是否有如此糟糕的交易系统，以至于生成初始群体所需的评估数量已经失控，在这种情况下，我们采取紧急退出。如果没有，我们计算这种失败的次数。如果达到了一个很大的数字(500 在这里是硬编码的；随意更改)，我们重置失败计数器并降低最低交易要求，因为根据我的经验，这是最常见的失败原因，除非mintrades被设置得非常小。在任何情况下，这个个体的失败导致它被跳过，而成功重置失败计数器。因此，需要一个批次的失败来触发最小交易计数的减少。在采取这种激烈的行动之前，情况必须非常糟糕。
if (value <= 0.0) { // If this individual is totally worthless
if (n_evals > max_evals) // Safety escape should ideally never happen
goto FINISHED ;
--ind ; // Skip it entirely
if (++failures >= 500) { // This many in a row
failures = 0 ;
mintrades = mintrades * 9 / 10 ;
if (mintrades < 1)
mintrades = 1 ;
}
continue ;
}
else
failures = 0 ;
我们保持最佳、最差和一般的表现。后两个是严格的进度报告，如果用户不会更新进度，最差和平均计算可以省略。
if (value > grand_best) { // Best ever
memcpy ( best , popptr , dim * sizeof(double) ) ;
grand_best = value ;
}
if (value < worstf)
worstf = value ;
avgf += value ;
如果我们已经找到了popsize个个体，我们就进入了过度初始化。在现有人群中搜索最差的个体。如果这个新的过度初始化的个体优于最差的，用它替换最差的，这样就改善了基因库。回想一下，性能存储在参数之后，所以它在索引[nvars]处。同样，我们只对用户更新保持平均性能；它在优化算法中不起作用。
if (ind >= popsize) { // If we finished pop1, now doing overinit
avgf = 0.0 ;
minptr = NULL ; // Not needed. Shuts up 'use before define'
for (i=0 ; i<popsize ; i++) { // Search pop1 for the worst
dtemp = (pop1+i*dim)[nvars] ;
avgf += dtemp ;
if ((i == 0) || (dtemp < worstf)) {
minptr = pop1 + i * dim ;
worstf = dtemp ;
}
} // Searching pop1 for worst
if (value > worstf) { // If this is better than the worst, replace worst with it
memcpy ( minptr , popptr , dim * sizeof(double) ) ;
avgf += value - worstf ; // Account for the substitution
}
} // If doing overinit
} // For all individuals (population and overinit)
此时，我们已经完全生成了初始群体。找到最好的执行者，因为我们偶尔会在上面做一点爬山(除非用户禁止这样做，这通常是一个糟糕的举动)。然后将点设置为旧(源)代和新(目的)代，并将收敛计数器归零。我们将使用n_tweaked来控制爬山。
ibest = n_tweaked = 0 ;
value = pop1[nvars] ;
for (ind=1 ; ind<popsize ; ind++) {
popptr = pop1 + ind * dim ; // Point to the slot in pop1
if (popptr[nvars] > value) {
value = popptr[nvars] ;
ibest = ind ;
}
}
old_gen = pop1 ; // This is the old, parent generation
new_gen = pop2 ; // The children will be produced here
bad_generations = 0 ; // Counts contiguous generations with no improvement of best
我们有嵌套循环，代是外部循环，代中的个体是内部循环。仅针对可选用户更新跟踪平均和最差情况。变量improved标记最佳个体是否在这一代中的任何一点上有所改进。这用于表示收敛。主父代parent1来自源群体，我们将创建的子代将进入目标群体。
for (generation=1 ; ; generation++) {
worstf = 1.e60 ;
avgf = 0.0 ;
improved = 0 ; // Will flag if we improved in this generation
for (ind=0 ; ind<popsize ; ind++) { // Generate all children for this generation
parent1 = old_gen + ind * dim ; // Pure (and tested) parent
dest_ptr = new_gen + ind * dim ; // Winner goes here for next gen
我们随机选择第二个父代和两个差异。这些必须不同于主要的亲本并且彼此不同。
do { i = (int) (unifrand() * popsize) ; }
while ( i >= popsize || i == ind ) ;
do { j = (int) (unifrand() * popsize) ; }
while ( j >= popsize || j == ind || j == i ) ;
do { k = (int) (unifrand() * popsize) ; }
while ( k >= popsize || k == ind || k == i || k == j ) ;
parent2 = old_gen + i * dim ; // Parent to mutate
diff1 = old_gen + j * dim ; // First differential vector
diff2 = old_gen + k * dim ; // Second differential vector
下面的代码负责变异和交叉来创建一个新的子节点。我们将遍历每个参数，随机决定每个参数是否变异和交叉。如果我们到了最后一个参数，但还没有这样做，我们就这样做，以确保至少有一个变化。我们随机选择一个起始参数，这样当我们到达终点时，我们不会总是在同一个地方进行最后的动作。这种变异很容易将参数推到合法范围之外。根据需要修复此问题。ensure_legal()例程将在后面讨论。
do { j = (int) (unifrand() * nvars) ; }
while ( j >= nvars ) ; // Pick a starting parameter
used_mutated_parameter = 0 ; // We must act at least once; we haven’t yet
for (i=nvars-1 ; i>=0 ; i--) {
if ((i == 0 && ! used_mutated_parameter) || (unifrand() < pcross)) {
dest_ptr[j] = parent2[j] + mutate_dev * (diff1[j] - diff2[j]) ;
used_mutated_parameter = 1 ;
} // We mutated this variable
else // We did not mutate this variable, so copy old value
dest_ptr[j] = parent1[j] ;
j = (j + 1) % nvars ; // Rotate through all variables
}
ensure_legal ( nvars , nints , low_bounds , high_bounds , dest_ptr ) ;
评估这个新创建的孩子的表现。如果优于初级亲本，放入目的群体。否则，将主要父代复制到目标群体中。记录有史以来最好的个人，它最终会返回给呼叫者。通过improved表明我们这一代有所进步，所以我们还不准备放弃。变量n_tweaked将很快与爬山结合使用。
value = criter ( dest_ptr , mintrades ) ;
if (value > parent1[nvars]) { // If the child is better than parent1
dest_ptr[nvars] = value ; // Get the child's value (The vars are already there)
if (value > grand_best) { // And update best so far
grand_best = value ;
memcpy ( best , dest_ptr , dim * sizeof(double) ) ;
ibest = ind ;
n_tweaked = 0 ;
improved = 1 ; // Flag that the best improved in this generation
}
}
else { // Else copy parent1 and its value
memcpy ( dest_ptr , parent1 , dim * sizeof(double) ) ;
value = parent1[nvars] ;
}
我们现在开始可选的(但非常有用的)爬山步骤。下面的代码是决定是否爬山和爬什么的逻辑。我们将在下一页讨论它。
if (pclimb > 0.0 &&
((ind == ibest && n_tweaked < nvars) || (unifrand() < pclimb))) {
if (ind == ibest) { // Once each generation tweak the best
++n_tweaked ; // But quit if done all vars
k = generation % nvars ; // Cycle through all vars
}
else { // Randomly choose an individual
k = (int) (unifrand() * nvars) ; // Which var to optimize
if (k >= nvars) // Safety only
k = nvars - 1 ;
}
如果用户指定pclimb =0，那么爬山(这里称为调整)将永远不会完成。假设我们可以做到这一点，检查两个条件，其中任何一个都将允许对这个个体进行单次攀爬操作，这个个体可能是新创建的子体，也可能是主要父体的副本。如果个体是目前为止最好的，并且它的所有变量都没有被调整，我们就调整它。回想一下，每次 grand best 改变时,n_tweaked被重置为零。如果这是最好的，我们计算这种调整，并根据代选择变量。最优秀的个体在连续多代中保持不变是很常见的，这种参数选择会导致调整在参数之间轮换，从而避免重复。
如果第一次测试失败(要么这不是最好的个体，要么它的所有参数都已经被调整)，那么我们掷骰子，随机决定是否调整这个个体中随机选择的参数。
整数和实数参数的调整方式不同，前者更简单、更快。下面是整数代码的一半:
if (k < nints) { // Is this an integer?
ivar = ibase = (int) dest_ptr[k] ;
ilow = (int) low_bounds[k] ;
ihigh = (int) high_bounds [k] ;
success = 0 ;
while (++ivar <= ihigh) {
dest_ptr[k] = ivar ;
test_val = criter ( dest_ptr , mintrades ) ;
if (test_val > value) {
value = test_val ;
ibase = ivar ;
success = 1 ;
}
else {
dest_ptr[k] = ibase ;
break ;
}
}
我们在ibase中保存了这个参数的当前值，因此如果没有发现改进，我们可以恢复它。我们将在当前值附近改变ivar以寻求改进。对其整个合法范围进行全面的全球搜索通常是浪费时间。如果我们发现任何改进，变量success就会标记。我们向上移动参数，直到它达到上限或者性能没有提高。(为了保持快速搜索，这里忽略了平坦性能之后进行改进的可能性。)只要我们在改进，我们就不断更新ibase和改进后的值。当性能没有提高时，这可能发生在第一次测试中，我们将参数设置为ibase，并停止前进。
如果增加参数没有成功，我们会尝试减少参数。该算法本质上与向上搜索算法相同，因此在这里讨论它没有意义。
if (! success) {
ivar = ibase ;
while (--ivar >= ilow) {
dest_ptr[k] = ivar ;
test_val = criter ( dest_ptr , mintrades ) ;
if (test_val > value) {
value = test_val ;
ibase = ivar ;
success = 1 ;
}
else {
dest_ptr[k] = ibase ;
break ;
}
} // While searching downward
} // If the upward search did not succeed
} // If k < nints (this parameter is an integer)
处理实参数的代码稍微复杂一点。如下一页所示，我们从将性能计算所需的信息复制到静态变量开始，所有这些都以local_开始。该技术允许参数优化例程是通用的，调用仅引用被优化的参数的标准函数。
else { // This is a real parameter
local_criter = criter ;
local_ivar = k ; // Pass it to criterion routine
local_base = dest_ptr[k] ; // Preserve orig var
local_x = dest_ptr ;
local_nvars = nvars ;
local_nints = nints ;
local_low_bounds = low_bounds ;
local_high_bounds = high_bounds ;
old_value = value ;
优化分两步完成。首先，我们调用一个粗略的全局搜索例程glob_max()(GLOB _ MAX 中的源代码。CPP)，它测试一个范围内的一些离散点，并找到具有最大函数值的点。如果该值在一个端点处增加，它会一直前进，直到找到一个峰值。然后在brentmax()(Brent max 中的源代码)中使用 Brent 的算法细化这个最大值。CPP)。不幸的是，这可能是一个昂贵的操作。但是回报往往是巨大的，特别是当差异进化已经让我们接近全局最大值，而我们所需要的只是最佳个体的精确最大化。
我们在参数的当前值附近开始粗略的全局搜索:
lower = local_base - 0.1 * (high_bounds[k] - low_bounds[k]) ;
upper = local_base + 0.1 * (high_bounds[k] - low_bounds[k]) ;
if (lower < low_bounds[k]) {
lower = low_bounds[k] ;
upper = low_bounds[k] + 0.2 * (high_bounds[k] - low_bounds[k]) ;
}
if (upper > high_bounds[k]) {
upper = high_bounds[k] ;
lower = high_bounds[k] - 0.2 * (high_bounds[k] - low_bounds[k]) ;
}
k = glob_max ( lower , upper , 7 , 0 , c_func ,
&x1 , &y1 , &x2 , &y2 , &x3 , &y3 ) ;
此时，我们有三个点，使得中心点具有最大函数值。对此进行细化，并调用ensure_legal()来确保参数在其合法范围内。很可能是这种情况，或者至少非常接近，因为当超过合法界限时，标准函数应用巨大的惩罚，并且最大化例程将对该惩罚做出有力的响应。如果性能已经得到改善，即使在强制合法性之后(这几乎总是会发生的情况)，保存高级参数并更新 grand best。最后，更新最差和平均性能，严格用于用户更新(不是算法的一部分)。
brentmax ( 5 , 1.e-8 , 0.0001 , c_func , &x1 , &x2 , &x3 , y2 ) ;
dest_ptr[local_ivar] = x2 ; // Optimized var value
ensure_legal ( nvars , nints , low_bounds , high_bounds , dest_ptr ) ;
value = criter ( dest_ptr , mintrades ) ;
if (value > old_value) {
dest_ptr[nvars] = value ;
}
else {
dest_ptr[local_ivar] = local_base ; // Restore original value
value = old_value ;
}
if (value > grand_best) { // Update best so far
grand_best = value ;
memcpy ( best , dest_ptr , dim * sizeof(double) ) ;
ibest = ind ;
n_tweaked = 0 ;
improved = 1 ; // Flag that the best improved in this generation
}
} // If optimizing real parameter
} // If doing hill-climbing step
if (value < worstf)
worstf = value ;
avgf += value ;
} // Create all children in this generation
我们差不多完成了。如果这一代发现最佳个体没有改进，则增加收敛计数器，如果达到用户指定的计数，则退出。但是如果我们真的得到了改善，重置计数器。然后针对源和目的地生成群体颠倒pop1和pop2的角色。剩下的少量代码只是清理工作，这里省略了。
if (! improved) {
++bad_generations ;
if (bad_generations > max_bad_gen)
goto FINISHED ;
}
else
bad_generations = 0 ;
if (old_gen == pop1) {
old_gen = pop2 ;
new_gen = pop1 ;
}
else {
old_gen = pop1 ;
new_gen = pop2 ;
}
} // For all generations
确保合法性的例程只是对照用户指定的限制检查每个参数，计算超出限制的严格惩罚(仅用于实际参数调整)，并强制执行限制。对于整数，我们分别对待正值和负值，以确保正确的截断。回想一下，突变通常会导致整数参数获得非整数值，所以我们在这里首先解决这个问题。
static double ensure_legal ( int nvars , int nints , double *low_bounds , double
*high_bounds , double *params )
{
int i, j, varnum, ilow, ihigh ;
double rlow, rhigh, penalty, dtemp ;
penalty = 0.0 ;
for (i=0 ; i<nvars ; i++) {
if (i < nints) { // Is this an integer parameter?
if (params[i] >= 0)
params[i] = (int) (params[i] + 0.5) ;
else if (params[i] < 0)
params[i] = -(int) (0.5 - params[i]) ;
}
if (params[i] > high_bounds[i]) {
penalty += 1.e10 * (params[i] - high_bounds[i]) ;
params[i] = high_bounds[i] ;
}
if (params[i] < low_bounds[i]) {
penalty += 1.e10 * (low_bounds[i] - params[i]) ;
params[i] = low_bounds[i] ;
}
}
return penalty ;
}
由glob_max()和brentmax()调用的程序是被优化的单个参数的简单函数。将适当的参数设置为试验值，并调用ensure_legal()来执行合法性并计算超出界限的可能惩罚。然后调用性能计算例程来计算交易系统的性能，并且从性能中减去惩罚(如果有的话)。
static double c_func ( double param )
{
double penalty ;
local_x[local_ivar] = param ;
penalty = ensure_legal ( local_nvars , local_nints , local_low_bounds ,
local_high_bounds , local_x ) ;
return local_criter ( local_x , mintrades ) - penalty ;
}
一个完整的程序来演示这个算法。CPP 将在第 112 页介绍。
扫一扫
576
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
