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 acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Browse other questions tagged or ask your own question.