Rookie questions about solving for small numbers and others
- To: mathgroup at smc.vnet.net
- Subject: [mg131017] Rookie questions about solving for small numbers and others
- From: Samuel Mark Young <sy81 at sussex.ac.uk>
- Date: Tue, 4 Jun 2013 01:58:23 -0400 (EDT)
- Delivered-to: l-mathgroup@mail-archive0.wolfram.com
- Delivered-to: l-mathgroup@wolfram.com
- Delivered-to: mathgroup-outx@smc.vnet.net
- Delivered-to: mathgroup-newsendx@smc.vnet.net
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 (=1B$B&R=1B(Btemp - provided that the value of sigma found is close to =1B$B&R=1B(Btemp, this works well) and solve numerically. Similarly, the easiest way to handle the in tegrations 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;
=1B$B&R=1B(Btemp = 0.2; The values here have been picked arbitrarily
zeta = x + f (x^2 - =1B$B&R=1B(B^2);
xCritical = x /. NSolve[1 == zeta /. =1B$B&R=1B(B -> =1B$B&R=1B(Btemp, x, Reals];
yCritical = xCritical/=1B$B&R=1B(B; This calculates the critical values to integrate between
=1B$B&R=1B(Btemp = =1B$B&R=1B(B /.
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]], =1B$B!g=1B(B}]), {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), {=1B$B&R=1B(B, 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 oppos=
e 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 mis=
sing something obvious.
3) For negative values of f, there are no solutions to x + f (x^2 - =1B$B&=
R=1B(Btemp^2)=1 if =1B$B&R=1B(Btemp is too small. The problem then, is th=
at, when =1B$B&R=1B(Btemp is increased, unless there is remarkable fine tun=
ing, the final value of sigma found is not similar to =1B$B&R=1B(Btemp. The=
only way I've found to handle this is to very slowly increment =1B$B&R=1B(=
Btemp until there are real solutions, then use FindRoot to find a value for=
sigma, and compare sigma to =1B$B&R=1B(Btemp to see if they match (to 4 s.=
f. is fine). However, this takes a very long time when I want to do this re=
peatedly.
Many thanks in advance for taking the time to read this, and any help is ve=
ry well appreciated. I think I have included all the information needed, bu=
t please ask if you need more. Please feel free to contact me directly at s=
y81 at sussex.ac.uk with any questions.
Regards,
Sam Young
- Follow-Ups:
- Re: Rookie questions about solving for small numbers and
- From: Samuel Mark Young <sy81@sussex.ac.uk>
- Re: Rookie questions about solving for small numbers and others
- From: Peter Klamser <klamser@googlemail.com>
- Re: Rookie questions about solving for small numbers and