14

I have a dataset of which I'm attempting to determine the full width at half maximum. A picture of the data is given below, and the dataset itself can be found here (pastebin). enter image description here

Now, this set has a minimum, which can be found rather easily. Say we call the data yy, then by using Min[yy] we've found the minimum. To find the full width half max, it is then probably easiest to look for a point closest to 1-(1-Min[yy])*0.5 in yy. The first thing that came to mind was using Nearest[yy,1-(1-Min[yy])*0.5,1] which indeed returns a value of yyclosest to half of the FWHM. But I'm having some trouble finding the actual position of this point in the data. Ironically, using Position[yy,Nearest[yy,1-(1-Min[yy])*0.5,1]] seems to come up empty. I'm not sure why exactly this is, but maybe there is some rounding going on, or incompatible outputs. I've tried thinking of a solution to this, but it doesn't seem obvious to me.

My question is therefore if anyone seems a workaround, or an all around better idea. My plan was to find this position of the half-max point, find the corresponding xx, subtract the xx element of the position of the max, and multiply by two. I don't need a 100% perfect estimation; this should be good enough.

Community
  • 1
user129412
  • 1,339
  • 8
  • 19
  • I found a solution actually, but it's quite simple. What should I do, answer my own question, or simply delete it? I used the answer of http://mathematica.stackexchange.com/questions/38570/efficient-way-to-give-nearest-neighbour-and-next-nearest-neighbour-of-every-poin where Rahul wrote nf = Nearest[Table[ps[[i]] -> i, {i, Length[ps]}]] (* so it returns the index of the site *) which is exactly what I needed to get Nearest working in my function. – user129412 Jun 15 '15 at 00:09
  • 1
    Answering your own question is encouraged,http://mathematica.stackexchange.com/help/self-answer. – image_doctor Jun 15 '15 at 00:43

2 Answers2

15

Here is an approach. You could also search around the midpoint. In the following data was imported using ImportString from pastebin link.

A non-linear model was developed. The starting points were 'guessed' with the assistance of manipulate.

nlm = NonlinearModelFit[N@ToExpression[data], 
  d - c/(1 + a (x - b)^2), {{a, 150000.}, {b, 7.62}, {c, 0.064}, {d, 
    1.}}, x]
Show[Plot[nlm[x], {x, 7.6, 7.65}, PlotRange -> {0.90, 1}, 
  PlotStyle -> Red], 
 ListPlot[data, PlotMarkers -> {Style[\[FilledDiamond], Black], 5}]]

enter image description here The fit (visually is not unreasonable, adjusted $R^2$ 0.999989).

Now to determine width.

bf = nlm["BestFitParameters"]
ypeak = nlm[x] /. First@Quiet[Solve[D[nlm[x], x] == 0, x]];
hmw = ypeak + ((d /. bf) - ypeak)/2;
{xr, xl} = x /. NSolve[nlm[x] == hmw, x];
xr - xl
Show[Plot[nlm[x], {x, 7.6, 7.65}, PlotRange -> {0.90, 1}, 
  PlotStyle -> Red, 
  Epilog -> {Green, PointSize[0.02], Point[{#, nlm@#} & /@ {xl, xr}]},
   GridLines -> {{xl, xr}, {d /. bf, hmw, ypeak}}, 
  GridLinesStyle -> Directive[Gray, Dashed]], 
 ListPlot[data, PlotMarkers -> {Style[\[FilledDiamond], Black], 5}]]

The x values: {7.61527, 7.61774}. The width: 0.00247168.

enter image description here

This can adpated, improved and I look forward to better answers.

ubpdqn
  • 60,617
  • 3
  • 59
  • 148
  • Interesting approach! I'll soon have time to post my answer, which if I remember correctly finds about the same value. Requires a little bit less manual input than fitting it, but certainly a solid option if one only has to do it for a few datasets and probably a bit more accurate than mine. – user129412 Jun 15 '15 at 11:24
  • May I ask how the nonlinear model is selected in the first place? – Kattern Jun 20 '15 at 01:26
  • 2
    @Kattern as Box is quoted: 'all models are wrong but some are useful' and others (such as Nassim Taleb) have extended, '...and some are dangerous', I merely guessed a symmetric 'bell shaped' model and tweaked parameters to fit...in the absence of any understanding of the data, any important considerations of plausibility of such a model for process producing the data, it was just a a shape matching exercise to facilitate measurement...which is what I understood as the objective... – ubpdqn Jun 20 '15 at 05:06
  • Thanks for you explanation! I find the part guessing a model is not an easy task. Thus, I asked about the method of choosing the objective function. – Kattern Jun 20 '15 at 05:10
4

My approach, based on part of the code in this post. Lets call the list of (x,y) data data. Then it is actually only two lines of code

half = Nearest[Table[data[[i, 2]] -> i, {i, Length[data]}], 
1 - 0.5*(1 - Min[data[[All, 2]]]), 1];
FWHM = First[Abs[data[[half, 1]] - data[[First[Position[data[[All, 2]],
Min[data[[All, 2]]]]], 1]]*2];

It is not all that elegant, but it works. half looks for the index of the point in the y data that is closest to half the minimum, and then FWHM calculates the difference between the x coordinate belonging to half the minimum and the actual minimum, and then multiplies by two. I have to use First because some of the output comes in angular brackets. Now, the whole construction is not very nice, I'm sure I could have used the commands in a more efficient way, so feel free to point that out.

The value that I get for my data: 0.00266667, which we can compare to ubpdqn's answer, who found 0.00247168. Same order of magnitude for sure. One point about this specific approach is that we can not simply take the half minimum on either side of the minimum, because it is not guaranteed that there aren't two points on the right side of the minimum that are closer to half of the minimum, than the closest point on the left. We therefore only take 1 point and multiply the distance by two, assuming symmetry. Moreover, this is prone to noise/inaccurate data, as the point we find might just by accident be close to the half minimum. On the other hand, one does not need to go through NonlinearModelFit, which requires initial conditions to be picked by hand.

Community
  • 1
user129412
  • 1,339
  • 8
  • 19
  • Have you tried adding random errors to your points and see how stable is your method of calculating FWHM? – BlacKow Jun 15 '15 at 15:31
  • No, I'll try and experiment with that later today. – user129412 Jun 15 '15 at 15:33
  • Generally speaking there are two ways of doing that: 1) fit with predefined function (good candidates at http://mathworld.wolfram.com/FullWidthatHalfMaximum.html) and then use FWHM calculated from fitted parameters - this is what @ubpdqn did or 2) if you have no idea about the function, you can cutoff all points outside of your peak and calculate standard deviation statistically. – BlacKow Jun 15 '15 at 15:37
  • Hm, I might misunderstand your original comment then. Couldn't I just add some gaussian noise and see how that influences the result? On a different note, I have 8 of these datasets and the value should be the same throughout all them. Some initial testing does seem to confirm that they are close. – user129412 Jun 15 '15 at 16:46
  • I'm just saying that your method is unstable and won't give robust results. I suggested you to prove it yourself by adding random noise to your dataset. After you try it you can use fitting and see if it's better. Also keep in mind that blurring out noise will affect (increase) FWHM, so you can't just get rid of noise the same way we discussed in your previous question about peak finding. – BlacKow Jun 15 '15 at 16:52
  • I see what you mean. In that case I will accept the first answer, as you make a good point on the stability. For me my method will be preferable, as I simply use it to get an initial value for a nonlinear fit. If I have to go through a nonlinear fit of the Peak to get an initial value, well, that's not going to be incredibly efficient. I didn't write this in my main post as I'm trying to keep it general. – user129412 Jun 15 '15 at 16:57