Find Intersection Point between Two Gaussians

(deleted)

....just realised my response was not correct!

Try this as a starter:

I don't know whether the log should be base e or base 10 - I changed it to the natural log in the code above.

When I run this with intGauss(2.5,5,1,1) I get 3.75.

Hope this helps.

Kurt

Sorry, I tried again with Find Roots - simply reversing the coefficients appears to do the trick:

Running this with intGauss(2.5,5,1,1) gives the output   W_polyRoots[0]= {cmplx(3.75,0)}.

Do check on the log base though before you use this in a real application!

Regards,

Kurt

You need to reverse the order of the coefficients in coefWave

(in the code in the original post - I didn't notice that Kurt's replies had already corrected that!)

I suggest that a completely analytic solution is possible if the baselines are the same: ln(Gaussian1)=ln(Gaussian2) produces a quadratic equation in 'x'. Why not solve it directly?

Vmmr5596 seems to be implying that the amplitudes and base lines are the same, and uses an unfortunate example with equal standard deviations. This gives a simple solution where the intersection is the average of the means, (mu1+mu2)/2 = 3.75, a degenerate double-valued root.

The quadratic method also allows for unequal Gaussian amplitudes, means, and standard deviations, but the general case may cause the roots to be complex (no real-valued intersection, easily testable), or have two real values (sometimes degenerate).

@s.r.chinn

surely FindRoots /P=coefWave does implement the analytic solution for a quadratic?

@tony:

I don't know what WaveMetrics does under the hood here. FindRoots should also have tolerance flag to ensure that precision requirements are maintained (done implicitly if the standard quadratic eqn solution with default DP variables is used).

@others:

The log() in above posts should be interpreted as IP natural logarithm ln(), since ln's of an exponential are used in the underlying derivation. I don't know the cribbed Python syntax for log().

The poster probably should have mentioned up front that the Normal function ( a special case of a Gaussian) is used.

Note that two solutions will be possible. You might want to add a search range to FindRoots if you know which single root is desired. Picture a tall skinny peak vs a laterally displaced shallow peak.The choice of computation method is entirely up to the user. My personal preference is the algebraic method.

@tony For a taste of how the quadratic formula can go wrong, see this: https://people.eecs.berkeley.edu/~wkahan/Qdrtcs.pdf

I'm checking the source to see if there's a test to use the quadratic formula; I don't think there is.

So is findroots using this surprisingly complicated formula, just in case the quadratic is ill-conditioned? That would explain why it's 10 times slower than using the high-school version of the formula to calculate real roots (with a check for degeneracy). I thought that maybe it was extra overhead of allocating all the variables and wave, or maybe that findroots really was using a numerical algorithm.

I don't know why it would be slow- a check of the source shows that the code I'm using (which isn't mine :) is using a careful implementation of the quadratic formula for degree 2. The code goes the route of finding the better-conditioned root using the quadratic formula and then using C/A to get the other root. And there is a comment, "COMPUTE DISCRIMINANT AVOIDING OVERFLOW". Clearly, this used to be FORTRAN code :)

What about using the analytical equation for the derivative of the Gaussian and finding the point(s) where that curve is disjoint in its slope? Would this be more robust and/or faster?

Can a simple algebra operation be done in FFT space?

then the speed difference has to be the cost of doing it right!

I just looked again at the code. It is originally FORTRAN that was run through f2c sometime probably in the previous millennium. The f2c-generated code appears to be grossly inefficient. I will look into doing something about that.