[SciPy-User] How to optimize 2D data with leastsq in Python?

[Matt Newville]

Nxkr,

On Sat, Jul 2, 2016 at 4:02 AM, nxkryptor nxkr <nxkryptor at gmail.com> wrote:

> I have data
> <https://drive.google.com/open?id=0BwwhEMUIYGyTQXdVR2dvLXhfZnc> of the
> form given below:
>
> X   Y1  Y2  Y3  Y4  Y5  Y61.42857 4.83    4.58    4.43    4.31    4.22    4.141.40845 3.87    3.63    3.49    3.38    3.3 3.231.38889 3.17    2.93    2.79    2.69    2.62    2.561.36986 2.65    2.41    2.27    2.17    2.11    2.051.35135 2.27    2.02    1.88    1.78    1.72    1.67
>     :     :       :       :       :        :      :1.01010 2.26    1.6 1.21    0.97    0.79    0.661.00000 2.01    1.44    1.1 0.88    0.73    0.62
>
> I am trying to optimize a function of the form:
>
> [image: enter image description here] <http://i.stack.imgur.com/xfX0Y.png> [image:
> enter image description here] <http://i.stack.imgur.com/aokSR.png>
>
> Now my y variable has 6 columns, I can optimize one column (in MWE #3)
> which I have done in the MWE. How can I optimize the function to get the
> parameters A, n and B for all the 6 column values? In other words how can I
> get the value one value of A, n and B for all 6 values of y.
>
> *MWE*
>
> import numpy as npfrom scipy.signal import argrelextremaimport scipy.interpolatefrom scipy.optimize import leastsq
>
> with open('./exp_fit.dat', "r") as data:
>     while True:
>         line = data.readline()
>         if not line.startswith('#'):
>             break
>     data_header = [i for i in line.strip().split('\t') if i]
>     _data_ = np.genfromtxt(data, names = data_header, dtype = None, delimiter = '\t')
> _data_.dtype.names = [j.replace('_', ' ') for j in _data_.dtype.names]
> data = np.array(_data_.tolist())
>
> x = _data_['X']
> x_interp = np.linspace(x[-1], x[0], 100)
>
> y = np.zeros(shape = (len(x), 6))
> y_interp = np.zeros(shape = (len(x_interp), 6))for i in range(6):
>     y[:, i] = data[:, i + 1]
>     tck = scipy.interpolate.splrep(x[::-1], y[::-1, i], s=0)
>     y_interp[::-1, i] = scipy.interpolate.splev(x_interp, tck, der = 0)
>
> def residuals(p, y, x):
>     A_lt, n_lt, B_lt, A_it, n_it, B_it, A_ht, n_ht, B_ht = p
>     err = y - (1 / (1 / ((A_it * (15 ** (-n_it)) * np.exp(B_it * x / 1000)) + \
>         (A_lt * (15 ** (-n_lt)) * np.exp(B_lt * x / 1000))) + \
>         1 / (A_ht * (15 ** (-n_ht)) * np.exp(B_ht * x / 1000))))
>     return err
> p1 = [2.789E-05, 1.47, 9.368E+03, 2.789E-05, 1.47, 9.368E+03, 2.789E-05, 1.47, 9.368E+03]
> plsq = leastsq(residuals, p1, args=(y_interp[:, 2], x_interp))
>
> The question is also available on SO <http://stackoverflow.com/questions/38155014/how-to-optimize-2d-data-with-leastsq-in-python>.
>
> nxkr
>
>
>
I strongly recommend writing a function for your power-law/exponential
function, and calling that 3 times.   I assume your code does what you
want, but it's far too messy to actually read.    Also, it seems somewhat
odd that you are interpolating your data onto a finer x-grid (but
backwards??) rather than simply fitting the data you actually have.   Maybe
there's a good reason for that.  Finally, I believe you can read in your
data much more simply with

     data = np.loadtxt('./exp_fit.dat', skiprows=3)

I'm afraid I do not actually understand what you are trying to do.  You ask
  "How can I optimize the function to get the parameters A, n and B for all
the 6 column values? In other words how can I get the value one value of A,
n and B for all 6 values of y".

Do you mean that you want to do 6 separate fits, one for each column?
Or perhaps you mean that you want to fit 9 parameter values (3 each A, n, B
for 3 of you exponential functions) to all of the data columns
simultaneously?   Or perhaps you mean something else entirely (for example
that you mean one value for A, but a different value for n and B for each
column).    Sorry, but I can't tell.

Finally, just so you're aware, fitting multiple damped exponential
functions is often very difficult.

--Matt Newville
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20160703/8a29fdc2/attachment.html>