Re: Trouble with a system of equations
- To: mathgroup at smc.vnet.net
- Subject: [mg77504] Re: Trouble with a system of equations
- From: Ray Koopman <koopman at sfu.ca>
- Date: Mon, 11 Jun 2007 04:19:40 -0400 (EDT)
- References: <f4gn49$htt$1@smc.vnet.net>
On Jun 10, 4:26 am, Yaroslav Bulatov <yarosla... at 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]] ].