       Precision/FullSimplify

• To: mathgroup at smc.vnet.net
• Subject: [mg56230] Precision/FullSimplify
• From: leigh pascoe <leigh at cephb.fr>
• Date: Wed, 20 Apr 2005 05:30:01 -0400 (EDT)
• Sender: owner-wri-mathgroup at wolfram.com

```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)

```