julia - Julia 中的二维曲线拟合

Question

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

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

How can I fit a Gaussian to a 2D array in 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])