Re: Rookie questions about solving for small numbers

• To: mathgroup at smc.vnet.net
• Subject: [mg131105] Re: Rookie questions about solving for small numbers
• From: Ray Koopman <koopman at sfu.ca>
• Date: Tue, 11 Jun 2013 02:33:52 -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

```It's better to compute Erfc[a]-Erfc[b] as Erf[a,b],
especially when the answer is small.

When w is quadratic in x then Pr(w > 1) will always be either
Pr(x < x1) + Pr(x > x2)  or  Pr(x1 < x < x2), but will Pr(w > 1)
always have one of those forms when you include higher-order terms?
In other words, does w = 1 always have either two or no real roots?

----- Samuel Mark Young <sy81 at sussex.ac.uk> wrote:
> Thank you for the advice. What I need is something that will do that automatically (without, for example, plotting a log-log graph to find approximate solutions) - so that I can quickly run this for a large range of values of f. I also need it to work for higher order equations (i.e. past 4th order, the equation cannot be solved symbolically).
>
> Essentially, I'm solving for a taylor series type expansion, where a variable (call it w) is expanded in terms of a gaussian variable x. i.e.
> w=x+f*(x^2-s^2)+g*x^3+h*(x^4 - 3*s^4)+... (Equation 1).
>
> And then solving Pr(w>1)=e (Equation 2) for s. The range of values of x which give w>1 depends on the coefficients f, g, h, etc, as well as s (the standard deviation of the Gaussian variable x). As seen for the quadratic case, x + f*(x^2 - s^2), if f is negative we have:
> Pr(w>1)=Pr(x>x1)-Pr(x>x2)
> I've written the RHS that way for 2 reasons: it was easier to write the code that way, and this gives a sum of 2 ERFCs to solve for.
>
> The code I included in the original email works for the majority of cases - but the problems I run into are the same as those which occur for just the quadratic case:
>
> 1) If the highest order term is even, then typically for negative values of the coefficients there are no solutions to Equation 1 unless s becomes large. For example, there may be no solutions unless s>0.7, so to solve I use a temporary value for s, say stemp=0.8, and I use NSolve to find the solutions for x. I then use find root to find a value for s by solving Equation 2 - which can come out as, for example, 0.1. However, s=0.1 would have meant there were no solutions for w=1. At the moment, my best method is to use a trial and error approach until stemp matches s to within 3sf (an approximation to 3sf is fine).
>
> This doesn't take too long to do once, around 10s or so on my laptop, but I want the code to run for multiple values of f, g, h etc, so even running for only 10 values of each coefficient quickly means I have to run this thousands of times. Ideally, I'd like Mathematica to have a streamlined way of doing the trial and error, or a way for FindRoot to solve the system of equations simultaneously (so that I don't need to use a temporary value for s to solve equation 1).
>
> 2) Underflow occurs a lot - you seem to have got around that using WorkingPrecision? Where do I need to place that in the code?
>
>
> Many thanks,
> Sam
> ________________________________________
> From: Ray Koopman [koopman at sfu.ca]
> Sent: 08 June 2013 08:29
> To: Samuel Young
> Cc: mathgroup at smc.vnet.net
> Subject: RE: Re: Rookie questions about solving for small numbers
>
> Writing s instead of sigma, and letting e denote the probability
> that is variously specified as 10^-5 or 10^-20, we have
>
> Pr[x + f*(x^2 - s^2) > 1] = e,
>
> where x is Normal with mean 0 and standard deviation s. Rewrite that as
>
> Pr[z + f*s*(z^2 - 1) > 1/s] = e, where z is standard normal.
>
> If f = 0 then we have Pr[z > 1/s] = e.
>
> Plot log e as a function of log s to identify approximate solutions,
> then use FindRoot.
>
> f = 0, e = 10^-5  gives s -> .234473
> f = 0, e = 10^-20 gives s -> .107964
>
> If f > 0 then  z + f*s*(z^2 - 1) = 1/s  has two real roots, say z1 < z2.
> The coefficient of z is > 0, so we want Pr[z < z1] + Pr[z > z2].
> Again do a log-log plot, then use FindRoot.
>
> f = .5, e = 10^-5  gives s -> .173683
> f = .5, e = 10^-20 gives s -> .0792307
>
> If -1/4 < f < 0 then the roots are real for all s > 0.
> The coefficient of z^2 is < 0, so we want Pr[z1 < z < z2].
>
> f = -.125, e = 10^-5  gives s -> .271650
> f = -.125, e = 10^-20 gives s -> .126184
>
> If f < -1/4 then the roots are real only for s > sqrt[-1-4f]/(-2f).
> f = -1/2 requires s > 1.
>
> f = -1/2, e = 10^-5, WorkingPrecision->30 gives s ->
> 1.000000000213493355543893763569933088415336334927`30.
>
> f = -1/2, e = 10^-20, WorkingPrecision->60 gives s ->
> 1.0000000000000000000000000000000000000002134933555514410720359161572336783883`60.
>
> In general, s -> smin as e -> 0 with f < -1/4.
> How much precision do you need? Would an approximation suffice?
>
> ----- Samuel Mark Young <sy81 at sussex.ac.uk> wrote:
>> Yes that's correct. The code I've included does this fine for positive values of f, but can't handle negative values very well. I also run into problems when I want to solve for Pr[x + f*(x^2 - sigma^2) > 1] = 10^-20 instead.
>>
>> Sam
>> ________________________________________
>> From: Ray Koopman [koopman at sfu.ca]
>> Sent: 07 June 2013 04:19
>> To: mathgroup at smc.vnet.net
>> Subject: Re: Rookie questions about solving for small numbers
>>
>> It's still not clear exactly what you want. Is this it?
>>
>> X is a random variable whose distribution is normal with mean 0
>> and standard deviation sigma. Given a constant f, solve for sigma
>> such that Pr[x + f*(x^2 - sigma^2) > 1] = 10^-5.

```

• Prev by Date: Re: Applications and Packages, WRI Strikes Out!
• Next by Date: System Modeler and Mathematica 9
• Previous by thread: Re: Rookie questions about solving for small numbers
• Next by thread: Panel Data on RLink