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

Re: Precision/FullSimplify

• To: mathgroup at smc.vnet.net
• Subject: [mg56280] Re: Precision/FullSimplify
• From: "Jens-Peer Kuska" <kuska at informatik.uni-leipzig.de>
• Date: Thu, 21 Apr 2005 05:35:59 -0400 (EDT)
• Organization: Uni Leipzig
• References: <d458t7$i0o$1@smc.vnet.net>
• Sender: owner-wri-mathgroup at wolfram.com

Hi,

what many Chop[] do ??

Regards
  Jens

"leigh pascoe" <leigh at cephb.fr> schrieb im 
Newsbeitrag news:d458t7$i0o$1 at smc.vnet.net...
> 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)
>