title: fitting a 5th order Bézier Curve to a data set
body: <p>This is a slightly specific problem, so a bit of knowledge of R and of <code>Bézier</code> curves is required to be of help... (thanks if you do!!)</p>
  
  <p>So I need some help with my R code: I have a series of discretely sampled observations and I am trying to fit a <code>Bézier</code> Curve of the 5th order through these points with simple <code>LSS</code> regression. I have some limitations on the position of the 6 control points:</p>
  
  <ul>
  <li>A &amp; B have the same Y-axis coordinate</li>
  <li>B &amp; C have the same X-axis coordinate</li>
  <li>C &amp; D have the same Y-axis coordinate</li>
  <li>D &amp; E have the same X-axis coordinate</li>
  <li>E &amp; F have the same Y-axis coordinate</li>
  <li>A is located on the observation 2 turning points ago from the last
  observation </li>
  <li>The X-axis coordinate of the last observation is
  somewhere between the X-axis coordinates of E and F</li>
  </ul>
  
  <p>Like this image:</p>
  
  <p><img src="https://i.stack.imgur.com/gyroG.jpg" alt="this"></p>
  
  <p>Say I have these data:</p>
  
  <pre><code>-0.01105 
  -0.01118 
  -0.01271 
  -0.01479
  -0.01729
  -0.01996
  -0.02250
  -0.02473
  -0.02554
  -0.02478
  -0.02207
  -0.01788
  -0.01319
  -0.00956
  </code></pre>
  
  <p>They have a "curvy" shape so a Bézier curve would fit: the result of my code is this image: the data are in red, the 5th order Bézier and its control points with their restrictions in blue:</p>
  
  <p>Like this image:</p>
  
  <p><img src="https://i.stack.imgur.com/4ilfT.jpg" alt="this"></p>
  
  <p>So you see that I have some kind of solution, but this is the problem:
  The X-axis location of right-most control point is always to the right of the last input data point, and to get an appropriate fit, I had to require a value of t (t goes from 0 to 1 in a Bézier) where t is at if the input data end (the "limit" variable in my code). How do I rewrite it so I don't have to do that anymore, and the horizontal spread of the t-values remains constant, also outside of the input data?
  (given the restrictions on the control points, and maximizing the fit of the part of the curve that overlaps with the input data)</p>
  
  <p>If you can help, please take a look at this R code, any help is .. much much appreciated and happy holidays!!
  ps: what I call exampledata.csv in my code is just the data above.</p>
  
  <p><div class="snippet" data-lang="js" data-hide="false" data-console="false" data-babel="false">
  <div class="snippet-code">
  <pre class="snippet-code-html lang-html prettyprint-override"><code>getT &lt;- function(x){
   # Calculates length from origin of each point in the path.
   # args:
   #  x : a one dimensional vector
   # Returns:
   #  out : a vector of distances from the origin, as a percent of end point - start point distance
   
   out &lt;- cumsum(abs(diff(x)))
   out &lt;- c(0, out/ out[length(out)])
   return(out)
  }
  
  
  cost_f &lt;- function(X,Y,K){
   pred &lt;-K%*%X
   c &lt;- Y- pred
   out &lt;- list(loss= as.vector(t(c)%*%c), pred = pred)
   return(out)
  }
  
  df &lt;- read.csv('exampledata.csv')
  T &lt;- nrow(df)
  df['d'] = 1:T
  
  
  # # identify all turning points:
  # turn_point &lt;- c(1)
  # for(i in 2:(T-1)){
   # if( ( (df[i,'x'] &lt; df[i-1,'x']) &amp; (df[i,'x'] &lt; df[i+1,'x'])) | ( (df[i,'x'] &gt; df[i-1,'x']) &amp; (df[i,'x'] &gt; df[i+1,'x']))   ){
    # turn_point &lt;- c(turn_point, i)
   # }
   
  # }
  
  
  
  fit_last_piece &lt;- function(df){
   limit &lt;- .79
   turn_point &lt;- c(1)
   for(i in 2:(T-1)){
    if( ( (df[i,'x'] &lt; df[i-1,'x']) &amp; (df[i,'x'] &lt; df[i+1,'x'])) | ( (df[i,'x'] &gt; df[i-1,'x']) &amp; (df[i,'x'] &gt; df[i+1,'x']))   ){
     turn_point &lt;- c(turn_point, i)
    }
    
   }
   nk &lt;- length(turn_point) # number of turning points
   data &lt;- df[turn_point[nk-1]:nrow(df),]
   end_x &lt;- data$d[1]
   end_y &lt;- data$x[1]
    
   constr_x &lt;- matrix(c(1,0,0,0,0,0,  # remember data is input column to column
         0,1,1,0,0,0,
         0,0,0,1,1,0,
         0,0,0,0,0,1),nrow = 6, ncol = 4)
   constr_y &lt;- matrix(c(1,1,0,0,0,0,
         0,0,1,1,0,0,
         0,0,0,0,1,1),nrow = 6, ncol = 3)
   
   M = matrix(c(-1,5,-10,10,-5,1,
       5,-20,30,-20,5,0,
       -10,30,-30,10,0,0,
       10,-20,10,0,0,0,
       -5,5,0,0,0,0,
       1,0,0,0,0,0),nrow = 6, ncol = 6)
   
   t_x = getT(data$d)*limit
   T_x = cbind(t_x^5, t_x^4 ,t_x^3, t_x^2, t_x,rep(1,length(t_x)))
   
   in_par &lt;- ( tail(data$d,1)-data$d[1])*c(2/5,4/5,6/5) + data$d[1]   # initial values of the intermediate x levels are at 1/3 and 2/3 midpoints 
   res_x &lt;- optim(par = in_par, fn = function(par){cost_f(c(data$d[1], par[1],par[2], par[3]), data$d, T_x%*%M%*%constr_x)$loss})
   
   
   #res_x &lt;- optimize(f = function(par){cost_f(c(df$d[1],par,df$d[nrow(df)]), df$d, T_x%*%M%*%constr_x)$loss}, interval = c(df$d[1],df$d[nrow(df)]),tol = .Machine$double.eps^0.25)
   optim_x &lt;- c(data$d[1],res_x$par)
   pred_x &lt;- cost_f(optim_x, data$d, T_x%*%M%*%constr_x)$pred
   
   
   t_y = getT(data$x)*limit
   T_y = cbind(t_y^5, t_y^4,t_y^3, t_y^2, t_y,rep(1,length(t_y)))
   in_par &lt;- c()
   res_y &lt;- optim(par = c(data$x[floor(nrow(data)/2)],tail(data$x,1)), fn = function(par){cost_f(c(data$x[1],par[1],par[2]), data$x, T_y%*%M%*%constr_y)$loss})
   optim_y &lt;- c(data$x[1],res_y$par[1],res_y$par[2])
   #pred_y &lt;- cost_f(res_y$par, df$x, T_y%*%M%*%constr_y)$pred
   pred_y &lt;- cost_f(optim_y, data$x, T_y%*%M%*%constr_y)$pred
  
   t_x_p &lt;- c(t_x,seq(tail(t_x,1),1,length.out = 10))
   T_x_p &lt;- cbind(t_x_p^5, t_x_p^4 ,t_x_p^3, t_x_p^2, t_x_p,rep(1,length(t_x_p)))
   t_y_p &lt;- c(t_y,seq(tail(t_y,1),1,length.out = 10))
   T_y_p &lt;- cbind(t_y_p^5, t_y_p^4 ,t_y_p^3, t_y_p^2, t_y_p,rep(1,length(t_y_p)))
   pred_x &lt;- T_x_p%*%M%*%constr_x%*%optim_x
   pred_y &lt;- T_y_p%*%M%*%constr_y%*%optim_y
  #  this part is new:
   plot(pred_x,pred_y, ylim = c(min(c(data$x, pred_y,res_y$par)), max(c(data$x, pred_y,res_y$par))),col="blue",type="b")
   points(data$d,data$x,col = 'red',type="b")
   points(pred_x[1],pred_y[1],pch=20,col='blue')
   points(res_x$par[1],pred_y[1],pch=20,col='blue')
   points(res_x$par[1],res_y$par[1],pch=20,col='blue')
   points(res_x$par[2],res_y$par[1],pch=20,col='blue')
   points(res_x$par[2],res_y$par[2],pch=20,col='blue')
   points(res_x$par[3],res_y$par[2],pch=20,col='blue')
   segments(pred_x[1],pred_y[1],res_x$par[1],pred_y[1],lty=3,col='blue')
   segments(res_x$par[1],pred_y[1],res_x$par[1],res_y$par[1],lty=3,col='blue')
   segments(res_x$par[1],res_y$par[1],res_x$par[2],res_y$par[1],lty=3,col='blue')
   segments(res_x$par[2],res_y$par[1],res_x$par[2],res_y$par[2],lty=3,col='blue')
   segments(res_x$par[2],res_y$par[2],res_x$par[3],res_y$par[2],lty=3,col='blue')
  }
  
  fit_last_piece(df)</code></pre>
  </div>
  </div>
  </p>