Is there a way to fit a piecewise linear function to a list?

Question

I need to find a piecewise linear function that is a proper approximation for sunrise/sunset times of each day of year. Due to some limitations, I cannot use any trigonometric functions. So the sine-looking function has to be approximated by a piecewise linear one. Here is the code to generate data:

location = GeoPosition[{30, 60}];
year = DateRange["2022/3/21", "2023/3/21"];
sunrises = DateDifference[#, Sunrise[location, #], "Minute"] & /@ year;
sunsets = DateDifference[#, Sunset[location, #], "Minute"] & /@ year;

The sunrise data looks like this

I want to approximate this function by $8$ straight lines. Thus the proposed function is defined as

r[x_] := Piecewise[{
 {393 + ((-393 + b)*(-1 + x))/(-1 + x2), x <= x2},
 {b + ((-b + c)*(x - x2))/(-x2 + x3), x2 < x && x <= x3},
 {c + ((-c + d)*(x - x3))/(-x3 + x4), x3 < x && x <= x4},
 {d + ((-d + e)*(x - x4))/(-x4 + x5), x4 < x && x <= x5},
 {e + ((-e + f)*(x - x5))/(-x5 + x6), x5 < x && x <= x6},
 {f + ((-f + g)*(x - x6))/(-x6 + x7), x6 < x && x <= x7},
 {g + ((393 - g)*(x - x7))/(-x7 + x8), x7 < x && x <= x8}},
 393 + ((-393 + i)*(x - x8))/(366 - x8)]

Now the problem is, I cannot find such fit using any functions that I have tried. The closest thing that I have found was This question, whose proposed method failed with this error:

FindFit::nrlnum: The function value {0.,0.0204082,<<364>>} is not a list of real numbers with dimensions {366} at {b,c,d,e,f,g,h,x2,x3,x4,<<4>>} = {370.,350.,<<12>>}.

FWIW, I did set the initial values, as suggested there

init = {{b, 340}, {c, 320}, {d, 340}, {e, 390}, {f, 430}, {g, 450}, {h, 430},
 {x2, 50}, {x3, 80}, {x4, 120}, {x5, 200}, {x6, 270}, {x7, 300}, {x8, 330}};
NonlinearModelFit[First /@ sunrises, r[x], init, x]

Is there any way to do this?

Answer 1

You were on the right track.

I fitted the simular function data = Table[-70 Sin[x] + 400 // N, {x, 0, 2 Pi, 2 Pi/365}]; because didn't get sunrise data.

data = Table[-70 Sin[x] + 400 // N, {x, 0, 2 Pi, 2 Pi/365}]; 

ListPlot[data]

r[x_] = Piecewise[{{a + ((-a + b)*(-1 + x))/(-1 + x2), 
0 < x <= x2}, {b + ((-b + c)*(x - x2))/(-x2 + x3), 
x2 < x && x <= x3}, {c + ((-c + d)*(x - x3))/(-x3 + x4), 
x3 < x && x <= x4}, {d + ((-d + e)*(x - x4))/(-x4 + x5), 
x4 < x && x <= x5}, {e + ((-e + f)*(x - x5))/(-x5 + x6), 
x5 < x && x <= x6}, {f + ((-f + g)*(x - x6))/(-x6 + x7), 
x6 < x && x <= x7}, {g + ((h - g)*(x - x7))/(-x7 + x8), 
x7 < x && x <= x8}}, h + ((-h + i)*(x - x8))/(366 - x8)]

init = {{a, 380}, {b, 340}, {c, 320}, {d, 340}, {e, 390}, {f, 
   430}, {g, 450}, {h, 390}, {i, 390}, {x2, 50}, {x3, 80}, {x4, 
   120}, {x5, 200}, {x6, 270}, {x7, 300}, {x8, 330}}

nlfit = NonlinearModelFit[data, r[x], init, x]

nlfit["BestFitParameters"]

(*   {a -> 398.58, b -> 343.842, c -> 327.174, d -> 344.802, e -> 451.384, 
 f -> 470.609, g -> 467.262, h -> 445.643, i -> 400.682, 
 x2 -> 51.9572, x3 -> 91.9915, x4 -> 134.41, x5 -> 229.068, 
 x6 -> 264.087, x7 -> 295.086, x8 -> 326.01}   *)

Plot[nlfit[x], {x, 0, 366}, 
 Epilog -> {PointSize[.005], Point@Transpose[{Range[366], data}]}, 
 PlotStyle -> {Red, Thickness[.015]}]

It is interesting to observe, you do not get absoute symmetry here, although this data are symmetric. Seems Fit doesn't find the global minimum of leastsquares. Think, it depends on initial values you give.

Answer 2

Here is an example with 8 straight lines. It still needs some fine tuning, but this you can do yourself:

x[0] = 1; x[1] = 60; x[2] = 80; x[3] = 110; x[4] = 200; x[5] = 280; 
x[6] = 300; x[7] = 320;x[8] = 360;
fit[x_] = 
 Piecewise[
  Transpose[{Fit[
       Transpose[{Table[x, {x, x[#[[1]]], x[#[[2]]]}], 
         sunrises[[x[#[[1]]] ;; x[#[[2]]]]]}], {1, x}, x], 
      x[#[[1]]] <= x <= x[#[[2]]]}] & /@ Partition[Range[0, 8], 2, 1]]
Show[Plot[fit[x], {x, 1, 350}, PlotStyle -> Red], ListPlot[sunrises]]