Re: Precision/FullSimplify
- To: mathgroup at smc.vnet.net
- Subject: [mg56288] Re: [mg56230] Precision/FullSimplify
- From: leigh pascoe <leigh at cephb.fr>
- Date: Thu, 21 Apr 2005 05:36:05 -0400 (EDT)
- References: <200504200930.FAA17830@smc.vnet.net> <42665C97.2020502@wolfram.com>
- Sender: owner-wri-mathgroup at wolfram.com
Daniel Lichtblau wrote:
> leigh pascoe wrote:
>
>> Gentle beings,
>>
>> I hope the method of cutting and pasting from Mathematica shown below
>> is acceptable. If not I would appreciate instructions on how to do it
>> in a more acceptable fashion.
>> Consider the likelihood equations for two parameters x and y, based
>> on experimental data (a,b,c,d,e,f,g) that are observed numbers
>> (integers). The likelihood and derivatives are defined in Ma by:
>>
>> L:=(b+d+f)*Log[x]+(e+g)*Log[y]-(a+b)*Log[1+x]-(f+g)*Log[x+y]-(c+d+e)*Log[1+2*x+y]
>>
>> D[L,x]
>> D[L,y]
>>
>> Yielding the likelihood equations:
>>
>> eq1=(b+d+f)/x-(a+b)/(1+x)-2*(c+d+e)/(1+2*x+y)-(f+g)/(x+y)\[Equal]0
>> eq2=(e+g)/y-(c+d+e)/(1+2*x+y)-(f+g)/(x+y)\[Equal]0
>>
>> which I would like to solve for the maximum likelihood estimates of
>> the parameters and then evaluate for different data sets. Firstly
>> define the information matrix:
>>
>> m11:=-D[D[L,x],x]
>> m12:=-D[D[L,x],y]
>> m21:=-D[D[L,y],x]
>> m22:=-D[D[L,y],y]
>> m:={{m11,m12},{m21,m22}}//MatrixForm
>> m
>>
>> which can be inverted when the equations are solved to give the
>> variances of the estimates. To solve the equations I tried
>>
>> Sols=Solve[{eq1,eq2},{x,y}]
>>
>> Which gives me the answers I want, however the expressions are quite
>> large. By inspection I find that the last solution is the one that
>> corresponds to the maximum likelihood. However I run into problems
>> evaluating the resulting expression. If I try to evaluate the last
>> solution numerically I get a complex number, I assume due to rounding
>> errors. e.g.
>>
>> N[Sols] /.{a->50,b->50,c\[Rule]2,d\[Rule]4,e->100,f\[Rule]2,g->100}
>>
>> produces
>>
>> \!\({{y \[Rule] \(-5.336646740150529`\) +
>> 1.2546145287986817`*^-14\ \[ImaginaryI],
>> x \[Rule] \(\(3.8558775093812896`\)\(\[InvisibleSpace]\)\) -
>> 2.220446049250313`*^-16\ \[ImaginaryI]}, {y \[Rule] \
>> \(-1.067199413695643`\) - 3.6306168234410265`*^-14\ \[ImaginaryI],
>> x \[Rule] \(-0.4135698170735974`\) -
>> 6.661338147750939`*^-16\ \[ImaginaryI]}, {y \[Rule]
>> \(\(50.`\)\(\
>> \[InvisibleSpace]\)\) + 1.1823875212257917`*^-14\ \[ImaginaryI],
>> x \[Rule] \(\(1.`\)\(\[InvisibleSpace]\)\) +
>> 6.661338147750939`*^-16\ \[ImaginaryI]}}\)
>>
>> The last solution is the correct one but without the addition of the
>> small complex part. I then tried
>>
>> N[Sols, $MaxPrecision]
>> /.{a->50,b->50,c\[Rule]2,d\[Rule]4,e->100,f\[Rule]2,
>> g->100}
>>
>> Which gives me presumably exact expressions, but it is still not
>> clear that they are real. So I tried
>>
>> FullSimplify[ N[Sols, $MaxPrecision]
>> /.{a->50,b->50,c\[Rule]2,d\[Rule]4,e->100,f\[Rule]2, g->100}]
>> N[%]
>>
>> {{y\[Rule]-5.33665,x\[Rule]3.85588},{y\[Rule]-1.0672,
>> x\[Rule]-0.41357},{y\[Rule]50.,x\[Rule]1.}}
>>
>> and the last solution appears to be the one I want, but it takes a
>> long time. Now let's try this with a different data set.
>>
>> FullSimplify[
>> N[Sols, $MaxPrecision] /.{a\[Rule]42,b\[Rule]90,c\[Rule]5,d\[Rule]29,
>> e\[Rule]49,f\[Rule]6,g\[Rule]17}]
>>
>> This one seems to hang forever. In any case there must be a better
>> way of evaluating the solutions for different data sets. Should I be
>> telling Mathematica to look for real solutions, using NSolve,
>> substituting in eq1 and eq2 before solving?? It is preferable to
>> generate all the solutions before choosing the one that gives the
>> maximum likelihood. Any suggestions, observations, help will be
>> greatly appreciated.
>>
>> Leigh
>> (Inexperienced User)
>
>
>
> If you first do
>
> SetOptions[Roots,Cubics->False,Quartics->False];
>
> you will get solutions in a form that is not subject to vagaries of
> cancellation of small imaginary parts.
>
> Here is a related recent discussion in this forum.
>
> http://forums.wolfram.com/mathgroup/archive/2005/Feb/msg00811.html
>
>
> Daniel Lichtblau
> Wolfram Research
>
>
Thank you that works marvellously.
Leigh
- References:
- Precision/FullSimplify
- From: leigh pascoe <leigh@cephb.fr>
- Precision/FullSimplify