Julia中的二维曲线拟合

Question

我在 Julia 中有一个数组 Z，它表示二维高斯函数的图像。 IE。 Z[i,j] 是像素 i,j 处的高斯高度。我想通过某种曲线拟合来确定高斯的参数（均值和协方差）。

我已经研究了适合Z 的各种方法：我首先尝试了Distributions 包，但它是为稍微不同的情况而设计的（随机选择的点）。然后我尝试了LsqFit 包，但它似乎是为 1D 拟合量身定制的，因为当我尝试拟合 2D 数据时它会抛出错误，而且我找不到任何文档来引导我找到解决方案。

如何在 Julia 中将高斯拟合到二维数组？

Answer 1

最简单的方法是使用 Optim.jl。这是一个示例代码（它没有针对速度进行优化，但它应该向您展示如何处理该问题）：

using Distributions, Optim

# generate some sample data    
true_d = MvNormal([1.0, 0.0], [2.0  1.0; 1.0 3.0])
const xr = -3:0.1:3
const yr = -3:0.1:3
const s = 5.0
const m = [s * pdf(true_d, [x, y]) for x in xr, y in yr]

decode(x) = (mu=x[1:2], sig=[x[3] x[4]; x[4] x[5]], s=x[6])

function objective(x)
    mu, sig, s = decode(x)
    try # sig might be infeasible so we have to handle this case
        est_d = MvNormal(mu, sig)
        ref_m = [s * pdf(est_d, [x, y]) for x in xr, y in yr]
        sum((a-b)^2 for (a,b) in zip(ref_m, m))
    catch
        sum(m)
    end
end

# test for an example starting point
result = optimize(objective, [1.0, 0.0, 1.0, 0.0, 1.0, 1.0])
decode(result.minimizer)

或者，您可以使用约束优化，例如像这样：

using Distributions, JuMP, NLopt

true_d = MvNormal([1.0, 0.0], [2.0  1.0; 1.0 3.0])
const xr = -3:0.1:3
const yr = -3:0.1:3
const s = 5.0
const Z = [s * pdf(true_d, [x, y]) for x in xr, y in yr]

m = Model(solver=NLoptSolver(algorithm=:LD_MMA))

@variable(m, m1)
@variable(m, m2)
@variable(m, sig11 >= 0.001)
@variable(m, sig12)
@variable(m, sig22 >= 0.001)
@variable(m, sc >= 0.001)

function obj(m1, m2, sig11, sig12, sig22, sc)
    est_d = MvNormal([m1, m2], [sig11 sig12; sig12 sig22])
    ref_Z = [sc * pdf(est_d, [x, y]) for x in xr, y in yr]
    sum((a-b)^2 for (a,b) in zip(ref_Z, Z))
end

JuMP.register(m, :obj, 6, obj, autodiff=true)
@NLobjective(m, Min, obj(m1, m2, sig11, sig12, sig22, sc))
@NLconstraint(m, sig12*sig12 + 0.001 <= sig11*sig22)

setvalue(m1, 0.0)
setvalue(m2, 0.0)
setvalue(sig11, 1.0)
setvalue(sig12, 0.0)
setvalue(sig22, 1.0)
setvalue(sc, 1.0)

status = solve(m)
getvalue.([m1, m2, sig11, sig12, sig22, sc])

Answer 2

原则上，你有一个损失函数

loss(μ, Σ) = sum(dist(Z[i,j], N([x(i), y(j)], μ, Σ)) for i in Ri, j in Rj)

其中x 和y 将您的索引转换为轴上的点（您需要知道网格距离和偏移位置），Ri 和Rj 将索引的范围。 dist 是您使用的距离度量，例如。平方差。

您应该能够通过将 μ 和 Σ 打包到单个向量中来将其传递给优化器：

pack(μ, Σ) = [μ; vec(Σ)]
unpack(v) = @views v[1:N], reshape(v[N+1:end], N, N)
loss_packed(v) = loss(unpack(v)...)

在你的情况下N = 2。 （也许解包需要一些优化来摆脱不必要的复制。）

另一件事是我们必须确保Σ 是正半定数（因此也是对称的）。一种方法是对打包损失函数进行不同的参数化，并优化一些下三角矩阵L，例如Σ = L * L'。在N = 2的情况下，我们可以写成

unpack(v) = v[1:2], LowerTriangular([v[3] zero(v[3]); v[4] v[5]])
loss_packed(v) = let (μ, L) = unpack(v)
    loss(μ, L * L')
end

（这当然容易进一步优化，例如将乘法直接扩展为loss）。另一种方法是将条件指定为优化器的约束。

要使优化器工作，您可能必须得到loss_packed 的导数。要么必须找到手动计算它（通过dist 的一个很好的选择），或者通过使用对数转换更容易（如果你很幸运，你会找到一种方法将其减少为线性问题......） .或者，您可以尝试找到一个进行自动微分的优化器。