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>