Rescaling data to make a "continuous" function

Question

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:

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}}};

UPDATE

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*)
For[i=2,i<=Length[SortedData],i++,
OverlappingPoints=Select[SortedData[[i]],#[[1]]<=Max[Result[[All,1]]]&];
(*Find overlapping points of next set to final set*)
Scaling=If[OverlappingPoints=={}, 
NArgMin[(Interpolation[Result][SortedData[[i,1,1]]]-s*SortedData[[i,1,2]])^2+(s*Interpolation[SortedData[[i]]][Result[[-1,1]]]-Result[[-1,2]])^2,s],
(*If no points overlap, extrapolate and fit the nearest points at each end*)
NArgMin[Total[(Interpolation[Result][#[[1]]]-s*#[[2]])^2&/@OverlappingPoints],s]];
(*If there is overlap, then only use that to fit*)
Result=Sort[Mean/@GatherBy[Join[Result,{1,Scaling}*#&/@SortedData[[i]]],First]]] 
(*Collect rescaled data together*)
ListLinePlot[Result,PlotStyle->Black]

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.

SECOND UPDATE

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.

Stitch[d_]:=Module[{ss,sd,flat,smooth,jerk,errors,fit},
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]]]}}]

Answer 1

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}]];
(*Discretizing*)

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;

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

Answer 2

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]