如何在 OpenTURNS 中使用 BIC 对分布列表进行排名？

Question

我在 OpenTURNS 中有一个 Sample，我想在它上面安装一个发行版。为了考虑参数的数量，我想使用贝叶斯信息标准（BIC）。贝叶斯信息准则 (BIC) 根据加权最大似然准则对模型列表进行排名，该准则考虑了样本大小和每个分布的参数数量。 BIC 分数越低越好。

我知道FittingTest.BestModelBIC 返回最适合数据的模型。但是，我希望看到的不仅仅是最合适的：也许第二好的 BIC 对我来说有更多的物理意义？

如何在 OpenTURNS 中执行此操作？

PS 这是BestModelBIC 的示例。

Answer 1

我建议使用GetContinuousUniVariateFactories 来获取连续分布工厂的列表，并将其与计算分布的BIC 分数的FittingTest.BIC 函数结合起来。最后，可以使用Sample的sortAccordingToAComponent方法对表格进行排序。

这里是这个方法的详细例子。我们首先生成一个样本。

import openturns as ot
import tqdm
import pandas as pd
sample = ot.Normal().getSample(100)

GetContinuousUniVariateFactories 静态方法返回用于连续分布的所有可用工厂的列表。我们可以使用这个列表而无需进一步处理，但直方图会在排名中排在第一位，因为它是专门为此目的而设计的。因此，我们不将其包含在 BIC 分数的计算中。

factories = ot.DistributionFactory.GetContinuousUniVariateFactories()
marginalFactories = []
for factory in factories:
    if str(factory) != "HistogramFactory":
        print(factory)
        marginalFactories.append(factory)
number_of_factories = len(marginalFactories)
print("Number of selected factories:", number_of_factories)

打印出来：

ArcsineFactory
BetaFactory
BurrFactory
[...]
Number of selected factories: 29

在以下脚本中，我们对之前创建的列表中的所有工厂执行for 循环。稍后我们将按升序对 BIC 分数进行排序。这就是我们将 BIC 分数和边际指数存储在 score_array 样本中的原因。对于某些分布，计算可能会很长。因此我们使用tqdm 模块来打印进度条。最后，某些发行版并非基于此特定样本。为了避免破坏 for 循环，我们将对 BIC 方法的调用包装到 try/except 中。如果分布拟合失败，我们将 BIC 分数设置为浮点数的最大有限值（即MaxScalar），大约等于 $10^{308}$。

score_array = ot.Sample(number_of_factories, 2)
for i in tqdm.tqdm(range(number_of_factories)):
    try:
        factory = marginalFactories[i]
        fitted_dist, bic_score = ot.FittingTest.BIC(sample, factory)
        score_array[i] = [i, bic_score]
    except TypeError:
        print("Cannot build ", factory)
        score_array[i] = [i, ot.SpecFunc.MaxScalar]

关键步骤是对包含 BIC 分数的数组进行排序。

sorted_BIC_scores = score_array.sortAccordingToAComponent(1)

可能有超过 30 个分布可以构建到样本上。在这里，我们将列表限制为 BIC 分数最低的前 10 个分布。我们将使用 Pandas 来很好地打印 BIC 分数。为此，我们创建了BIC_data 列表，其中包含工厂名称和相应的 BIC 分数。这是使用sorted_BIC_scores 第一列中的分布索引的地方。但是，Sample 存储 floats：我们必须将它们转换为整数，然后才能将其用作索引。

BIC_data = []
rank = list(range(min(number_of_factories, 10)))
for i in rank:
    distribution_index = int(sorted_BIC_scores[i, 0])
    factory = marginalFactories[distribution_index]
    BIC_score = sorted_BIC_scores[i, 1]
    BIC_data.append([factory, BIC_score])

现在是最简单的部分，我们最终使用 Pandas 的DataFrame。

columns = ["BIC", "Factory"]
df = pd.DataFrame(BIC_data, index=rank, columns=columns,)
print(df)

打印出来的内容如下：

Rank	BIC	Factory
0	NormalFactory	2.902476
1	WeibullMaxFactory	2.933988
2	WeibullMinFactory	2.938930
3	TruncatedNormalFactory	2.939140
4	LogisticFactory	2.945102
5	LogNormalFactory	2.948479
6	StudentFactory	2.948733
7	TriangularFactory	2.968305
8	BetaFactory	3.033244
9	RayleighFactory	3.051117

我们看到，幸运的是，NormalFactory 适合高斯样本。截断正态工厂排在第 4 位，在两种 Weibull 分布之后。