6
$\begingroup$

I am trying to fit a data set to another one. There must be a couple of threads that treat this issue but I have either not found them or was unable to adapt the examples to my specific case.

Here is the drill: I collect data from a spectrometer in front of which I have a slit - of which I control the width. I take the data for the smallest slit width and smooth it with a Gaussian (via a convolution product) and subsequently use this "smoothed data set" as my theory curve. In other words, I want to fit this smoothed data set to the sets for wider slit width by controlling the width and height of the Gaussian function used to smooth it.

I don't really want to find the best fit for each data set and compare them, but rather adapt the smoothing of the first data set to best fit the following ones.

My data sets are two dimensional lists with the wavelength in one column and the photon counts in another one (all real numbers).

The narrowest slit data is:

In[242]:= ramp001mm
Out[242]= {{190.999,560},{191.019,482},{191.039,569},{191.058,532},{191.078,564},{191.098,594},{191.118,561},{191.139,592},{191.159,557},{191.179,556},{191.199,566},{191.219,579},{191.239,581},{191.258,584},{191.278,576},{191.298,532},{191.318,570},{191.338,513},{191.359,528},{191.379,537},{191.399,585},{191.419,620},{191.439,633},{191.459,606},{191.478,673},{191.498,590},{191.518,591},{191.538,600},{191.559,550},{191.579,485},{191.599,466},{191.619,485},{191.639,455},{191.659,427},{191.678,400},{191.698,427},{191.718,393},{191.738,368},{191.759,368},{191.779,275},{191.799,301},{191.819,294},{191.839,261},{191.859,248},{191.879,226},{191.898,252},{191.918,180},{191.938,199},{191.958,200},{191.979,197},{191.999,172},{192.019,188},{192.039,188},{192.059,189},{192.079,183},{192.098,165},{192.118,175},{192.138,184},{192.158,187},{192.178,172},{192.199,154},{192.219,141},{192.239,130},{192.259,122},{192.279,103},{192.299,106},{192.318,96},{192.338,94},{192.358,78},{192.378,79},{192.399,72},{192.419,84},{192.439,92},{192.459,65},{192.479,85},{192.499,73},{192.519,61},{192.538,55},{192.558,57},{192.578,55},{192.598,59},{192.619,65},{192.639,59},{192.659,60},{192.679,77},{192.699,59},{192.719,55},{192.738,70},{192.758,65},{192.778,58},{192.798,62},{192.818,66},{192.839,64},{192.859,94},{192.879,62},{192.899,63},{192.919,66},{192.939,63},{192.959,58},{192.978,64},{192.998,84},{193.018,45},{193.038,67},{193.059,43},{193.079,58},{193.099,55},{193.119,80},{193.139,60},{193.159,50},{193.179,48},{193.198,62},{193.218,39},{193.238,45},{193.258,49},{193.279,63},{193.299,64},{193.319,48},{193.339,50},{193.359,60},{193.379,57},{193.398,37},{193.418,59},{193.438,61},{193.458,63},{193.478,46},{193.499,47}}

I use ListConvolve to smooth the above data set with a Gaussian in the following way:

ListConvolve[h Table[ Exp[-s^2/\[Sigma]^2]/Sqrt[ 2Pi],{s,-17,17}],ramp001mm[[All,2]]]

