首页 文章

Julia中的2D曲线拟合

提问于
浏览 次
2

我在朱莉娅有一个数组 Z ，代表2D高斯函数的图像 . 即 Z[i,j] 是像素i，j处的高斯的高度 . 我想确定高斯参数（均值和协方差），可能是通过某种曲线拟合 .

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

如何在朱莉娅中将高斯分布到2D阵列？

julia-lang curve-fitting data-fitting

2 回答

  • 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])
    回复于
  • 1

    原则上，您有一个损失功能

    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 ） . 另一种方法是将条件指定为优化器中的约束 .

    为了让Optimzer工作，你可能必须获得 loss_packed 的衍生物 . 要么必须找到手动计算它（通过 dist 的良好选择），或者通过使用日志转换更容易（如果你很幸运，你会找到一种方法将其减少为线性问题......） . 或者，您可以尝试找到一个自动区分的优化器 .

    回复于

相关问题