I have acquired many sets of data, which all represent a single function, but which are randomly rescaled by a constant (due to the measurement specifics). I'm looking to effectively stitch them together as a continuous sort of function by rescaling each data set, however this has proven difficult since their ranges don't always overlap. Ideally something like:

Rescaling to continuous data

Where the resulting absolute scale doesn't matter, but the structural features are important.

The obvious solution is to interpolate/extrapolate nearby curves, and minimize the differences between neighbors. However, I haven't been able to make this work very well, as I'm not sure if there's a good way to select which curves should be paired/minimized together. Any suggestions?

Example={{{2.04,3.94},{2.46,3.81},{2.89,3.56},{3.1,3.18},{3.44,2.81},{3.75,2.42},{3.91,2.03},{4.12,1.75},{4.59,1.44},{5.,1.28},{5.14,1.17}},{{0.23,5.26},{0.4,6.02},{0.65,6.81},{0.96,7.47},{1.3,7.86},{1.68,7.96},{1.82,8.08},{2.15,7.84},{2.47,7.39},{2.78,6.78},{3.1,6.11},{3.43,5.33},{3.86,4.61},{4.1,3.81}},{{3.21,7.62},{3.43,6.8},{3.72,5.7},{4.04,4.81},{4.32,3.99},{4.67,3.39},{4.94,2.97},{5.29,2.85},{5.51,2.77},{5.95,3.16},{6.05,3.36}},{{6.79,2.11},{6.98,2.32},{7.2,2.6},{7.66,2.62},{7.83,2.71},{8.21,2.63},{8.5,2.55},{8.62,2.34},{8.97,2.04}},{{7.63,4.03},{7.93,4.18},{8.2,4.02},{8.49,3.87},{8.77,3.46},{9.22,3.13},{9.35,2.51},{9.61,2.21},{9.95, 1.86}}};


flinty suggested one technique, whereby data could be attached in order (say from left-to-right), and I've attempted a quick and dirty rendition of this:

SortedData=SortBy[Example,First];(*Sort by minimum x position*)
Result=SortedData[[1]];(*Rescaled Final Data is initially the first dataset*)
(*Find overlapping points of next set to final set*)
(*If no points overlap, extrapolate and fit the nearest points at each end*)
(*If there is overlap, then only use that to fit*)
(*Collect rescaled data together*)

flinty suggestion

This result does a pretty good job, although it has two possible issues:

  1. Fitting one additional curve at a time has trouble with regions where more than two curves overlap. This can be seen in the region around (x=5), where there is more noise compared to the same region fit by eye.

  2. Interpolation requires nonduplicate input, so data with the same x values cannot be interpolated together. I have gotten around this by simply averaging the scaled y-value when x is the same, but I expect that this may not be the best option.


aooiiii had a great approach, and I modified it a bit as QuadraticOptimization is a newer function that I can't use at home. This uses NMinimize to minimize the error in scaling parameters (s) of the log-data, while regularizing the function (y) in several possible ways, using simple approximations of first ("flat"), second ("smooth") and third ("jerk") derivatives at neighboring points. The main difference is that while aooiiii used many y's spanning between gaps in data, this version uses the input x positions to assign y points. I found the best-looking results using the third derivative ("jerk"), so the other regularization terms are commented out.

