MathGroup Archive 2007

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

Search the Archive

Re: Re: Trouble with a system of equations

  • To: mathgroup at smc.vnet.net
  • Subject: [mg77535] Re: [mg77504] Re: Trouble with a system of equations
  • From: Ray Koopman <koopman at sfu.ca>
  • Date: Tue, 12 Jun 2007 01:23:21 -0400 (EDT)
  • Reply-to: koopman at sfu.ca

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
> 


  • Prev by Date: Re: licensing issues with V6.0 Financial Data
  • Next by Date: Re: 3D interpolation
  • Previous by thread: Re: Trouble with a system of equations
  • Next by thread: Re: Re: Trouble with a system of equations