Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I have discrete collection of data points (10 to 10^4). I want to describe them by a continuous function and find a x value z, where areas S1 and S2 bounded by {x, f[x], yBot, yTop} are equal. Furthermore, I think it is called Matano plane.

There is 1 important special case: if function values dips beyond {yBot, yTop}, then that area should be added to the other area.

Shaded regions on this plot should be equal (S1=blue,S2=yellow):

enter image description here

Example data for this plot: n = 3; and z is the divider what I am looking for

rc = RandomChoice[{Sech[2 x - 0.5]/1.5 + 1.5, 2}]
data = Transpose[{
   # + RandomReal[] & /@ Range[-10, 30, 2],
   Tanh[#] + rc /. x -> # & /@ Range[-4, 4, 0.4]}]
f = GetRLine3[{data}, 3][x];
{xLimits, yLimits} = {
  {Min@#1, Max@#1}, 
  {Mean@#2[[1 ;; n]], Mean@#2[[-n ;; -1]]}} & @@ Transpose[data]

Where GetRLine3 is:

GetRLine3[MMStdata_, IO_: 1][x_: x] :=ListInterpolation[#, InterpolationOrder -> IO,
  Method -> "Spline"][x] & /@ (({{#[[1]]}, #[[2]]}) & /@ # & /@ MMStdata);

I can write a adhoc solution that finds z, but then it takes about 5.953s to find.

z = step = xLimits[[1]] + (xLimits[[2]] - xLimits[[1]])/2;
While[{part1, part2} = Flatten[{Integrate[f - yLimits[[1]], {x, xLimits[[1]], z}],
     -Integrate[(f - yLimits[[2]]), {x, z, xLimits[[2]]}]}] // N; part1 != part2,
 If[part1 > part2, step /= 2; z -= step, step /= 2; z += step]]

In this case z=8.075996010318093.

I can't seem to be able to use solve to find z. I have tried other regression curves

  • Fit with Table[x^k,{k,0,100}]
  • NonlinearModelFit with a/2 Erf[(x-\[Mu])/(Sqrt[2]\[Sigma])]+y0

and they are about 10x faster, but they did't work for both small and large sets of data. Finding accurate model or set of functions seems to be beyond me.

Is there a better or faster way to solve for z in range xLimits.

share|improve this question
add comment

2 Answers

up vote 6 down vote accepted

Maybe I'm missing something but isn't the value for z just the solution of the equation

(z - xLimits[[1]])*yLimits[[1]] + (xLimits[[2]] - z)*yLimits[[2]] == area

where area is the total area under the graph, i.e.

area = NIntegrate[GetRLine3[{data}, 3][x], {x, xLimits[[1]], xLimits[[2]]}][[1]]

Therefore, z is equal to

z = (xLimits[[2]]*yLimits[[2]] - xLimits[[1]]*yLimits[[1]] - area)/
      (yLimits[[2]] - yLimits[[1]])
share|improve this answer
    
This seems to be correct and it's x24 times faster. –  Margus Apr 28 '12 at 8:21
add comment

For a start it seems to me that you can use NIntegrate:

z = step = xLimits[[1]] + (xLimits[[2]] - xLimits[[1]])/2;

While[
  {part1, part2} = 
   Flatten[{
    NIntegrate[f - yLimits[[1]], {x, xLimits[[1]], z}],
   -NIntegrate[(f - yLimits[[2]]), {x, z, xLimits[[2]]}]
   }]; 
  part1 != part2, 
  If[part1 > part2, step /= 2; z -= step, step /= 2; z += step]
] // Timing

z
{0.14, Null}

10.3795
share|improve this answer
    
You are right, change form Integrate + //N to NIntegrate does give 2.5x speed up in time. –  Margus Apr 27 '12 at 21:58
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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