Re: Re: Trouble with a system of equations
- To: mathgroup at smc.vnet.net
- Subject: [mg77538] Re: [mg77504] Re: Trouble with a system of equations
- From: DrMajorBob <drmajorbob at bigfoot.com>
- Date: Tue, 12 Jun 2007 01:24:57 -0400 (EDT)
- References: <5738359.1181600553361.JavaMail.root@m35>
- Reply-to: drmajorbob at bigfoot.com
So the solution for k=3 is...
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 <koopman at sfu.ca> wrote:
> On Mon, 11 Jun 2007 07:38:00 -0500 drmajorbob at 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=
ks.
>>
>> 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 <koopman at sfu.ca> wrot=
e:
>>
>>> 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]] ].
>>>
>>>
>>>
>>
>>
>>
>> --
>> DrMajorBob at bigfoot.com
>>
>
-- =
DrMajorBob at bigfoot.com