The $h$ and $\sigma$ parameters are the height, and width, of the Gaussian, respectively. I have arbitrarily set the number of points the Gaussian spans to 34 (s runs from -17 to 17). I want both $h$ and $\sigma$ to be (real valued) irrational numbers (i.e. I don't necessarily want them to be integers).

My final aim is to fit the above smoothed data to the following data (for example):

In[243]:= ramp002mm
Out[243]= {{190.999,1686},{191.019,1693},{191.039,1656},{191.058,1645},{191.078,1645},{191.098,1626},{191.118,1583},{191.139,1586},{191.159,1644},{191.179,1578},{191.199,1658},{191.219,1783},{191.239,1728},{191.258,1650},{191.278,1613},{191.298,1573},{191.318,1566},{191.338,1527},{191.359,1562},{191.379,1606},{191.399,1687},{191.419,1823},{191.439,1841},{191.459,1909},{191.478,1872},{191.498,1697},{191.518,1709},{191.538,1618},{191.559,1571},{191.579,1502},{191.599,1366},{191.619,1305},{191.639,1237},{191.659,1222},{191.678,1193},{191.698,1196},{191.718,1134},{191.738,1082},{191.759,992},{191.779,886},{191.799,874},{191.819,853},{191.839,741},{191.859,724},{191.879,666},{191.898,591},{191.918,582},{191.938,520},{191.958,526},{191.979,502},{191.999,491},{192.019,518},{192.039,514},{192.059,605},{192.079,545},{192.098,528},{192.118,486},{192.138,481},{192.158,515},{192.178,387},{192.199,388},{192.219,360},{192.239,278},{192.259,311},{192.279,295},{192.299,269},{192.318,245},{192.338,232},{192.358,217},{192.378,205},{192.399,206},{192.419,184},{192.439,175},{192.459,190},{192.479,170},{192.499,175},{192.519,169},{192.538,147},{192.558,166},{192.578,161},{192.598,155},{192.619,160},{192.639,156},{192.659,161},{192.679,150},{192.699,138},{192.719,135},{192.738,143},{192.758,126},{192.778,165},{192.798,123},{192.818,116},{192.839,149},{192.859,139},{192.879,143},{192.899,151},{192.919,123},{192.939,141},{192.959,118},{192.978,122},{192.998,141},{193.018,146},{193.038,135},{193.059,146},{193.079,126},{193.099,125},{193.119,134},{193.139,117},{193.159,111},{193.179,136},{193.198,135},{193.218,115},{193.238,118},{193.258,121},{193.279,123},{193.299,116},{193.319,125},{193.339,112},{193.359,113},{193.379,116},{193.398,110},{193.418,115},{193.438,109},{193.458,103},{193.478,109},{193.499,102}}

And I don't really know how to ask Mathematica to do that for me. So far I have been doing the fit manually, by visually assessing the quality of the fit when plotting the smoothed data (ramp001mm smoothed) on top of the raw data (ramp002mm) and adjusting $h$ and $\sigma$ accordingly. I would want to either find a built-in function of Mathematica that would do this (NonlinearModelFit?) or write a function myself (that would minimise the least-squared values and residuals by exploring the parameter space {$h.\sigma$}.

EDIT:

Clearly, ListConvolve returns a list that is 34 elements shorter than the input ramp001mm (17 elements shorter from the beginning and from the end of the input list), that is a 92 elements long list. I thus try and fit this 'theory curve' to a shortened version of the other lists (i.e. to ramp002mm[[1 + 17 ;; Length[ramp002mm] - 17, 2]] in this example).

If you need any clarification, please let me know!

Improve this question
$\endgroup$
5
  • 1
    $\begingroup$ I might not fully grasp what you are after, so: you create a set of functions conv=ListConvolve[...] that depend on the same $h$ and $\sigma$ each and then want to find particular values of $h$ and $\sigma$ so that the conv will overlay the set of points ramp002mm best? Did you note that you have 92 points in conv and 126 in ramp002mm? $\endgroup$
    – corey979
    Sep 9, 2016 at 12:14
  • $\begingroup$ Yes @corey979, that's exactly what I want: explore the 2-D parameter space ($h$ and $\sigma$) to best fit conv to ramp002mm. The list Length changes because of the smoothing by the Gaussian - that is defined over 34 points: it deletes 17 points from the beginning and end of the original data set. $\endgroup$
    – MXJ
    Sep 9, 2016 at 12:18
  • 1
    $\begingroup$ conv gives a set of 92 points; ramp002mm consists of 126 points - these numbers don't match. ListConvolve gives a new list of points dependent on two parameters, but it's shorter than the original list. I don't see a way t o fit two datasets of unequal lengths; maybe someone will find a tick to do so. $\endgroup$
    – corey979
    Sep 9, 2016 at 12:34
  • $\begingroup$ Sorry, I see why this is unclear. As I said, ListConvolve returns a 92 points long list, so what I do is that I discard the 17 first and 17 last points of ramp002mm, thus fitting conv to ramp002mm[[1 + 17 ;; Length[ramp002mm] - 17, 2]]. At the moment I do this fit by plotting the two lists (conv and the shortened ramp002mm) and visually adjusting $h$ and $\sigma$. I would want Mathematica to adjust these two parameters automatically for me. $\endgroup$
    – MXJ
    Sep 9, 2016 at 13:37
  • $\begingroup$ At some point if you want to make inferences about $h$ and $\sigma$, you'll need to write down an explicit model that characterizes the relationship including the description as to how deviations (sometimes labeled as random errors) are modeled. The information about the underlying model you have in mind is not an intrinsic feature of the data. $\endgroup$
    – JimB
    Sep 9, 2016 at 13:48

2 Answers 2

Reset to default
7
$\begingroup$

Create a smoothed dataset dependent on the parameters $h$ and $\sigma$:

conv = ListConvolve[
h Table[Exp[-s^2/σ^2]/Sqrt[2 π], {s, -17, 17}], 
ramp001mm[[All, 2]]];

and a shortened version of the second dataset:

ndat = ramp002mm[[1 + 17 ;; Length[ramp002mm] - 17, 2]];

Then one can perform a minimization of the sum of squares in a straightforward manner:

NMinimize[Total[(conv - ndat)^2], {h, σ}]

{170287., {h -> 2.69303, σ -> 1.50493}}

which gives an answer in an instant. To directly assign values to $h$ and $\sigma$:

{h, σ} = {h, σ} /. 
NMinimize[Total[(conv - ndat)^2], {h, σ}][[2]]

{2.69303, 1.50493}

The output can be plotted:

enter image description here

where blue is ndat and red is the conv with fitted $h$ and $\sigma$.

For conv I took exactly the form provided in the question. However, to be a true Gaussian there should be a multiplicative $1/\sigma$, i.e. (σ Sqrt[2 π]) should replace Sqrt[2 π]. After this replacement, one gets

{h, σ} = {1.07352, 0.149686}

and the plot looks like this:

enter image description here

EDIT: Another approach might be minimizing a $\chi^2=\sum\limits_{i=1}^N\frac{(E_i-O_i)^2}{E_i}$ instead of the sum of squares. The true Gaussian smoothing is performed:

conv = ListConvolve[
   h Table[Exp[-s^2/σ^2]/(σ Sqrt[2 π]), {s, -17, 
      17}], ramp001mm[[All, 2]]];

min = NMinimize[Total[(conv - ndat)^2/conv], {{h, 0., 5.}, {σ, 0., 5.}}]

{483.825, {h -> 3.93629, σ -> 1.83513}}

chi2 = First@min
{h, σ} = {h, σ} /. min[[2]]

A global NMinimize doesn't converge here so an interval for $h$ and $\sigma$ is required. Moreover, convergence wasn't met when the intervals were {0., 10.} so I narrowed them to {0., 5.}. The parameters attained different values than before because a different function was minimized. The result is as follows:

enter image description here

The value $\chi^2=483.825$ can be used as a statistical diagnostic; it's common to report a $\chi^2/dof$ value, although I'm not sure what the number of degrees of freedom should be for this setup.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Brilliant, thanks a lot @corey979! I didn't manage to use NMinize this way, because I didn't understand what I had to minimize was the error function. A pity NonLinearModelFit cannot be used for the present purpose for I would have liked to use the diagnostics tools it offers. $\endgroup$
    – MXJ
    Sep 9, 2016 at 15:17
11
$\begingroup$

"A pity NonLinearModelFit cannot be used for the present purpose..."

Don't count out NonlinearModelFit (or its cousins) yet.

First a plot of the two datasets together:

ListPlot[{d2, d1}, Joined -> True,
 PlotLegends -> LineLegend[Automatic, {"Wider slit data", "Narrowist slit data"}]

Two datasets]

But this doesn't directly get at the relationship that exists between the two photon counts. A better figure for that purpose is the following:

ListPlot[Transpose[{d1[[All, 2]], d2[[All, 2]]}],
 AxesLabel -> {"Photon counts\nNarrowist slit data", 
   "Photon counts\nWider slit data"}]

Relationship of photon counts

This looks pretty linear to me. We use NonlinearModelFit (even though initially this is a bit of overkill) and examine the fit:

nlm = NonlinearModelFit[Transpose[{d1[[All, 2]], d2[[All, 2]]}],
   a + b PhotonCount, {a, b}, PhotonCount];
nlm["BestFitParameters"]
(* {a->-41.30773326085267`,b->2.9783305744372406`} *)

(* Predicted vs. residual *)
ListPlot[Transpose[{nlm["PredictedResponse"], nlm["FitResiduals"]}],
 Frame -> True, 
 FrameLabel -> {{"Residual", ""}, {"Predicted", 
    "Predicted vs. Residual"}}]
(* Quantile plot *)
QuantilePlot[nlm["FitResiduals"]]
(* Histogram of residuals *)
Show[Histogram[nlm["FitResiduals"], "FreedmanDiaconis", "PDF", 
  Frame -> True,
  FrameLabel -> {{"Probability density", ""}, {"Residual", 
     "Histogram"}}],
 Plot[PDF[
    NormalDistribution[0, StandardDeviation[nlm["FitResiduals"]]]][
   x], {x, -150, 150}]]

Residual plots

The residual plots don't look so great. But with count data sometimes taking the square root works better (by equalizing the variance about the fitted curve). So we fit a relationship using the square roots of the counts:

nlm = NonlinearModelFit[
   Transpose[{d1[[All, 2]]^0.5, d2[[All, 2]]^0.5}],
   a + b PhotonCount, {a, b}, PhotonCount];
nlm["BestFitParameters"]
(* {a->-2.1299900144167685`,b->1.7962540226784693`} *)

(* Predicted vs. residual *)
ListPlot[Transpose[{nlm["PredictedResponse"], nlm["FitResiduals"]}],
 Frame -> True, 
 FrameLabel -> {{"Residual", ""}, {"Predicted", "Predicted vs. Residual"}}]
(* Quantile plot *)
QuantilePlot[nlm["FitResiduals"],
 FrameLabel -> {{"Empirical Quantiles", ""}, {"Theoretical Quantiles",
     "QQ Plot"}}]
(* Histogram of residuals *)
Show[Histogram[nlm["FitResiduals"], "FreedmanDiaconis", "PDF", 
  Frame -> True,
  FrameLabel -> {{"Probability density", ""}, {"Residual", 
     "Histogram"}}],
 Plot[PDF[
    NormalDistribution[0, StandardDeviation[nlm["FitResiduals"]]]][
   x], {x, -4, 4}]]

Residuals from square root transformation

These residual plots look much better. Here's the original "wider slit data" and the prediction:

ListPlot[{d2, Transpose[{d2[[All, 1]], nlm["PredictedResponse"]^2}]},
 Joined -> {False, True}, Frame -> True,
 FrameLabel -> {{"Photon count", ""}, {"Wavelength", ""}}, 
 PlotLegends -> LineLegend[Automatic, {"Wider slit data", 
    "Predicted from narrower slit data"}]]

Data and fit with square root transformation

I'm not convinced that you need to do any smoothing (or without an expected theoretical curve why even do it?) But you could smooth either before or after if that was really needed. (It's possible that GeneralizedLinearModelFit might do a bit better as it does seem that you are attempting to predict a Poisson mean from Poisson counts. The implies that the predictor variable is random rather than fixed as assumed by the regression approach.)

Improve this answer
$\endgroup$
4
  • 1
    $\begingroup$ beautiful exposition +1 o/c :) $\endgroup$
    – ubpdqn
    Sep 10, 2016 at 7:34
  • $\begingroup$ Wow! Awesome answer. Very clear helpful explanation. $\endgroup$
    – Jack LaVigne
    Sep 11, 2016 at 14:07
  • $\begingroup$ Thanks @Jim! :-) This is a very detailed and clear explanation. The whole point of my procedure is to smooth the first data set with a Gaussian because I expect this function models the resolution of my spectrometer well - in other words, the $\sigma$ of the Gaussian is the spectral resolution of my spectrometer. Fitting the smoothed narrowest slit data to data for wider slits opening allows me to find the spectral resolution of the spectrometer as a function of the slit opening. $\endgroup$
    – MXJ
    Sep 12, 2016 at 9:03
  • $\begingroup$ As the others said, great post! $\endgroup$
    – Anton Antonov
    Oct 25, 2016 at 11:00

Your Answer

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

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