Services & ResourcesWolfram ForumsMathGroup Archive

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

Re: Trouble with a system of equations

• To: mathgroup at smc.vnet.net
• Subject: [mg77500] Re: [mg77488] Trouble with a system of equations
• From: Andrzej Kozlowski <akoz at mimuw.edu.pl>
• Date: Mon, 11 Jun 2007 04:17:34 -0400 (EDT)
• References: <200706101120.HAA18117@smc.vnet.net>

On 10 Jun 2007, at 20:20, Yaroslav Bulatov wrote:

> Hi, I'm trying to solve a certain kind of system of equations, and
> while they are solvable by hand, Mathematica 6.0 has problems solving
> it
>
> Here's an example
>
> eqns = {a + b + c + d == 4*m0, b + d == 4*m1, c + d == 4*m2, d ==
> 4*m3} /. {a -> t0/(1 + t0), b -> (t0*t1)/(1 + t0*t1), c -> (t0*t2)/(1
> + t0*t2), d -> (t0*t1*t2*t3)/(1 + t0*t1*t2*t3)}
> Solve[eqns, {t0, t1, t2, t3}]
>
> The solution can be found by hand and verified below
>
> sol = {t0 -> a/(1/4 - a), t1 -> (b/(1/4 - b))*((1/4 - a)/a), t2 -> (c/
> (1/4 - c))*((1/4 - a)/a), t3 -> (m3/(1/4 - m3))*(a/(1/4 - a))*((1/4 -
> b)/b)*((1/4 - c)/c)} /. {a -> m0 - m1 - m2 + m3, b -> m1 - m3, c -> m2
> - m3}
> eqns /. sol // Simplify
>
> This is an example of estimating equations for a saturated logistic
> regression model with 2 independent variables. I'd like to see if
> formulas also exist for more variables, but they are too cumbersome to
> solve by hand. Are there any Mathematica tricks I can use to answer
> this question?
>
> Here's the procedure that generates the system of equations for d
> variables (d=2 produces the system above)
>
> logeq[d_] := Module[{bounds, monomials, params,
> partition,derivs,sums},
>    xs = (Subscript[x, #1] & ) /@ Table[i, {i, 1, d}];
>     monomials = Subsets[xs]; monomials = (Prepend[#1, 1] & ) /@
> monomials;
>     monomials = (Times @@ #1 & ) /@ monomials;
>     params = (Subscript[th, #1] & ) /@ Table[i, {i, 0, 2^d - 1}];
>     monomials = (Times @@ #1 & ) /@ Thread[{params, monomials}];
>     partition = Log[1 + Exp[Plus @@ monomials]];
>     derivs = (D[partition, Subscript[th, #1]] & ) /@
>       Table[i, {i, 0, 2^d - 1}]; bounds = ({#1, 0, 1} & ) /@ xs;
>     sums = (Table[#1, Evaluate[Sequence @@ bounds]] & ) /@ derivs;
>     sums = (Plus @@ #1 & ) /@ (Flatten[#1] & ) /@ sums;
>     Thread[sums == Table[Subscript[m, i], {i, 0, 2^d - 1}]]]
>
>


Here is a simple way to solve your original system. I have not tried  
any of the more general ones.

  eqns = {a + b + c + d == 4*m0, b + d == 4*m1, c + d == 4*m2, d ==  
4*m3} /. {a -> t0/(1 + t0), b -> (t0*t1)/(1 + t0*t1),
      c -> (t0*t2)/(1 + t0*t2), d -> (t0*t1*t2*t3)/(1 + t0*t1*t2*t3)};

  g = {t0/(t0 + 1) + (t1*t0)/(t0*t1 + 1) + (t2*t0)/(t0*t2 + 1) +
        (t1*t2*t3*t0)/(t0*t1*t2*t3 + 1) == 4*m0,
      (t0*t1)/(t0*t1 + 1) + (t0*t2*t3*t1)/(t0*t1*t2*t3 + 1) == 4*m1,
      (t0*t2)/(t0*t2 + 1) + (t0*t1*t3*t2)/(t0*t1*t2*t3 + 1) == 4*m2,
      (t0*t1*t2*t3)/(t0*t1*t2*t3 + 1) == 4*m3} /. Equal -> Subtract;

g1 = Map[Numerator[Together[#]] &, g];

gr = GroebnerBasis[g1, {t0, t1, t2, t3}];

Simplify[Solve[gr == 0, {t0, t1, t2, t3}]]
{{t0 -> -((4*(-m0 + m1 + m2 - m3))/(-4*m0 + 4*m1 + 4*m2 - 4*m3 + 1)),
    t1 -> ((-4*m0 + 4*m1 + 4*m2 - 4*m3 + 1)*(m3 - m1))/((-m0 + m1 +  
m2 - m3)*
       (-4*m1 + 4*m3 + 1)), t2 -> ((-4*m0 + 4*m1 + 4*m2 - 4*m3 + 1)* 
(m2 - m3))/
      ((4*m2 - 4*m3 - 1)*(-m0 + m1 + m2 - m3)),
    t3 -> ((4*m1 - 4*m3 - 1)*(-m0 + m1 + m2 - m3)*m3*(-4*m2 + 4*m3 +  
1))/
      ((-4*m0 + 4*m1 + 4*m2 - 4*m3 + 1)*(m1 - m3)*(m3 - m2)*(4*m3 -  
1))}}


Andrzej Kozlowski

• References:
  • Trouble with a system of equations
    • From: Yaroslav Bulatov <yaroslavvb@gmail.com>