Re: [Help-gsl] Non-linear fitting and complex data

[Francesco Abbate]

2010/1/28 Saurav Pathak <address@hidden>

> Hi,
>
> I am trying to fit real parameters using a complex function.  That is, the
> GSL example:
> (
> http://www.gnu.org/software/gsl/manual/html_node/Example-programs-for-Nonlinear-Least_002dSquares-Fitting.html
> )
>
> int expb_f (const gsl_vector * x, void *data,            gsl_vector * f)
>
>
> f in my case is complex.  This means, my Jacobian too is complex.
> It is not clear to me from the example how I may specify this.   Could
> someone please help?
>

Hi,

GSL does not have direct support for non-linear fit of complex data. You can
still fit your data by working on the real and imaginary part. If the
fitting parameters are real numbers this is quite is to arrange, tipically
you have to pack the data in a gsl_vector by alternating the real and
imaginary part. The Jacobian will be splitted with the same logic. I can
give you more details if it can helps.

Otherwise, if you want to fit complex data of complex parameters the affair
will get more complicated because some components of the jacobian will be
related by the cauchy-riemann equations and you need to take ths into
account. I strongly advise in this case you use only real parameters because
otherwise it gets too much complicated.

If you are interested my program, gsl shell :

http://www.nongnu.org/gsl-shell/index.html

does have direct support for non linear fit of complex data with real
parameters. In the download directory you can find the latest release
gsl-shell-0.9.6-2 with the windows binary. Here a working example for fit of
complex data:

------------------
n = 50
p0 = vector {0.55, -4, 16, 0.23}
function fmodel(x, t)
   return x[1] * exp((x[2] + 1i*x[3]) * t + 1i * x[4])
end

r = rng() -- random number generator
y = cnew(n, 1, |k,j| fmodel(p0, (k-1)/n) + 0.01 * rnd.gaussian(r))

function cexpf(x, f, J)
   for k=1, n do
      local t, y = (k-1)/n, y[k]
      if f then f:set(k, 1, fmodel(x,t) - y) end
      if J then
 local e = exp((x[2] + 1i*x[3]) * t + 1i * x[4])
 J:set(k, 1, e)
 J:set(k, 2, t * x[1] * e)
 J:set(k, 3, 1i * t * x[1] * e)
 J:set(k, 4, 1i * x[1] * e)
      end
   end
end

function print_state(s)
      print ("x: ", s.x:row_print())
      print ("chi square: ", cmul(h(s.f), s.f)[1])
end

s = csolver {fdf= cexpf, n= n, p= 4, x0= vector {2.1, -2.8, 18, 0}}
repeat
      print_state (s)
      local status = s:iterate()
until status ~= 'continue'

print_state (s)
------------------------

You can also plot the results with the following commands:
------------------
pl = plot()
lr = path(0, real(y[1])) -- the data at t=0
lf = path(0, real(fmodel(p0, 0))) -- the model at t = 0
for k=1, y:dims() do
   lr:line_to((k-1)/n, real(y[k]))
   lf:line_to((k-1)/n, real(fmodel(p0, (k-1)/n)))
end

pl:add(lr, 'black', {{'stroke'}, {'marker'}})
pl:add_line(lf)
pl:show()
------------------
Best regards,
Francesco