Take the 2-minute tour ×
Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It's 100% free, no registration required.

In order to find the minimax solution of a certain problem I have to determine the constant $c$ that provides the root for the following equation:

$$l(c)=3\times \left [ 1-\Phi \left( c-75 \right) \right]-\Phi \left(c-78 \right)=0$$

where $\Phi (.) $ is as usual the CDF of a standard normal random variable. Using Newton's algorithm then we have to iterate the following:

$$c_1=c_0-\frac{l(c_0)}{l \prime (c_0)}$$

where $l \prime (c)=-3 \times \phi \left(c-75 \right)-\phi \left(c-78 \right)$

and the corresponding R-code for say, 100 iterations:

for(i in 1:100){

l=3* (1-pnorm(c-75))-pnorm(c-78)

lprime=-3* dnorm(c-75)-dnorm(c-78)

c=c-(l/lprime)

}

c

I can see from the book that the solution should be $76.8$.

The problem is that the algorithm reaches that value only for a small neighbourhood around $c=75$ and does not converge otherwise (try it for $c=60$ or $c=80$). The algorithm is then essentially useless as in most cases I would not have knowledge of the solution and therefore I would have to pick a very rough initial estimate.

Is there anything I can do to improve the algorithm? If not, is there perhaps another algorithm that converges for a larger set of initial values?

I know this is more of a mathematical question but I would be very interested in hearing the opinions of people with experience in statistics and numerical optimization so please do not move it, at least not right away.

Thank you.

share|improve this question
add comment

1 Answer 1

up vote 6 down vote accepted

1. Getting a good starting point

Note that $\Phi(c-k)$ is increasing in $c$

Consequently consider $3\times \left [ 1-\Phi \left( c-75 \right) \right]=\Phi \left(c-78 \right)$

The LHS is a decreasing function of $c$ and the RHS is an increasing function of $c$.

The intersection point will necessarily be to the right of $3(1-\Phi(c-75))=1$; it will also be to the right of $1-\Phi(c-75)=\Phi(c-78)$ -- because the factor of 3 (which is >1) pushes the intersection further right than if it was 1.

In this case, the first of those is about 75.43 and the second is 76.5; so we should look to the right of 76.5 in this case (in fact I think it will always be the case as long as the multiplier - 3 in this case - is $>1$). So we could start at 76.5. This sort of approach should get a good starting point for a wide variety of similar cases.

2. Transforming the problem

We can do better if we transform as follows:

$\left [ 1-\Phi \left( c-75 \right) \right]=\Phi(c-78)/3$

$\Phi^{-1} \left[ 1-\Phi \left( c-75 \right) \right]=\Phi^{-1}[\Phi(c-78)/3]$

The LHS is linear (it's just $75-c$), while the right is asymptotically linear (c-78) for small values, but has a horizontal asymptote ($\Phi^{-1}(\frac{1}{3})\approx -.43$) for large values. This equality should be much more stable for a solver over a wide range of start values and we already have a pretty good start value (76.5).

3. Use a stable algorithm

Newton's method isn't stable. There are more stable options for rootfinding problems, often available in prepackaged routines. The R function uniroot uses a stable but fast algorithm*, which is quite good; you should work as far as possible with algorithms that start with upper and lower bounds on the zero and work by moving the bounds that will contain the root in a stable way. Even something basic like binary section will at least not head off to infinity, but with problems that should be reasonably near linear** regula falsi often converges a lot faster.

I think 78 (or whatever number is in its place) works as an upper limit even if the multiple of 3 increases or decreases a lot.

* see the description in the code

** as the transformed version of the problem should usually be


All of this suggests that for problems of the form $a\times \left [ 1-\Phi \left( c-l \right) \right]=\Phi \left(c-u \right)$ where $a>1$ and $u>l$, you could solve the problem

$l-c =\Phi^{-1}[\Phi(c-u)/a]\quad$ (both these function values will be negative at the intersection)

with left bound of $(l+u)/2$ and a right bound of $u$, using an algorithm that respects the bounds (perhaps regula falsi, or better something that can fall back on binary section when things don't progress).

It may well be possible to find a much better upper limit than that, but unless you're solving these problems a lot, it may not be worth the effort.

share|improve this answer
    
Thank you very much. I'll study these algorithms further. –  JohnK Apr 22 at 10:07
    
"The intersection point will necessarily be to the right of 3(1−Φ(c−75))=1;". Sorry if it is trivial but would you mind explaining that a bit? –  JohnK Apr 22 at 11:04
    
The LHS of that equation is decreasing. The RHS there is an upper bound for Φ(c−78), so where that upper bound crosses the LHS will be further left than where the actual function does (since the real crossing point will be lower down on the LHS function it must also be further right). Try drawing it! –  Glen_b Apr 22 at 11:28
    
I get it now, thanks. –  JohnK Apr 22 at 12:06
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.