Services & ResourcesWolfram ForumsMathGroup Archive

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

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)

• Follow-Ups:
  • Re: Precision/FullSimplify
    • From: yehuda ben-shimol <bsyehuda@gmail.com>