MathGroup Archive 2005

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: A nasty 2x2 system of equations?

  • To: mathgroup at smc.vnet.net
  • Subject: [mg54144] Re: [mg54098] A nasty 2x2 system of equations?
  • From: Daniel Lichtblau <danl at wolfram.com>
  • Date: Fri, 11 Feb 2005 03:34:54 -0500 (EST)
  • References: <200502100747.CAA16654@smc.vnet.net>
  • Sender: owner-wri-mathgroup at wolfram.com

Bruce Colletti wrote:

>Re Mathematica 5.1.
>
>The code below keeps running without getting an answer.  Is the code flawed or is this a really nasty system of 2-equations in 2-unknowns?
>
>Here's background:  X and Y are independent beta-distributed random variables and Z is a convex combination of X and Y, i.e., Z = cX + (1 - c)Y.  
>
>Although Z need not be beta-distributed, let's pretend it is and in turn, solve for its parms (a3 and b3) in terms of c and the known parms of X and Y.  
>
>Thanks.
>
>Bruce
>
>-------------
>
>
>X = BetaDistribution[a1, b1]; 
>Y = BetaDistribution[a2, b2]; 
>Z = BetaDistribution[a3, b3]; 
>
>Solve[{Mean[Z] == c*Mean[X] + (1 - c)*Mean[Y], 
>   Variance[Z] == c^2*Variance[X] + (1 - c)^2*Variance[Y]}, {a3, b3}]
>  
>

I'm not sure that it should be quite so hard but certainly it's not 
pretty once you take into account the complexity of coefficients.. You 
can solve it by breaking into steps, removing denominators, finding a 
"good" generating set allowing for "generic' coefficients that can 
themselves be rational functions in the other variables. This can be 
done as below.

eqns = {a3/(a3 + b3) == (a2*(1 - c))/(a2 + b2) + (a1*c)/(a1 + b1),
 (a3*b3)/((a3 + b3)^2*(1 + a3 + b3)) ==
  (a2*b2*(1 - c)^2)/((a2 + b2)^2*(1 + a2 + b2)) +
   (a1*b1*c^2)/((a1 + b1)^2*(1 + a1 + b1))};
rats = Map[#[[1]]-#[[2]]&, eqns];
polys = Numerator[Together[rats]]

In[19]:= Timing[gb = GroebnerBasis[polys, {a3,b3},
  CoefficientDomain->RationalFunctions];]
Out[19]= {0.44 Second, Null}

This is now straightforward for Solve to handle.

In[20]:= Timing[soln = Solve[gb==0, {a3,b3}];]
Out[20]= {0.18 Second, Null}

In[21]:= InputForm[Timing[soln2 = Simplify[soln]]]

Out[21]//InputForm=
{0.92999999999999886*Second, {{a3 -> 0, b3 -> 0}, {a3 -> 0, b3 -> 0},
  {a3 -> -(((-(a2*b1*(-1 + c)) + a1*(a2 + b2*c))*
       (a2*b1*(1 + b1)*(-1 + c)*(b2 + c + a2*c) + a1^2*b2*(-1 + c)*
         (a2 + c + b2*c) - a1*(a2^2*b1*c + b2*(1 + b2)*(1 + b1 - c)*c +
          a2*(b1*c + b2*(1 - c + 2*b1*(1 - c + c^2))))))/
      (a1^3*a2*b2*(-1 + c)^2 + a2*b1^2*(1 + b1)*b2*(-1 + c)^2 +
       a1^2*a2*(1 + 3*b1)*b2*(-1 + c)^2 +
       a1*b1*(a2^3*c^2 + b2^2*(1 + b2)*c^2 + a2^2*(1 + 3*b2)*c^2 +
         a2*b2*(2 + 3*b1*(-1 + c)^2 - 4*c + (4 + 3*b2)*c^2)))),
   b3 -> (a1^3*b2^2*(-1 + c)^2*(a2 + c + b2*c) - a2*b1^2*(1 + b1)*(-1 + c)*
       (b2^2 + a2*(1 + a2)*c^2 + b2*(c + 2*a2*c)) -
      a1^2*b2*(-1 + c)*(2*a2^2*b1*c + b2*(1 + b2)*(1 + 2*b1 - c)*c +
        a2*(b1*c*(1 + c) + b2*(1 - c + b1*(3 - 2*c + 3*c^2)))) +
      a1*b1*(b2^2*(1 + b2)*(1 + b1 - c)*c + a2^3*b1*c^2 +
        a2^2*c*(b1*c + b2*(2 - 3*c + c^2 + b1*(4 - 4*c + 3*c^2))) +
        a2*b2*(c*(1 - c + b1*(2 - c + c^2)) + b2*(2 - 3*c + 2*c^2 - c^3 +
            b1*(3 - 4*c + 4*c^2)))))/(a1^3*a2*b2*(-1 + c)^2 +
      a2*b1^2*(1 + b1)*b2*(-1 + c)^2 + a1^2*a2*(1 + 3*b1)*b2*(-1 + c)^2 +
      a1*b1*(a2^3*c^2 + b2^2*(1 + b2)*c^2 + a2^2*(1 + 3*b2)*c^2 +
        a2*b2*(2 + 3*b1*(-1 + c)^2 - 4*c + (4 + 3*b2)*c^2)))}}}

Now you might want to substitue into the original equations, or take 
limits, to see if any of these cause problems with vanishing denominators.


Daniel Lichtblau
Wolfram Research




  • Prev by Date: Re: finding package in ExtraPackages`Enhancements`
  • Next by Date: Re: Re: finding package in ExtraPackages`Enhancements`
  • Previous by thread: A nasty 2x2 system of equations?
  • Next by thread: Re: A nasty 2x2 system of equations?