[R] Finding solution for non-linear equations
peter dalgaard
pdalgd at gmail.com
Wed Oct 16 12:04:16 CEST 2013
On Oct 15, 2013, at 17:18 , Berend Hasselman wrote:
>
> On 15-10-2013, at 15:24, Ron Michael <ron_michael70 at yahoo.com> wrote:
>
>> Hi,
>>
>> I need to solve following simultaneous equations for A, B, Y1, Y2:
>>
>> B * Phi(Y1 - A) + (1-B) * Phi(Y1 + A) = 0.05
>> B * Phi(Y2 - A) + (1-B) * Phi(Y2 + A) = 0.01
>>
>> Y1 <= -1.65
>> Y2 >= -2.33
>>
>> 0 <= B <=1
>>
>> Phi is CDF for standard normal
>>
>> If there is no unique solution, then I should be able to get some feassible solution(s)
>>
>> Is there any way that using R I can achieve that?
>
>
> You cannot solve a system of 2 equations with 4 unknowns (variables).
> You can only try to find 4 values that get as close as possible (in whatever sense) to solving the system.
Think again....
A=0, B=whatever, Y1=qnorm(.05), Y2=qnorm(.01)
(4 equations in 2 unknowns is generally harder.)
There are also solutions for A+Y1 = qnorm(.05), A+Y2 = qnorm(.01), B=0, etc. --- except that they will violate the constraints on at least one of Y1, Y2 (assuming that the constraints are actually the two qnorm()s rounded to two decimal places)!
The 65536$ question is whether there are nontrivial solutions, with Y1+/-A on different sides of qnorm(.05), and consequently 0 < B < 1. I kind of suspect that there aren't any; for instance
> Y1 <- qnorm(.05)
> Y2 <- qnorm(.01)
> A <- .01
> B <- (pnorm(Y1+A)-0.05)/(pnorm(Y1+A)-pnorm(Y1-A))
> B*(pnorm(Y1-A))+(1-B)*pnorm(Y1+A)
[1] 0.05
> B*(pnorm(Y2-A))+(1-B)*pnorm(Y2+A)
[1] 0.01000091
> B*(pnorm(Y2+.001-A))+(1-B)*pnorm(Y2+.001+A)
[1] 0.01002759
Notice that we can get the first equation satisfied exactly, but then the 2nd one overshoots by an amount which is increasing in Y2. We could fix that by decreasing Y2 but that violates the constraint.
--
Peter Dalgaard, Professor
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: pd.mes at cbs.dk Priv: PDalgd at gmail.com
More information about the R-help
mailing list