title: Fitting piecewise function in Python
body: <p>I'm trying to fit a piecewise defined function to a data set in Python. I've searched for quite a while now, but I haven't found an answer whether it is possible or not.</p>
  
  <p>To get an impression of what I am trying to do, look at the following example (which is not working for me). Here I'm trying to fit a shifted absolute value function (f(x) = |x-p|) to a dataset with p as the fit parameter. </p>
  
  <pre class="lang-python prettyprint-override"><code>import scipy.optimize as so
  import numpy as np
  
  def fitfunc(x,p):
     if x&gt;p:
        return x-p
     else:
        return -(x-p)
  
  fitfunc = np.vectorize(fitfunc) #vectorize so you can use func with array
  
  x=np.arange(1,10)
  y=fitfunc(x,6)+0.1*np.random.randn(len(x))
  
  popt, pcov = so.curve_fit(fitfunc, x, y) #fitting routine that gives error
  </code></pre>
  
  <p>Is there any way of accomplishing this in Python? </p>
  
  <p>A way of doing this in R is :</p>
  
  <pre class="lang-R prettyprint-override"><code># Fit of a absolute value function f(x)=|x-p|
  
  f.lr &lt;- function(x,p) {
      ifelse(x&gt;p, x-p,-(x-p))
  }
  x &lt;- seq(0,10)  #
  y &lt;- f.lr(x,6) + rnorm (length(x),0,2)
  plot(y ~ x)
  fit.lr &lt;- nls(y ~ f.lr(x,p), start = list(p = 0), trace = T, control = list(warnOnly = T,minFactor = 1/2048))
  summary(fit.lr)
  coefficients(fit.lr)
  p.fit &lt;- coefficients(fit.lr)["p"]
  x_fine &lt;- seq(0,10,length.out=1000)
  lines(x_fine,f.lr(x_fine,p.fit),type='l',col='red')
  lines(x,f.lr(x,6),type='l',col='blue')
  </code></pre>
  
  <p>After even more research I found a way of doing it. In this solution, I don't like the fact that I have to define the error function myself. Further I'm not really sure why it has to be in this lambda-style. Therefore any kind of suggestions or more sophisticated solutions are very welcome.</p>
  
  <pre class="lang-python prettyprint-override"><code>import scipy.optimize as so
  import numpy as np
  import matplotlib.pyplot as plt
  
  def fitfunc(p,x): return x - p if x &gt; p else p - x 
  
  def array_fitfunc(p,x):
      y = np.zeros(x.shape)
      for i in range(len(y)):
          y[i]=fitfunc(x[i],p)
      return y
  
  errfunc = lambda p, x, y: array_fitfunc(p, x) - y # Distance to the target function
  
  x=np.arange(1,10)
  x_fine=np.arange(1,10,0.1)
  y=array_fitfunc(6,x)+1*np.random.randn(len(x)) #data with noise
  
  p1, success = so.leastsq(errfunc, -100, args=(x, y), epsfcn=1.) # -100 is the initial value for p; epsfcn sets the step width
  
  plt.plot(x,y,'o') # fit data
  plt.plot(x_fine,array_fitfunc(6,x_fine),'r-') #original function
  plt.plot(x_fine,array_fitfunc(p1[0],x_fine),'b-') #fitted version
  plt.show()
  </code></pre>