4

我想最小化一组方程，其中变量的不确定性是已知的。本质上，我想检验给定测量变量符合方程给出的公式约束的假设。这似乎是我应该能够用 scipy-optimize 做的事情。例如我有三个方程：

8 = 0.5 * x1 + 1.0 * x2 + 1.5 * x3 + 2.0 * x4  
4 = 0.0 * x1 + 0.0 * x2 + 1.0 * x3 + 1.0 * x4  
1 = 1.0 * x1 + 1.0 * x2 + 0.0 * x3 + 0.0 * x4

以及四个具有 1-sigma 不确定性的测量未知数： 

x1 = 0.246 ± 0.007  
x2 = 0.749 ± 0.010  
x3 = 1.738 ± 0.009  
x4 = 2.248 ± 0.007

寻找正确方向的任何指示。

4

3 回答 3

4

这是我的方法。假设x1-x4在每个均值（1-sigma 不确定性）周围近似正态分布，问题就变成了使用 3 个线性约束函数最小化误差平方和的问题。因此，我们可以使用攻击它scipy.optimize.fmin_slsqp()

In [19]:

def eq_f1(x):
    return (x*np.array([0.5, 1.0, 1.5, 2.0])).sum()-8
def eq_f2(x):
    return (x*np.array([0.0, 0.0, 1.0, 1.0])).sum()-4
def eq_f3(x):
    return (x*np.array([1.0, 1.0, 0.0, 0.0])).sum()-1
def error_f(x):
    error=(x-np.array([0.246, 0.749, 1.738, 2.248]))/np.array([0.007, 0.010, 0.009, 0.007])
    return (error*error).sum()
In [20]:

so.fmin_slsqp(error_f, np.array([0.246, 0.749, 1.738, 2.248]), eqcons=[eq_f1, eq_f2, eq_f3])
Optimization terminated successfully.    (Exit mode 0)
            Current function value: 2.17576389592
            Iterations: 4
            Function evaluations: 32
            Gradient evaluations: 4
Out[20]:
array([ 0.25056582,  0.74943418,  1.74943418,  2.25056582])
于 2014-04-22T23:01:07.810 回答
0

在我看来，我有一个非常相似的问题。我对 py 比较陌生，主要用它来排序和减少 pandas 的数据。

我有一组线性方程，我想在其中找到最佳拟合参数。但是，数据集具有已知的不确定性，需要在括号中考虑）。

x1*99(1)+x2*45(1)=52(0.2)
x1*1(0.5)+x2*16(1)=15(0.1)

此外还有一些限制：

x1>=0
x2>=0
x1+x2=1

我的方法是将方程视为约束并求解残差之和，如上面的示例所示。

在没有不确定性的情况下解决这个问题不是问题。我要求在找到最佳拟合参数的同时获得有关如何考虑不确定性的提示。

于 2020-04-30T12:07:33.803 回答
0

如给定的那样，问题没有解决方案。这是因为如果输入 x1、x2、x3 和 x4 是高斯的，那么输出：

y1 = 0.5 * x1 + 1.0 * x2 + 1.5 * x3 + 2.0 * x4 - 8.0
y2 = 0.0 * x1 + 0.0 * x2 + 1.0 * x3 + 1.0 * x4 - 4.0
y3 = 1.0 * x1 + 1.0 * x2 + 0.0 * x3 + 0.0 * x4 - 1.0

也是高斯的。假设 x1、x2、x3 和 x4 是独立的随机变量，这在OpenTURNS中很容易看出：

import openturns as ot
x1 = ot.Normal(0.246, 0.007)
x2 = ot.Normal(0.749, 0.010)
x3 = ot.Normal(1.738, 0.009)
x4 = ot.Normal(2.248, 0.007)
y1 = 0.5 * x1 + 1.0 * x2 + 1.5 * x3 + 2.0 * x4 - 8.0
y2 = 0.0 * x1 + 0.0 * x2 + 1.0 * x3 + 1.0 * x4 - 4.0
y3 = 1.0 * x1 + 1.0 * x2 + 0.0 * x3 + 0.0 * x4 - 1.0

以下脚本生成图表：

graph1 = y1.drawPDF()
graph1.setLegends(["y1"])
graph2 = y2.drawPDF()
graph2.setLegends(["y2"])
graph3 = y3.drawPDF()
graph3.setLegends(["y3"])
graph1.add(graph2)
graph1.add(graph3)
graph1.setColors(["dodgerblue3",
                   "darkorange1", 
                   "forestgreen"])
