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

Question

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

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

如何在OpenTURNS中执行此操作

PS 下面是BestModelBIC的一个示例

Answer 1

我建议使用GetContinuousUniVariateFactories来获取连续分发工厂的列表，并将其与FittingTest.BIC函数相结合，该函数计算分发的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个发行版。我们将使用熊猫来很好地打印BIC分数。为此，我们创建BIC_data列表，其中包含工厂名称和相应的BIC分数。这是使用sorted_BIC_scores第一列中分布的索引的地方。然而，Sample存储float：我们必须先将它们转换为整数，然后再将其用作索引

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)

这张照片看起来像这样：

^{tb1}$

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