I am trying to fit a time-dependent curve at each time step. I do so in minimizing along $x_c$ the quadratic error between the curve and a reference solution $ 1/(1 + \exp\left(\sqrt{S}(x-x_c)\right) $:

$$ \frac{1}{N+2} \sum_{i=0}^{N+1} \left(y_i - \frac{1}{1 + \exp\left(\sqrt{S}(x_i-x_c)\right)}\right)^2 $$

Because the calculations are easy, I calculated the gradient so as to use a conjugate gradient method implemented in optimize.minimize in scipy.

So: fairly easy stuff.

Below is my MWE and here is my data file:

import numpy as np
from matplotlib import pyplot as plt
from scipy import optimize

def energie(centre, x, y, S):
    erreur = 0.
    for i in range(0, y.size):
        erreur += (y[i] - 1./(1.+np.exp(np.sqrt(S)*(x[i]-centre))))**2.
    return erreur/y.size

def gradenergie(centre, x, y, S):
    grad = 0.
    for i in range(0, y.size):
        grad += - (y[i] - 1./(1.+np.exp(np.sqrt(S)*(x[i]-centre)))) * np.sqrt(S) * \
        np.exp(np.sqrt(S)*(x[i]-centre)) * (1.+np.exp(np.sqrt(S)*(x[i]-centre)))**-2.
    return 2.*grad/y.size

T = np.arange(0., 36., 0.5/np.sqrt(2))
x = np.arange(301) * 0.5
q = np.load('data/data.npy')
centres = np.zeros(q.shape[0])

for n in range(q.shape[0]):
    # initial guess
    centreq = centres[n-1]
    # call to optimize.minimize
    resultatq = optimize.minimize(energie, centreq, method = 'CG', jac = gradenergie, args = (x, q[n , : ], 0.1))
    # write to centres
    centres[n] = resultatq.x

dcentresqdt = np.gradient(centres, T[1]-T[0])
plt.plot(dcentresqdt[2 : ])

which produces the following graph: enter image description here with obvious irregularities. Any idea why?



in the code centre or centreq stands for $x_c$, each q[n, :] stands for the $(y_i)$, and the $(x_i)$ are the same at each time step, they are in the array x.

Indeed I plot the derivative of centres with respect to time, but if I plot centres itself I have the same kind of problem.



I realized my problem is so standard that there is actually a scipy function that does the job in no time (101 seconds for Nelder Mead vs 0.39 for the version below):

from scipy.optimize import curve_fit

def kink(x, centre):
    S = 0.1
    return 1./(1.+np.exp(np.sqrt(S)*(x-centre)))

and later on in the code:

for n in range(q.shape[0]):
    centreq = centres[n-1]
    popt, pcov = curve_fit(kink, x, q[n , : ], p0 = centreq)
    centres[n] = popt
  • 1
    $\begingroup$ I find it hard to tell from this question what the relationship is between $y_i,x_i,x_c$ in the formula and centres,q in the code. What exactly is being plotted here? It looks like there are multiple sets of data $\{(y^{(k)}_{0:n-1},x^{(k)}_{0:n-1})\}_k$ each leading to an $x_c^{(k)}$ and then $x_c^{(k+1)}-x_c^{(k)}$ is plotted as a function of $k$? It probably makes sense to you, but could you add a mathematical description that more precisely matches the code, please? $\endgroup$ – Kirill Jul 13 '18 at 16:14
  • 1
    $\begingroup$ I just changed the method of minimization to Nelder-Mead and obtained a smooth curve. Although, I agree with @Kirill that it's not clear what the plot is about. $\endgroup$ – nicoguaro Jul 13 '18 at 17:39
  • $\begingroup$ @Kirill you are right, I'm adding details. $\endgroup$ – bela83 Jul 14 '18 at 21:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.