graph1.setXTitle("Y")

前面的脚本产生以下输出。

输出 Y 的分布

鉴于 0.0 在此分布中的位置，我会说求解方程在数学上是不可能的，但在物理上与数据一致。

实际上，我猜您为 x1, ..., x4 给出的高斯分布是根据数据估计的。所以我宁愿将问题重新表述如下：

给定 x1、x2、x3、x4 的观测值样本，e1、e2、e3 的值是多少，使得：

y1 = 0.5 * x1 + 1.0 * x2 + 1.5 * x3 + 2.0 * x4 - 8 + e1 = 0
y2 = 0.0 * x1 + 0.0 * x2 + 1.0 * x3 + 1.0 * x4 - 4 + e2 = 0
y3 = 1.0 * x1 + 1.0 * x2 + 0.0 * x3 + 0.0 * x4 - 1 + e3 = 0

这将问题变成了一个反演问题，可以通过校准e1，e2，e3来解决。此外，给定 x1, ..., x4 的有限样本大小，我们可能想要生成 e1, e2, e3 的分布。这可以通过引导输入/输出对 (x, y) 来完成：e1、e2、e3 的分布反映了这些参数的可变性，具体取决于手头的样本。

首先，我们必须从发行版中生成一个样本（我想你有这个样本，但到目前为止还没有发布）：

distribution = ot.ComposedDistribution([x1, x2, x3, x4])
sampleSize = 10
xobs = distribution.getSample(sampleSize)

然后我们定义模型：

formulas = [
    "y1 := 0.5 * x1 + 1.0 * x2 + 1.5 * x3 + 2.0 * x4 + e1 - 8.0",
    "y2 := 0.0 * x1 + 0.0 * x2 + 1.0 * x3 + 1.0 * x4 + e2 - 4.0",
    "y3 := 1.0 * x1 + 1.0 * x2 + 0.0 * x3 + 0.0 * x4 + e3 - 1.0"
]
program = ";".join(formulas)
g = ot.SymbolicFunction(["x1", "x2", "x3", "x4", "e1", "e2", "e3"],
                        ["y1", "y2", "y3"], 
                        program)

并设置观察到的输出，这是一个零样本：

yobs = ot.Sample(sampleSize, 3)

我们从初始值为零开始，并定义要校准的函数：

e1Initial = 0.0
e2Initial = 0.0
e3Initial = 0.0
thetaPrior = ot.Point([e1Initial,e2Initial,e3Initial])
calibratedIndices = [4, 5, 6]
mycf = ot.ParametricFunction(g, calibratedIndices, thetaPrior)

然后我们可以校准模型：

algo = ot.NonLinearLeastSquaresCalibration(mycf, xobs, yobs, thetaPrior)
algo.run()
calibrationResult = algo.getResult()
print(calibrationResult.getParameterMAP())

这打印：

[0.0265988,0.0153057,0.00495758]

这意味着误差 e1、e2、e3 相当小。我们可以计算一个置信区间：

thetaPosterior = calibrationResult.getParameterPosterior()
print(thetaPosterior.computeBilateralConfidenceIntervalWithMarginalProbability(0.95)[0])

这打印：

[0.0110046, 0.0404756]
[0.00921992, 0.0210059]
[-0.00601084, 0.0156665]

第三个参数 e3 可能为零，但既不是 e1 也不是 e2。最后，我们可以得到错误的分布：

thetaPosterior = calibrationResult.getParameterPosterior()

并绘制它：

graph1 = thetaPosterior.getMarginal(0).drawPDF()
graph2 = thetaPosterior.getMarginal(1).drawPDF()
graph3 = thetaPosterior.getMarginal(2).drawPDF()
graph1.add(graph2)
graph1.add(graph3)
graph1.setColors(["dodgerblue3",
                  "darkorange1", 
                  "forestgreen"])
graph1

这会产生：

e1、e2、e3的分布

这表明，考虑到观察到的输入 x1、...、x4 的可变性，e3 可能为零。但 e1 和 e2 不能为零。该样本的结论是第三个方程可以通过 x1、...、x4 的观测值近似求解，但不是第一个方程的两个。

于 2020-04-30T22:07:17.543 回答