Thank you for the replies - some of the formatting got lost via email last time. Here is a version which you should be able to read!
Hello,
I'm trying to solve a problem involving various integrals. Essentially, I'm perturbing a gaussian distribution, and then trying to find a value for sigma (the standard deviation) for which there is a 10^-5 chance of being greater than 1 (i.e. What value of sigma gives a value of 10^-5 when the pdf is integrated from 1 to infinity). The aim is to find how sigma changes with the different perturbations. The code below is a shortened version of what I'm currently using - you may not need to know all the relevant details.
Explanation: Here I'm adding a quadratic perturbation with coefficient f to the Gaussian variable x to make a new variable zeta. The aim here is to plot a graph (eventually) showing how sigma changes with f. Because I want to use this code to work with higher order equations (for which there is no analytic solution), I use a temporary value for sigma (sigmatemp - provided that the value of sigma found is close to sigmatemp, this works well) and solve numerically. Similarly, the easiest way to handle the integrations is just to find the relevant values of x for which zeta>1, and integrate the Gaussian distribution over those values. I then use FindRoot to find a value of sigma which satisfies the equation.
f = 0.5;
sigmatemp = 0.2; The values here have been picked arbitrarily
zeta = x + f (x^2 - sigma^2);
xCritical = x /. NSolve[1 == zeta /. sigma -> sigmatemp, x, Reals];
yCritical = xCritical/sigma; (* This calculates the critical values to integrate between *)
sigmatemp =sigma /.
FindRoot[
Sum[(1/
Sqrt[2*Pi]) (If[(Abs[D[zeta, x]] /. x -> xCritical[[n]]) >
0, If[(D[zeta, x] /. x -> xCritical[[n]]) > 0, 1, -1],
0]) (Integrate[Exp[-(y^2)/2], {y, yCritical[[n]], INFINITY}]), {n,
Length[xCritical]}] The Sum[] command generates a series of ERFC's (* this code will throw up errors for certain values of f, but the full code doesn't *)
+ If[(Simplify[D[zeta, x] /. x -> xCritical[[1]]]) < 0, 1, 0] = 10^(-5), {sigma, 0.00001, 2}] (* FindRoot searches for a value of sigma between 0 and 2 (the solution should always lie in this range - though putting in zero exactly results in divide by zero errors *)
There are currently 3 problems I'm having:
1) Underflow occurs in the computation a lot for certain values of f
2) I want to be able to solve for when the integrals equal 10^-20 (as oppose to 10^-5) - which is greater than machine precision. I've tried fiddling with settings like AccuracyGoal, PrecisionGoal and WorkingPrecision but can't find anything that makes it work reliably (instead of a smooth curve, I end up with jagged spikes). Its entirely possible, even likely, that I'm missing something obvious.
3) For negative values of f, there are no solutions to x + f (x^2 - sigmatemp^2)=1 if sigmatemp is too small. i.e. x-2(x^2-0.1^2)=1 has no solutions, but x-(x^2-1^2)=1 does.
The problem then, is that, when sigmatemp is increased, unless there is remarkable fine tuning, the final value of sigma found is not similar to sigmatemp. The only way I've found to handle this is to very slowly increment sigmatemp until there are real solutions, then use FindRoot to find a value for sigma, and compare sigma to sigmatemp to see if they match (to 4 s.f. is fine). However, this takes a very long time when I want to do this repeatedly for multiple values of f.
Many thanks in advance for taking the time to read this, and any help is very well appreciated. I think I have included all the information needed, but please ask if you need more. Please feel free to contact me directly at sy81 at sussex.ac.uk with any questions.
Regards,
Sam Young