Trouble with a system of equations
Yaroslav Bulatov
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}]]]
Andrzej Kozlowski
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
Ray Koopman
> 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),
> 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}]]]
Your're making the problem more complicated that it needs to be.
A saturated model for k dichotomous predictors has 2^k cells and
2^k parameters. The parameter estimates are given by
LinearSolve[x,Log[p/(1-p)],
where x is the design matrix, Dimensions[x] = {2^k,2^k},
and p is the vector of observed proportions.
For a closed-form solution you need to prespecify the dummy coding,
the cell order, and the parameter order. If the dummy coding is {0,1},
the cell order is as given by Tuples[{0,1},k],
and the parameter order is as given by Subsets[Range[k]],
then the closed-form solution is
m.Log[p/(1-p)],
where
m = Inverse[ Subsets[Times@@#]& /@ Tuples[{0,1},k]] ].
dh
HI Yaroslav,
your equation in terms of a,b,c,d is linear and easy to solve. However,
by transformations like a->a -> t0/(1 + t0) you are creating rational
equations with poles. Now, as far as I can see, there is nothing that
prevents you to use a,b,c... as variables. Why do you want to make your
problem harder than necessary?
hope this helps, Daniel
dimitris
5.2
In[21]:=
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)}
Reduce[eqns, {t0, t1, t2, t3}]
Out[21]=
{t0/(1 + t0) + (t0*t1)/(1 + t0*t1) + (t0*t2)/(1 + t0*t2) +
(t0*t1*t2*t3)/(1 + t0*t1*t2*t3) == 4*m0,
(t0*t1)/(1 + t0*t1) + (t0*t1*t2*t3)/(1 + t0*t1*t2*t3) == 4*m1,
(t0*t2)/(1 + t0*t2) + (t0*t1*t2*t3)/(1 + t0*t1*t2*t3) == 4*m2,
(t0*t1*t2*t3)/(1 + t0*t1*t2*t3) == 4*m3}
Out[22]=
(m3 == 0 && m2 == 0 && m1 == 0 && m0 == 0 && t0 == 0) || (m3 == 0 &&
m1 == 0 && -1 + 4*m0 - 4*m2 != 0 &&
t0 == -((4*(m0 - m2))/(-1 + 4*m0 - 4*m2)) && t1 == 0 && (m0 -
m2)*(-1 + 4*m2) != 0 &&
t2 == (-m2 + 4*m0*m2 - 4*m2^2)/((m0 - m2)*(-1 + 4*m2)) && m2 != 0)
|| (m3 == 0 && m2 == 0 && -1 + 4*m0 - 4*m1 != 0 &&
t0 == -((4*(m0 - m1))/(-1 + 4*m0 - 4*m1)) && (m0 - m1)*(-1 + 4*m1) !
= 0 &&
t1 == (-m1 + 4*m0*m1 - 4*m1^2)/((m0 - m1)*(-1 + 4*m1)) && t2 == 0
&& m0 != 0) ||
(-1 + 4*m0 - 4*m1 - 4*m2 + 4*m3 != 0 && t0 == -((4*(m0 - m1 - m2 +
m3))/(-1 + 4*m0 - 4*m1 - 4*m2 + 4*m3)) &&
(-1 + 4*m1 - 4*m3)*(m0 - m1 - m2 + m3) != 0 &&
t1 == (m1 - 4*m0*m1 + 4*m1^2 + 4*m1*m2 - m3 + 4*m0*m3 - 8*m1*m3 -
4*m2*m3 + 4*m3^2)/
((-1 + 4*m1 - 4*m3)*(-m0 + m1 + m2 - m3)) && (m0 - m2)*(-1 + 4*m2
- 4*m3) != 0 &&
t2 == (-m2 + 4*m0*m2 - 4*m1*m2 - 4*m2^2 + m3 - 4*m0*m3 + 4*m1*m3 +
8*m2*m3 - 4*m3^2 - m2*t1 + 4*m1*m2*t1 + m3*t1 -
4*m1*m3*t1 - 4*m2*m3*t1 + 4*m3^2*t1)/((m0 - m2)*(-1 + 4*m2 -
4*m3)) && (m1 - m3)*(m2 - m3)*(-1 + 4*m3) != 0 &&
t3 == ((-m3)*t0 + 4*m1*m3*t0 + 4*m2*m3*t0 - 16*m1*m2*m3*t0 -
8*m3^2*t0 + 16*m1*m3^2*t0 + 16*m2*m3^2*t0 - 16*m3^3*t0)/
(4*(m1 - m3)*(m2 - m3)*(-1 + 4*m3))) || (m0 == m2 && 1 + 4*m1 -
4*m3 != 0 && t0 == -((4*(m1 - m3))/(1 + 4*m1 - 4*m3)) &&
-1 + 4*m1 - 4*m3 != 0 && t1 == (1 + 4*m1 - 4*m3)/(-1 + 4*m1 - 4*m3)
&& (-1 + 4*m2 - 4*m3)*(m1 - m3) != 0 &&
t2 == (m2 + 4*m1*m2 - m3 - 4*m1*m3 - 4*m2*m3 + 4*m3^2)/((-1 + 4*m2
- 4*m3)*(m1 - m3)) && (m2 - m3)*(-1 + 4*m3) != 0 &&
t3 == (m3 - 4*m2*m3 + 4*m3^2 + 2*m3*t0 - 8*m2*m3*t0 + 8*m3^2*t0)/
((m2 - m3)*(-1 + 4*m3))) ||
(m3 == 0 && m2 == 0 && m0 == 0 && 1 + 4*m1 != 0 && t0 == -((4*m1)/(1
+ 4*m1)) && -1 + 4*m1 != 0 &&
t1 == (1 + 4*m1)/(-1 + 4*m1) && t2 == 0 && m1 != 0)
Dimitris
DrMajorBob
*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)}
Here's another way of stating the same definition:
eq = {a + b + c + d == 4*m0, b + d == 4*m1, c + d == 4*m2,
d == 4*m3};
rules = {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)};
eqns = eq /. rules
{t0/(1 + t0) + (t0 t1)/(1 + t0 t1) + (t0 t2)/(1 + t0 t2) + (
t0 t1 t2 t3)/(1 + t0 t1 t2 t3) ==
4 m0, (t0 t1)/(1 + t0 t1) + (t0 t1 t2 t3)/(1 + t0 t1 t2 t3) ==
4 m1, (t0 t2)/(1 + t0 t2) + (t0 t1 t2 t3)/(1 + t0 t1 t2 t3) ==
4 m2, (t0 t1 t2 t3)/(1 + t0 t1 t2 t3) == 4 m3}
The result is something Solve can't solve for the t variables (for
whatever reason). But you didn't have to complicate things that way.
Instead, you can solve the problems separately:
s1 = Solve[eq, {a, b, c, d}]
{{a -> 4 (m0 - m1 - m2 + m3), b -> 4 (m1 - m3), c -> 4 (m2 - m3),
d -> 4 m3}}
s2 = Solve[rules /. Rule -> Equal, {t0, t1, t2, t3}]
{{t3 -> (a (1 - b - c + b c) d)/((-1 + a) b c (-1 + d)),
t1 -> (-b + a b)/(a (-1 + b)), t2 -> ((-1 + a) c)/(a (-1 + c)),
t0 -> -a/(-1 + a)}}
and combine the two solutions:
s2[[1]] /. s1
{{t3 -> ((1 - 4 (m1 - m3) - 4 (m2 - m3) +
16 (m1 - m3) (m2 - m3)) m3 (m0 - m1 - m2 + m3))/((m1 - m3) (m2 -
m3) (-1 + 4 m3) (-1 + 4 (m0 - m1 - m2 + m3))),
t1 -> (-4 (m1 - m3) + 16 (m1 - m3) (m0 - m1 - m2 + m3))/(
4 (-1 + 4 (m1 - m3)) (m0 - m1 - m2 + m3)),
t2 -> ((m2 - m3) (-1 + 4 (m0 - m1 - m2 + m3)))/((-1 +
4 (m2 - m3)) (m0 - m1 - m2 + m3)),
t0 -> -(4 (m0 - m1 - m2 + m3))/(-1 + 4 (m0 - m1 - m2 + m3))}}
eqns /. % // Simplify
{{True, True, True, True}}
Always take advantage of what you know about the problem structure.
Bobby
On Sun, 10 Jun 2007 06:20:40 -0500, Yaroslav Bulatov
<yaros...@gmail.com> 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)}
Ray Koopman
> There's an unmatched bracket in
>
>> m = Inverse[ Subsets[Times@@#]& /@ Tuples[{0,1},k]] ].
>
> and I haven't found a way (so far) to correct it so that the code works.
>
> Bobby
Ah, the joys and perils of posting code without testing it first.
Maybe someday I'll get Mathematica for my machine at home.
Aside from the extra ], that must must have been teleported from
LinearSolve[x,Log[p/(1-p)], where a ] is missing, the problem is
that I simply copied the form of an example in the Subsets online
documentation, with h changed to Times, not realizing that Times@@#
would be evaluated before Subsets got to it.
Here are two ways to get x:
ReleaseHold@Subsets[Hold@Times@@#]& /@ Tuples[{0,1},k]
or (preferably, I think)
Times@@@Subsets@#& /@ Tuples[{0,1},k].
With[{k = 2}, Inverse[ Times@@@Subsets@#& /@ Tuples[{0,1},k] ]]
{{ 1, 0, 0, 0},
{-1, 0, 1, 0},
{-1, 1, 0, 0},
{ 1,-1,-1, 1}}
With[{k = 3}, Inverse[ Times@@@Subsets@#& /@ Tuples[{0,1},k] ]]
{{ 1, 0, 0, 0, 0, 0, 0, 0},
{-1, 0, 0, 0, 1, 0, 0, 0},
{-1, 0, 1, 0, 0, 0, 0, 0},
{-1, 1, 0, 0, 0, 0, 0, 0},
{ 1, 0,-1, 0,-1, 0, 1, 0},
{ 1,-1, 0, 0,-1, 1, 0, 0},
{ 1,-1,-1, 1, 0, 0, 0, 0},
{-1, 1, 1,-1, 1,-1,-1, 1}}
>
> On Mon, 11 Jun 2007 03:19:40 -0500, Ray Koopman <koo...@sfu.ca> wrote:
>
>> On Jun 10, 4:26 am, Yaroslav Bulatov <yarosla...@gmail.com> 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}]]]
>>
>> Your're making the problem more complicated that it needs to be.
>> A saturated model for k dichotomous predictors has 2^k cells and
>> 2^k parameters. The parameter estimates are given by
>>
>> LinearSolve[x,Log[p/(1-p)],
>>
>> where x is the design matrix, Dimensions[x] = {2^k,2^k},
>> and p is the vector of observed proportions.
>>
>> For a closed-form solution you need to prespecify the dummy coding,
>> the cell order, and the parameter order. If the dummy coding is {0,1},
>> the cell order is as given by Tuples[{0,1},k],
>> and the parameter order is as given by Subsets[Range[k]],
>> then the closed-form solution is
>>
>> m.Log[p/(1-p)],
>>
>> where
>>
>> m = Inverse[ Subsets[Times@@#]& /@ Tuples[{0,1},k]] ].
>>
>>
>>
>
>
>
> --
> DrMaj...@bigfoot.com
>
DrMajorBob
Clear[mi, obs, p, lp]
mi[k_] :=
mi[k] = Inverse@Apply[Times, Subsets /@ Tuples[{0, 1}, k], {2}];
obs[k_] := obs[k] = Array[p, {2^k}]
lp[k_] := lp[k] = Log[obs[k]/(1 - obs[k])]
mi[3].lp[3]
{Log[p[1]/(1 - p[1])], -Log[p[1]/(1 - p[1])] +
Log[p[5]/(1 - p[5])], -Log[p[1]/(1 - p[1])] +
Log[p[3]/(1 - p[3])], -Log[p[1]/(1 - p[1])] + Log[p[2]/(1 - p[2])],
Log[p[1]/(1 - p[1])] - Log[p[3]/(1 - p[3])] - Log[p[5]/(1 - p[5])] +
Log[p[7]/(1 - p[7])],
Log[p[1]/(1 - p[1])] - Log[p[2]/(1 - p[2])] - Log[p[5]/(1 - p[5])] +
Log[p[6]/(1 - p[6])],
Log[p[1]/(1 - p[1])] - Log[p[2]/(1 - p[2])] - Log[p[3]/(1 - p[3])] +
Log[p[4]/(1 - p[4])], -Log[p[1]/(1 - p[1])] + Log[p[2]/(1 - p[2])] +
Log[p[3]/(1 - p[3])] - Log[p[4]/(1 - p[4])] +
Log[p[5]/(1 - p[5])] - Log[p[6]/(1 - p[6])] - Log[p[7]/(1 - p[7])] +
Log[p[8]/(1 - p[8])]}
?
Bobby
On Mon, 11 Jun 2007 16:45:38 -0500, Ray Koopman <koo...@sfu.ca> wrote:
> On Mon, 11 Jun 2007 07:38:00 -0500 drmaj...@bigfoot.com wrote:
>
>> There's an unmatched bracket in
>>
>>> m = Inverse[ Subsets[Times@@#]& /@ Tuples[{0,1},k]] ].
>>
>> and I haven't found a way (so far) to correct it so that the code wor=
>> On Mon, 11 Jun 2007 03:19:40 -0500, Ray Koopman <koo...@sfu.ca> wrot=
e:
>>
>>> On Jun 10, 4:26 am, Yaroslav Bulatov <yarosla...@gmail.com> 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}]]]
>>>
>>> Your're making the problem more complicated that it needs to be.
>>> A saturated model for k dichotomous predictors has 2^k cells and
>>> 2^k parameters. The parameter estimates are given by
>>>
>>> LinearSolve[x,Log[p/(1-p)],
>>>
>>> where x is the design matrix, Dimensions[x] = {2^k,2^k},
>>> and p is the vector of observed proportions.
>>>
>>> For a closed-form solution you need to prespecify the dummy coding,
>>> the cell order, and the parameter order. If the dummy coding is {0,1=
},
>>> the cell order is as given by Tuples[{0,1},k],
>>> and the parameter order is as given by Subsets[Range[k]],
>>> then the closed-form solution is
>>>
>>> m.Log[p/(1-p)],
>>>
>>> where
>>>
>>> m = Inverse[ Subsets[Times@@#]& /@ Tuples[{0,1},k]] ].
>>>
>>>
>>>
>>
>>
>>
>> --
>> DrMaj...@bigfoot.com
>>
>
-- =
Ray Koopman
>>> and I haven't found a way (so far) to correct it so that the code works.
>>>>> 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}]]]
>>>>
>>>> Your're making the problem more complicated that it needs to be.
>>>> A saturated model for k dichotomous predictors has 2^k cells and
>>>> 2^k parameters. The parameter estimates are given by
>>>>
>>>> LinearSolve[x,Log[p/(1-p)],
>>>>
>>>> where x is the design matrix, Dimensions[x] = {2^k,2^k},
>>>> and p is the vector of observed proportions.
>>>>
>>>> For a closed-form solution you need to prespecify the dummy coding,
>>>> the cell order, and the parameter order. If the dummy coding is {0,1},
>>>> the cell order is as given by Tuples[{0,1},k],
>>>> and the parameter order is as given by Subsets[Range[k]],
>>>> then the closed-form solution is
>>>>
>>>> m.Log[p/(1-p)],
>>>>
>>>> where
>>>>
>>>> m = Inverse[ Subsets[Times@@#]& /@ Tuples[{0,1},k]] ].
>>>>
>>>>
>>>>
>>>
>>>
>>>
>>> --
>>> DrMaj...@bigfoot.com
>>>
>>
>
>
>
> --
Yaroslav Bulatov
OK that works, thanks. One problem with the original system is that I
set up the estimation equations in terms of E[x_i y] (assuming uniform
distribution over x's), and not empirical odds ratios which apparently
make it much harder to solve (Andrzei's solution works fine for 2
variables, but takes too long for 3)
Andrzej Kozlowski
On 13 Jun 2007, at 20:42, Yaroslav Bulatov wrote:
> On Jun 11, 10:31 pm, Ray Koopman <koop...@sfu.ca> wrote:
> OK that works, thanks. One problem with the original system is that I
> set up the estimation equations in terms of E[x_i y] (assuming uniform
> distribution over x's), and not empirical odds ratios which apparently
> make it much harder to solve (Andrzei's solution works fine for 2
> variables, but takes too long for 3)
>
>
I think a variant of my method also works for three variables. I will
soon sent you a notebook, after adding to it sme comments.
Andrzej Kozlowski
Yaroslav Bulatov
becomes linear with a proper choice of variables. But now getting the
original variables is problematic.
In particular, if I replace all terms of the form Exp[a+b+c]/1+Exp[a+b
+c..] with single variables, I will have to solve the following to get
a,b,c..back
mapping[k_] := Module[{mm, params, tuples},
mm = KroneckerProduct @@ Table[{{1, 1}, {1, 0}}, {i, 1, k}];
params = (Subscript[t, #1] & ) /@ Range[1, 2^k];
tuples = (Plus @@ (params*#1) & ) /@ mm /. a_Plus :> Times @@ a;
terms = (#1/(#1 + 1) & ) /@ tuples;
{MapIndexed[#1 == Subscript[m, #2[[1]]] & , terms], params}]
Inverting it for 2 variables works
Solve @@ mapping[2]
But for 3 variables it takes too long
Solve @@ mapping[3]
Making equations polynomial and solving them gives solution in terms
of "InverseFunction" for 3 variables and takes too long for 4
Carl Woll
The mapping[k] equations are trivial to solve if one solves them for
t_1, t_2, etc in order. For example:
tsol = MapThread[Solve[#1, #2] &, {Reverse[mapping[3][[1]]],
mapping[3][[2]]}] //Flatten
tsol solves for the t_i in terms of m_i and lower order t_i. Next, we
replace the t_i in the right hand side with their known values:
solution = Thread[Rule[tsol[[All, 1]], tsol[[All, 2]] //. tsol]];
Finally, let's check that the solution works:
In[41]:= mapping[3][[1]] /. solution // Simplify
Out[41]= {True,True,True,True,True,True,True,True}
I tried this procedure for mapping[7], and it took about 12 seconds.
Carl Woll
Wolfram Research
Yaroslav Bulatov
equations in question, here's a step by step summary
http://www.yaroslavvb.com/research/reports/logistic-equations/logistic.nb
Basically the default Solve method doesn't detect that the equations
can be solved by eliminating variables in the right order, but if we
specify elimination order manually, it has no problem.
On Jun 16, 12:40 am, Carl Woll <c...@wolfram.com> wrote:
> Yaroslav Bulatov wrote:
> >As Ray Koopman and some others pointed out, the system ofequations
> >becomes linear with a proper choice of variables. But now getting the
> >original variables is problematic.
>
> >In particular, if I replace all terms of the form Exp[a+b+c]/1+Exp[a+b
> >+c..] with single variables, I will have to solve the following to get
> >a,b,c..back
>
> >mapping[k_] := Module[{mm, params, tuples},
> > mm = KroneckerProduct @@ Table[{{1, 1}, {1, 0}}, {i, 1, k}];
> > params = (Subscript[t, #1] & ) /@ Range[1, 2^k];
> > tuples = (Plus @@ (params*#1) & ) /@ mm /. a_Plus :> Times @@ a;
> > terms = (#1/(#1 + 1) & ) /@ tuples;
> > {MapIndexed[#1 == Subscript[m, #2[[1]]] & , terms], params}]
>
> >Inverting it for 2 variables works
> >Solve @@ mapping[2]
>
> >But for 3 variables it takes too long
> >Solve @@ mapping[3]
>
> >Makingequationspolynomial and solving them gives solution in terms
> >of "InverseFunction" for 3 variables and takes too long for 4
>
> The mapping[k]equationsare trivial to solve if one solves them for