[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: [Help-gsl] nonlinear fit
From: |
Francesco Abbate |
Subject: |
Re: [Help-gsl] nonlinear fit |
Date: |
Fri, 11 Mar 2011 10:44:23 +0100 |
2011/3/10 I.I. Stefan <address@hidden>:
> Hey,
>
> Thanks. I implemented that, but now the routine is not converging. This is
> the error I'm getting:
>
> status = the iteration has not converged yet
Hi,
it seems that the Jacobian is wrong. The sign should be changed in the line:
gsl_matrix_set (J, i, 1, (-(A/(sigma*sigma) * e * (t-mu))));
When this change is done it should work.
I include here a script for GSL Shell 1.1 to do the work. It will
generate the fake data, perform the linear fit and show a nice plot.
It can be useful for debugging a complex jacobian, you can easily do a
small error that spoil all the calculations.
Best regards,
Francesco
-------------------- 8< -------------------
-- work with gsl shell 1.1
local nb = 15
local nh = (nb-1)/2
local r = rng()
local sig, mu = 1.2, -2.3
local dx = 5*sig/nh
local xsmp= |i| mu + (i-1-nh)*dx
local data = {}
data.t = new(nb, 1, |i| xsmp(i))
data.y = new(nb, 1, |i| pdf.gaussian(xsmp(i) - mu, sig))
local function gaussfdf(x, f, J)
local sig, mu = x[1], x[2]
local A = 1/sqrt(2*pi*sig^2)
for k=1, nb do
local t, y = data.t[k], data.y[k]
local e = exp( -(t-mu)^2/(2*sig^2))
if f then f[k] = A * e - y end
if J then
J:set(k, 1, A/sig * e * ((t-mu)^2/sig^2 - 1))
J:set(k, 2, A * e * (t-mu)/sig^2)
end
end
end
s = cnlfsolver {fdf= gaussfdf, n= nb, p0= vector {0.8, 0.0}}
repeat
local status = s:iterate()
until status ~= 'continue'
print('FOUND', tr(s.p))
local p = plot('Simulated Gaussian Distribution')
local function get_bars(t, y)
local i = 0
return function()
i = i + 1
if i <= nb+1 then
local ti = i <= nb and t[i] or t[i-1] + dx
local yi = i <= nb and y[i] or y[nb]
return ti, yi
end
end
end
p:add(ibars(get_bars(data.t, data.y)), rgba(0,0,0.8, 0.6))
p:addline(fxline(|x| pdf.gaussian(x - s.p[2], s.p[1]), mu-5*sig, mu+5*sig))
p:show()
return p
-------------------- 8< -------------------