# Estimating the full width half maximum of a peak

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). 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.

• 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 mathematica.stackexchange.com/questions/38570/… 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. Jun 15, 2015 at 0:09

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}]] 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. • 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. Jun 15, 2015 at 11:24
• May I ask how the nonlinear model is selected in the first place? Jun 20, 2015 at 1:26
• @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... Jun 20, 2015 at 5: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. Jun 20, 2015 at 5:10

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.

• Have you tried adding random errors to your points and see how stable is your method of calculating FWHM? Jun 15, 2015 at 15:31
• No, I'll try and experiment with that later today. Jun 15, 2015 at 15:33
• Generally speaking there are two ways of doing that: 1) fit with predefined function (good candidates at 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. Jun 15, 2015 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. Jun 15, 2015 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. Jun 15, 2015 at 16:52