2
$\begingroup$

I have numeric data points (x,y) values from a peak. I constructed a interpolation to get an interpolating function.

Interpolation[data]

What I want to do is the following: Determine the width at 10% peak height. I looked for the peak maximum (y max value) and multiplied it by 0.1.

tenpercent = (datasets[[1]] // Max[#[[All, 2]]] &)*0.1

At this y value, I want to look for the two x values in the interpolating function that have this y value. From this I could easily determine the width of the peak (x2-x1). I tried using Reduce,FindRoot and NDSolve to find the two solutions of the interpolating function (x1,x2), but nothing worked. By searching through the threads, I found that my problem is more difficult than I thought. But I cannot believe that there is no easy way to solve this problem in Mathematica.

$\endgroup$
4
  • 2
    $\begingroup$ Please provide sample data to go with the rest of your code. $\endgroup$
    – Yves Klett
    May 8, 2015 at 13:24
  • $\begingroup$ I can't give you my sample data, but just imagine, that my interpolating function looks like a GaussNormalDistribution, e.g. Plot[PDF[NormalDistribution[5, 1.5]][x], {x, 0, 15}]; What I want is the width of the peak at 10% peak height $\endgroup$
    – Niki
    May 8, 2015 at 14:09
  • $\begingroup$ Let's say g is your interpolation. You might want to refine something on the lines of sampleddata = Table[{x, g[x]}, {x, a, b, step}]; Select[sampleddata, (Abs[(#[[2]] - tenpercent] < tol) &] $\endgroup$
    – Peltio
    May 8, 2015 at 14:24
  • 3
    $\begingroup$ Then please add that bit of code to generate the bogus data. This will make your question much more useful for visitors trying to evaluate your code, and should enhance chances to get good answers. $\endgroup$
    – Yves Klett
    May 8, 2015 at 14:26

3 Answers 3

5
$\begingroup$

It sounds like your problem would be resolved by defining the Interpolation itself as a function of x (or whatever variable you like.) This allows you to feed it into functions like FindMaximum and FindRoot much more easily. For example:

sampledata = Table[{x, PDF[NormalDistribution[5, 1.5]][x]}, {x, 0, 15, 0.5}]
interpfn[x_] = Interpolation[sampledata][x]
peakpos = FindMaximum[interpfn[x0], x0]
FindRoot[interpfn[x] == 0.1*peakpos[[1]], {x, (x0 /. peakpos[[2]]) + 0.05}]
FindRoot[interpfn[x] == 0.1*peakpos[[1]], {x, (x0 /. peakpos[[2]]) - 0.05}]

{0.265962, {x0 -> 5.}}
{x -> 8.21961}
{x -> 1.78039}

A quick check:

interpfn[x /. {x -> 8.21961}]
interpfn[x /. {x -> 1.78039}]

0.0265962
0.0265962

As desired, these are 10% of the peak value found above.

A few notes: first, the 0.05 values in the FindRoot functions are there to "push" the algorithm away from the peak in the positive or negative directions. You may want to tinker with this value for better results, depending on how quickly your data varies. Second, this code throws a surprising number of errors when I run it, largely due to FindRoot "guessing" values outside the range of the InterpolatingFunction. This means that it may not be terribly robust if your data doesn't vary as smoothly as a Gaussian. If this is the case, a method based on Select (as proposed by @Peltio in the comments) might be a better bet.

$\endgroup$
1
  • $\begingroup$ setting the InterpolationOrder to 1 to have a simple linear interpolation is one way to get rid of the surprising errors with FindRoot going outside the obvious range with FindMaximum $\endgroup$ May 8, 2015 at 15:49
3
$\begingroup$

Just for fun:

Using sampledata from other answers, you could approximate width from interpolation extracting mesh points (would not be robust for more complex continuous curves):

ip[x_] := Interpolation[sampledata][x];
mx = Quiet@FindMaximum[ip[x], x];
pl = Plot[ip[x], {x, 0, 15}, MeshFunctions -> (#2 &), 
   Mesh -> {{0.1 mx[[1]]}}];
Differences[(Part[pl[[1, 1]], ##] & @@ 
    Cases[pl, Point[x__] :> x, -1])][[1, 1]]

yielding approximate width: 6.44072 (vs 6.4379 with analytic and FindRoot). Could be useful in approximations for starting positions for numerical approximations for other data.

func[dist_, param_, p_] := 
 Quiet[{{##}, #2 - #1} & @@ (z /. 
     Solve[PDF[dist @@ param, z] == p PDF[dist @@ param, param[[1]]], 
      z])]
vis[d_, pm_, p_] := Module[{val, pts}, val = func[d, pm, p];
  pts = {#, PDF[d @@ pm, #]} & /@ val[[1]];
  Plot[PDF[d @@ pm, x], {x, 0, 1.5 val[[1, 2]]}, 
   Epilog -> {{Red, PointSize[0.02], Point[pts]}, 
     Arrowheads[{-0.05, 0.05}], Arrow[pts], 
     Text[val[[2]], Total@pts/2, {0, -1}]}, Frame -> True]]

So some play:

Manipulate[
 vis[a, {b1, b2}, 
  c], {{a, NormalDistribution, "distribution"}, {NormalDistribution, 
   CauchyDistribution, LogNormalDistribution}}, {{b1, 5, "location"}, 
  5, 7, Appearance -> "Labeled"}, {{b2, 1.5, "scale"}, 1, 2,
  Appearance -> "Labeled"}, {{c, 0.1, "fraction of peak"}, 0.1, 0.5, 
  Appearance -> "Labeled"}]

enter image description here

$\endgroup$
2
$\begingroup$

Micheal Seifert has already addressed the numerical approach. I wanted to add a different general approach that relies on the assumption that your peaks have a known shape, e.g. a Gaussian as you suggested.

You could then estimate the positions at a given fraction $n$ of the height as follows:

FullSimplify[
 PDF[NormalDistribution[mu, sigma], x] == 
  PDF[NormalDistribution[mu, sigma], mu]/n,
 Assumptions -> sigma > 0
 ]
(* Out: E^( (mu-x)^2 / (2 sigma^2) ) / n == 1 *)

PowerExpand@Log[ %[[1]] ] == Log[ %[[2]] ]
(* Out: (mu-x)^2 / (2 sigma^2) - Log[n] == 0 *)

FullSimplify[%, Assumptions -> sigma > 0]
(* (mu-x)^2 \[Equal] 2 sigma^2 Log[n] *)

pointpositions = Solve[%, x]
(* 
{
{x -> mu - Sqrt[2] sigma Sqrt[ Log[n] ]},
{x -> mu + Sqrt[2] sigma Sqrt[ Log[n] ]}
}
*)

The difference between these two locations is the width you are looking for in the case a Gaussian is a good approximation for your peak:

Differences[x /. %]
(* Out: {2 Sqrt[2] sigma Sqrt[Log[n]]} *)

In your case, with $mu=5$ and $sigma=1.5$, the 10% of the maximum height is achieved at the following points:

pointpositions /. {mu -> 5, sigma -> 1.5, n -> 10}
(* Out: { {x -> 1.78105}, {x -> 8.21895} } *)

You could use these values as starting points for FindRoot as shown in Michael Seifert's answer, or if estimates are good enough, you could calculate the width at 10% height directly:

Differences[x /. %]
(* Out: {6.4379} *)

Note: I have guided Mathematica through some of the steps, rather than giving it the full equation and letting it calculate away, because the solution it proposes to the general equation contains conditional expressions that can be more easily avoided in each step by the Assumptions given to the FullSimplify functions.

$\endgroup$

Your Answer

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

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