ss=Array[s,Length[d]];(*Scaling parameters*)
sd=Flatten[MapThread[{#[[All,1]],Log[#[[All,2]]]+#2}\[Transpose]&,{d,ss}],1];(*Changing to a log scale so scaling can't approach zero*)
xs=Union[sd[[All,1]]];(*List of unique x-values*)
ys=Array[y,Length[xs]];(*Corresponding y-function*)
(*flat=Total[Function[{x1,y1,x2,y2},((y2-y1)/(x2-x1))^2]@@@Flatten[Partition[{xs,ys}\[Transpose],2,1],{{1},{2,3}}]];(*Differences of nearby y-values*)*)
(*smooth=Total[Function[{x1,y1,x2,y2,x3,y3},(((x2-x1)(y3-y2)-(x3-x2)(y2-y1))/((x3-x2)(x3-x1)(x2-x1)))^2]@@@Flatten[Partition[{xs,ys}\[Transpose],3,1],{{1},{2,3}}]];(*Differences of nearby slopes*)*)
jerk=Total[Function[{x1,y1,x2,y2,x3,y3,x4,y4},(((x3(y1-y2)+x1(y2-y3)+x2(y3-y1))/((x1-x2)(x1-x3))-(x4(y2-y3)+x2(y3-y4)+x3(y4-y2))/((x4-x2)(x4-x3)))/((x2-x3) (x4+x3-x2-x1)))^2] @@@Flatten[Partition[{xs,ys}\[Transpose],4,1],{{1},{2,3}}]];(*Differences of nearby curvature*)
errors=Total[((sd[[All,1]]/.Rule@@@({xs,ys}\[Transpose]))-sd[[All,2]])^2];(*Differences of function to data*)
fit=NMinimize[(*flat/100+smooth/100+*)jerk/1000+errors/.s[1]->0,Join[ys,ss[[2 ;;]]]][[2]];(*Minimize all differences*)
stitched={xs,Exp[ys]}\[Transpose]/.fit;(*The optimized function*)
MapThread[{#[[All,1]],#[[All,2]]*#2}\[Transpose]&,{d,Exp[ss]}]/.s[1]->0/.fit(*Rescaled data*)]

Grid[{{"Initial Data","Final Scaled Data"},{ListLinePlot[Example,ImageSize->250],Show[ListLinePlot[Stitch[Example],ImageSize->250],ListPlot[stitched,PlotStyle->Directive[PointSize[0.02],Black]]]}}]

Final Rescaling

  • 1
    $\begingroup$ Your problem is solvable only if all measurements have overlapping parts I think. In your example it seems to be possible to group the three left ones and the two right ones. How to combine them "continous" ??? $\endgroup$ May 20, 2020 at 10:09
  • 3
    $\begingroup$ If you know the order the pieces fit together from left to right, and assuming the first one is fixed, then you could slide a piece along looking for maximum 'agreement' with the current fixed pieces, adding it to the fixed set, and repeat until all pieces are fixed. $\endgroup$
    – flinty
    May 20, 2020 at 10:33
  • $\begingroup$ Are they scaled via constant or do they have constants in x and y added to them aswell? $\endgroup$ May 20, 2020 at 11:02
  • 1
    $\begingroup$ Nice! Since advising is so much easier than coding, a couple more suggestions. 1 - make sure that your regularization term is approximately invariant to the grid, i.e. it's an approximation of some functional of the smooth function, independent of the x coordinate. 2 - with a uniform grid and interpolation in the error term, achieving (1) is probably the easiest. Linear is fine, but splines allow for a coarser grid. $\endgroup$
    – aooiiii
    May 26, 2020 at 11:54
  • 1
    $\begingroup$ Finally, you must have noticed that the estimated scaling varies a lot with the model and hyperparameter choice, especially near the gap. It would be nice to determine the hyperparameters automatically and compute uncertainties with a Bayesian approach. Sadly, Mathematica can't do this yet, but there's a Mathematica interface for Stan, which is simply the best for Bayesian modelling. Learning the basics of it takes about an evening. $\endgroup$
    – aooiiii
    May 26, 2020 at 11:56

2 Answers 2


A quick and dirty proof-of-concept implementation of my QuadraticOptimization idea. I haven't given it much thought, and the algorithm may require improvements, such as irregular grid, logarithmic scale, deciding how much and what type of smothness penalty is needed etc. The part I'm unsure about the most is requiring the smoothed curve to be above 1. There are probably better ways to prevent the optimizer from setting all of the scaling coefficients to 0, thus pointlessly achieving zero smoothness penalty and zero error.

data = Map[{Round[100 #[[1]]], #[[2]]} &, Example, {2}];
{min, max} = MinMax[Map[First, data, {2}]];

smoothness = Total@Table[(y[i] - 2 y[i + 1] + y[i + 2])^2, {i, min, max - 2}];
(*C2 smoothness penalty. One might combine several types of them here.*)

error = Total@Flatten@Table[
     (y[data[[i, j, 1]]] - s[i] data[[i, j, 2]])^2,
     {i, Length[data]},
     {j, Length[data[[i]]]}];

constr = Table[y[i] >= 1, {i, min, max}];

vars = Join[
   Table[y[i], {i, min, max}],
   Table[s[i], {i, Length[data]}]

sol = QuadraticOptimization[1000 smoothness + error, constr, vars];

patches = Table[{data[[i, j, 1]], data[[i, j, 2]] s[i]},
    {i, Length[data]},
    {j, Length[data[[i]]]}] /. sol;
smoothed = Table[{i, y[i]}, {i, min, max}] /. sol;

  ListPlot[patches, Joined -> True], 
  ListPlot[smoothed, Joined -> True, 
   PlotStyle -> {Opacity[0.1], Thickness[0.05]}]

enter image description here

  • $\begingroup$ It would seem that keeping the smoothed curve above 1 could just be achieved at the end of the process by finding the minimum of the smoothed curve and adding Boole[ymin < 1]*(1-ymin). $\endgroup$
    – JimB
    May 20, 2020 at 15:37
  • $\begingroup$ Very cool! I think this does exactly what I want, although unfortunately my home computer only has v11.3, while QuadraticOptimization is from v12.0. (I have access to v12 at an office though, so I should be able to test it) I think setting one of the s[i] coefficients to 1 (i.e. not scaling one set of data) may also offer a reasonable constraint for optimization. $\endgroup$
    – ChaSta
    May 21, 2020 at 8:45
  • 1
    $\begingroup$ @ChaSta Thanks! However, I must advise against fixing one coefficient. The optimizer is urged to bring the function as close to zero as possible: the smaller the scale, the lesser the penalties. Fixing the leftmost coefficient would result in a biased estimation, producing curves leaning to bottom-right. My original solution probably has some less pronounced bias too, but at least it's location independent. Switching to logarithmic coordinates would completely eliminate the bias, but make the lower regions more important, which may or may not be desirable. $\endgroup$
    – aooiiii
    May 21, 2020 at 16:44
  • 1
    $\begingroup$ Log scale also eliminates the need for constraints, reducing the optimization part to a trivial linear algebra problem. In Mathematica 11.3, you can construct the relevant matrix manually or use FindMinimum/NMinimize - they work much faster on unconstrained problems. $\endgroup$
    – aooiiii
    May 21, 2020 at 17:16
  • $\begingroup$ I have updated the above with a log-scale method that is a linear algebra problem. Thanks again! $\endgroup$
    – ChaSta
    May 26, 2020 at 8:39

Here is an approach that estimates the multiplicative constants by taking the log of the response variable and estimates the resulting additive constants.

(* Take the log of the response so that the adjustment is additive 
   and include the adjustments for each set of data *)
(* Force the last data set to have an adjustment of 0 *)
data2 = data;
n = Length[data];
adj[n] = 0;
data2[[All, All, 2]] = Log[data[[#, All, 2]]] + adj[#] & /@ Range[Length[data]];

(* Determine the binning parameters *)
{xmin, xmax} = MinMax[data[[All, All, 1]]];
nBins = 20;
width = (xmax - xmin)/nBins;

(* Calculate total of the variances *)
t = Total[Table[Variance[Select[Flatten[data2, 1], 
  -width/2 <= #[[1]] - xmin - (i - 1) width <= width/2 &][[All, 2]]] /. Abs[z_] -> z,
  {i, 1, nBins + 1}]] /. Variance[{z_}] -> 0;

(* Minimize the total of the variances and plot the result *)
sol = FindMinimum[t, Table[{adj[i], 0}, {i, n - 1}]]
(* {0.0518024, {adj[1] -> 0.510144, adj[2] -> -0.157574, adj[3] -> -0.352569, adj[4] -> 0.447345}} *)

(* Plot results on original scale *)
data3 = data2;
data3[[All, All, 2]] = Exp[data2[[All, All, 2]] /. sol[[2]]];
ListPlot[data3, Joined -> True, PlotLegends -> Automatic]

Adjusted data

  • $\begingroup$ This certainly does a great job for overlapping data, but doesn't adjust for the gap region around 6.5. $\endgroup$
    – ChaSta
    May 21, 2020 at 8:15
  • $\begingroup$ One can't adjust for the gap just based on the data. @UlrichNeumann already mentioned this in his comment. It's logically impossible. However, you have a curve form in your mind that would fill in the gap. But as good as some of the folks are here, they can't read your mind. You would need to give more information. That some might feel they can connect the two sets of segments as you want is only luck or wishful thinking without that additional information. $\endgroup$
    – JimB
    May 21, 2020 at 14:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.