Working with Interpolating Functions

I have two sets of data that are not regularly sampled, and I want to produce a quantitative measurement of how similar they are while fitting them with an $x$ shift and $x$ and $y$ scaling.

Following the Mathematica interpolation documentation, f = Interpolation[data] gives me nice graphs and allows me to access one interpolated point at a time via f[x].

What I'd like to do is something along the lines of:

experiment = Interpolation[expdata]
FindFit[theorydata, a*experiment[b*x + c], {a, b, c}, x]


where theorydata and expdata are non-regularly sampled lists of $x$ and $y$ values.

This successfully inputs $x$ values from the theorydata list, but just prints them out as FindFit[{theorydata}, {a InterpolatingFunction[{{0., 0.32878}}, <>][0.499755 b + c] ...]

Alternatively, if something like that is not an option, it would be useful to know how to generally work with interpolating functions. For example, I cannot seem to subtract two interpolating functions from each other:

experiment = Interpolation[expdata]
theory = Interpolation[theorydata]
f = experiment - theory


This outputs:

InterpolatingFunction[{{0., 0.32758}}, <>] - InterpolatingFunction[{{0., 0.32878}}, <>]


and, f[0.1], for example outputs:

InterpolatingFunction[{{0., 0.32758}}, <>] - InterpolatingFunction[{{0., 0.32878}}, <>][0.1]


I am, however, able to re-sample the interpolated data with:

result = Table[{x, experiment[x]}, {x, 0, 0.3, 0.01}]


This then gives me a regularly sampled list.

In summary, I suppose the biggest question here is how to work with InterpolatingFunction objects. Can they be added/subtracted? Can they be treated as a function and used in fitting algorithms? How exactly is an InterpolatingFunction stored in Mathematica?

EDIT: Fixed a typo in the output

Sample of data (small sample of a large number of points):

    expdata={{0, 0.032174}, {0.00497, 0.032446}, {0.010701, 0.032923}, {0.015671,
0.033402}, {0.021403, 0.034199}, {0.028477, 0.035362}, {0.033447,
0.036444}, {0.039179, 0.038074}, {0.044149, 0.039642}, {0.049119,
0.041422}}
theorydata={{0, 0.033955}, {0.00497, 0.034502}, {0.010701, 0.035312}, {0.015671,
0.036135}, {0.020641, 0.037109}, {0.026373, 0.038481}, {0.031343,
0.039824}, {0.036313, 0.041344}, {0.041283, 0.043037}, {0.047014,
0.045256}}


Note that that the datasets are very similar, but they are not sampled for all of the same x values.

• Could you share a sample of your two datasets? Also, your idea of fitting might not be along the lines of what FindFit does, i.e. you need to have an analytical model and experimental data points. In your case, you could generate a function that expresses the squared differences between your experimental and corresponding theoretical data points as a function of the shift parameters, and then minimize that function numerically to get the best-fit shift parameters for this the theoretical calculations fit the experimental findings. – MarcoB Sep 30 '15 at 21:43
• @MarcoB You're partially right,but remember that "I have two sets of data that are not regularly sampled," So probably no correspondence – Dr. belisarius Sep 30 '15 at 22:05
• @MarcoB I added the first 10 datapoints for the theory and experimental datasets. Since I can use the Table function to resample the interpolated data, I could do a manual sum of squares minimization loop, altering the shift/scaling of the original datasets with each iteration. I would really like to avoid having to to that though. – Snipatomic Sep 30 '15 at 22:19
• @belisarius That's an excellent point: I had overlooked that complicating factor. Thank you for pointing that out. – MarcoB Oct 1 '15 at 1:35
• @Snipatomic Thank you for posting the data; please take a look at my answer below. I hope it is a step in the right direction. – MarcoB Oct 1 '15 at 1:36

First let us use your theoretical and experimental data points to generate interpolating functions:

(* Calculate interpolating functions for the two datasets *)
expint = Interpolation[expdata];
theorint = Interpolation[theorydata];


Now we use the interpolating functions to "resample" your data points; we can then express the point-by-point squared differences between the two data sets as a function of stretch and offset parameters $a$, $b$, and $c$.

(* Generate sum of squared distances as function of a, b, c to be minimized *)
Clear[diff]
diff[a_, b_, c_] = Total[Table[(a expint[b x + c] - theorint[x])^2, {x, 0, 0.047, 0.005}]];


We can minimize the squared difference function diff numerically to find the values of the parameters $a$, $b$, and $c$ that make expdata best approximate theorydata:

(* Find a, b, c that minimize diff *)
Off[InterpolatingFunction::dmval]
paramsatminimum = NMinimize[diff[a, b, c], {a, b, c}][[2]]
On[InterpolatingFunction::dmval]


During minimization NMinimize tries out some values of the parameters that bring the argument of the interpolating functions outside of their domain, which forces extrapolation, and a lot of warnings are generated. Since no harm is done by those attempts, I temporarily turn Off the warning message, carry out the optimization, then turn it back On.

Finally, we can plot the results by plugging the best-fit values of the parameters into the interpolating function expint, and compare them with the theorydata points:

(* Plot resulting interpolating function together with theorydata points *)
Plot[
a expint[b x + c] /. paramsatminimum, {x, 0, 0.05},
Epilog -> {PointSize[0.02], Red, Point@theorydata},
Frame -> True, Axes -> False, PlotRangePadding -> Scaled[.05]
]

(* Out: {a -> 1.01519, b -> 0.861444, c -> 0.0159744} *)


• Thank you so much for your help! I am having one problem currently: even when I run your code as written, I end up with a result that is wildly off. However, if I set the WorkingPrecision->10 in the NMinimize loop, I end up with accurate results. Is there a default setting that is messed up here? – Snipatomic Oct 1 '15 at 5:50
• @Snipatomic When you say that you run the code as written, are you using the small sample data set, or your complete set? What version of MMA do you run? I am on MMA 10.2, Win7-64. I ran my calculations without specifying any particular precision value, so they ran at machine precision; you shouldn't need to set anything either. Do you get any errors/warnings? – MarcoB Oct 1 '15 at 5:54
• I am running 10.2, Win10-64, Student Edition. My comment originally applied to the small dataset - it spat out {a -> 0.88840, b -> 0.61518, c -> 0.03802}, which is quite different from what you posted. However, with WorkingPrecision->10, I get exactly what you wrote. The only errors in the output are dmval (when I allow the warnings that you turned off). I have played with it a bit more now (both with the small dataset and my full set), and the NMinimize function seems to be extremely sensitive to the WorkingPrecision setting. Is there another way I can make this more robust? – Snipatomic Oct 1 '15 at 6:52