使用scipy-optimize求解带有变量和不确定性的线性方程？

Question

我想要最小化一组方程式，其中的变量是已知的，但带有不确定性。简单来说，我想要检验一下这些测量得到的变量是否符合方程式所给出的公式限制。这听起来像是我可以用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  

希望能得到一些正确方向的指引。

Answer 1

我觉得我遇到了一个很相似的问题。我对Python还比较陌生，主要是用它来用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

我的想法是把这些方程当作约束条件，然后像上面的例子那样，解决残差的总和。

如果不考虑不确定性，解决这个问题并不难。我想请教一下，如何在寻找最合适的参数时考虑这些不确定性。

Answer 2

根据给出的情况，这个问题是没有解决方案的。这是因为如果输入的 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")

之前的脚本会产生以下输出。

根据这个分布中 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

这将产生：

这表明，考虑到观察到的输入 x1 到 x4 的变化性，e3 可能为零。但 e1 和 e2 不能为零。对于这个样本的结论是，第三个方程大致由观察到的 x1 到 x4 的值解决，但前两个方程则没有。

Answer 3

这是我的方法。假设 x1-x4 大致呈正态分布，围绕每个均值（1个标准差的不确定性），那么这个问题就变成了要最小化误差平方和，并且有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])