help-gsl
[Top][All Lists]
Advanced

[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< -------------------



reply via email to

[Prev in Thread] Current Thread [Next in Thread]