Re: Re: Trouble with a system of equations
- To: mathgroup at smc.vnet.net
- Subject: [mg77585] Re: [mg77504] Re: Trouble with a system of equations
- From: Ray Koopman <koopman at sfu.ca>
- Date: Wed, 13 Jun 2007 07:36:36 -0400 (EDT)
- Reply-to: koopman at sfu.ca
Yes, that's it.
On Mon, 11 Jun 2007 22:30:21 -0500 drmajorbob at bigfoot.com wrote:
> 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 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 <koopman at sfu.ca> wrote:
>>>
>>>> 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
>