11

Can anyone provide an example of providing a Jacobian to a least squares function in scipy?

I can't figure out the method signature they want - they say it should be a function, yet it's very hard to figure out what input parameters in what order this function should accept.

nbro
  • 15,395
  • 32
  • 113
  • 196
George Karpenkov
  • 2,094
  • 1
  • 16
  • 36

1 Answers1

16

Here's the exponential decay fitting that I got to work with this:

import numpy as np
from scipy.optimize import leastsq

def f(var,xs):
    return var[0]*np.exp(-var[1]*xs)+var[2]

def func(var, xs, ys):
    return f(var,xs) - ys

def dfunc(var,xs,ys):
    v = np.exp(-var[1]*xs)
    return [v,-var[0]*xs*v,np.ones(len(xs))]

xs = np.linspace(0,4,50)
ys = f([2.5,1.3,0.5],xs)
yn = ys + 0.2*np.random.normal(size=len(xs))
fit = leastsq(func,[10,10,10],args=(xs,yn),Dfun=dfunc,col_deriv=1)

If I wanted to use col_deriv=0, I think that I would have to basically take the transpose of what I return with dfunc. You're quite right though: the documentation on this isn't so great.

Justin Peel
  • 46,722
  • 6
  • 58
  • 80
  • It does work, yet ironically still fails on the example from my previous question :P okay, i probably should pick some other method – George Karpenkov Oct 19 '10 at 06:43
  • Yes, as tillsten said. Basically, it's the difference between fitting for exponential decay and exponential growth. That's a big difference. I think that you'll have to try some other method that uses the second derivative to have a chance of solving when you guess the wrong sign to start. It might need an additional momentum term of something like that too. – Justin Peel Oct 20 '10 at 01:56
  • Thanks, Justin. 2 Years later, and still helping people. – Geoff Oct 24 '12 at 23:18