场函数元模型 OpenTurns 1.16rc1

Meta-model of the field function OpenTurns 1.16rc1

将 Openturns 从 1.15 更新到 1.16rc1 后,我在构建字段函数的元模型时遇到以下问题:

减少计算负担:

ot.ResourceMap.SetAsUnsignedInteger("FittingTest-KolmogorovSamplingSize", 1)
algo = ot.FunctionalChaosAlgorithm(sample_X, outputSampleChaos)
algo.run()
metaModel = ot.PointToFieldConnection(postProcessing, algo.getResult().getMetaModel())

“FittingTest-KolmogorovSamplingSize”已从 OpenTurns 1.16rc1 中删除,当我尝试将拟合测试替换为:

ot.ResourceMap.SetAsUnsignedInteger("FittingTest-LillieforsMaximumSamplingSize", 10)

或与

ot.ResourceMap.SetAsUnsignedInteger("FittingTest-LillieforsMinimumSamplingSize", 1)

代码冻结。有解决办法吗?

解决方案是使用:

degree = 1
dimension_xi_X = 3
dimension_xi_Y = 450
enumerateFunction = ot.LinearEnumerateFunction(dimension_xi_X)
basis = ot.OrthogonalProductPolynomialFactory(
    [ot.HermiteFactory()] * dimension_xi_X, enumerateFunction)
basisSize =450 #enumerateFunction.getStrataCumulatedCardinal(degree)
adaptive = ot.FixedStrategy(basis, basisSize)
projection = ot.LeastSquaresStrategy(
    ot.LeastSquaresMetaModelSelectionFactory(ot.LARS(), ot.CorrectedLeaveOneOut()))
ot.ResourceMap.SetAsScalar("LeastSquaresMetaModelSelection-ErrorThreshold", 1.0e-7)
algo_chaos = ot.FunctionalChaosAlgorithm(sample_X, 
outputSampleChaos,basis.getMeasure(), adaptive, projection)
algo_chaos.run()
result_chaos = algo_chaos.getResult()
meta_model = result_chaos.getMetaModel()
metaModel1 = ot.PointToFieldConnection(postProcessing, 
algo_chaos.getResult().getMetaModel())

建议的解决方案是简单地使用另一个分布来为您的数据建模。您可以使用适当维度的任何其他多元连续分布。 IMO 这不是一个有效的答案,因为分布没有 link 到您的数据。

经检查,问题似乎与Lilliefors的测试无关。在 OT 1.15 中,我们使用此测试(以 Kolmogorov 的错误名称)自动 select 适合输入样本的分布,但我们切换到更复杂的 selection 算法(参见 MetaModelAlgorithm::BuildDistribution).它基于使用原始 Kolomgorov 测试的第一遍(因此忽略了参数已被估计的事实),然后 information-based 标准用于 select 最相关的模型(AIC、AICC、BIC 取决于关于 ResourceMap 中“MetaModelAlgorithm-ModelSelectionCriterion”键的值。问题是由 TrapezoidalFactory class 在 Kolmogorov 阶段引起的。我将尽快在 OpenTURNS master 中提供修复。与此同时,我有将建议的解决方案调整为更适合您的数据的内容:

degree = 6
dimension_xi_X = 3
dimension_xi_Y = 450
enumerateFunction = ot.HyperbolicAnisotropicEnumerateFunction(dimension_xi_X, 0.8)
basis = ot.OrthogonalProductPolynomialFactory(
    [ot.StandardDistributionPolynomialFactory(ot.HistogramFactory().build(sample_X[:,i])) for i in range(dimension_xi_X)], enumerateFunction)
basisSize = enumerateFunction.getStrataCumulatedCardinal(degree)
#basis = ot.OrthogonalProductPolynomialFactory(
#    [ot.HermiteFactory()] * dimension_xi_X, enumerateFunction)
#basisSize = 450#enumerateFunction.getStrataCumulatedCardinal(degree)
adaptive = ot.FixedStrategy(basis, basisSize)
projection = ot.LeastSquaresStrategy(
    ot.LeastSquaresMetaModelSelectionFactory(ot.LARS(), ot.CorrectedLeaveOneOut()))
ot.ResourceMap.SetAsScalar("LeastSquaresMetaModelSelection-ErrorThreshold", 1.0e-7)
algo_chaos = ot.FunctionalChaosAlgorithm(sample_X, 
outputSampleChaos,basis.getMeasure(), adaptive, projection)
algo_chaos.run()
result_chaos = algo_chaos.getResult()
meta_model = result_chaos.getMetaModel()
metaModel = ot.PointToFieldConnection(postProcessing, 
algo_chaos.getResult().getMetaModel())

我还实现了 L2 错误的快速粗略估计器:

# Meta_model validation
iMax = 5
# Input values
sample_X_validation = ot.Sample(np.array(month_1_parameters_MSE.iloc[:iMax,0:3]))
print("sample size=", sample_X_validation.getSize())
# sample_X = ot.Sample(month_1_parameters_MSE[['Rseries','Rsh','Isc']])

# output values
#month_1_simulated.iloc[0:1].transpose()
Field = ot.Field(mesh,np.array(month_1_simulated.iloc[0:1]).transpose())
sample_Y_validation = ot.ProcessSample(1,Field)
for k in range(1,iMax):
    sample_Y_validation.add( np.array(month_1_simulated.iloc[k:k+1]).transpose() )


# In[18]:


graph = sample_Y_validation.drawMarginal(0)
graph.setColors(['red'])
drawables = graph.getDrawables()
graph2 = metaModel(sample_X_validation).drawMarginal(0)
graph2.setColors(['blue'])
drawables = graph2.getDrawables()
graph.add(graph2)
graph.setTitle('Model/Metamodel Validation')
graph.setXTitle(r'$t$')
graph.setYTitle(r'$z$')
drawables = graph.getDrawables()
L2_error = 0.0
for i in range(iMax):
    L2_error = (drawables[i].getData()[:,1]-drawables[iMax+i].getData()[:,1]).computeRawMoment(2)[0]
print("L2_error=", L2_error)

之前的答案错误为 79.488,新提案错误为 1.3994。这是图形比较。

Comparison between test data & previous answer Comparison between test data & new proposal