Mathematica 9 is now available
Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2005
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2005

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

Search the Archive

Re: Precision/FullSimplify

  • To: mathgroup at smc.vnet.net
  • Subject: [mg56290] Re: Precision/FullSimplify
  • From: Peter Pein <petsie at arcor.de>
  • Date: Thu, 21 Apr 2005 05:36:08 -0400 (EDT)
  • References: <d458t7$i0o$1@smc.vnet.net>
  • Sender: owner-wri-mathgroup at wolfram.com

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

You first approximate exact numbers which appear in Sols and _after_
that replace the variables by exact numbers. This can't be the optimum
strategy...

N[Sols]/.{a->50,...}//Chop should work here, but it is IMHO better to do
exact calculations as long as possible:

subs = {a -> 50, b -> 50, c -> 2, d -> 4, e -> 100, f -> 2, g -> 100};
sols = Solve[eq1 && eq2, {x, y}];
nsol = Simplify[RootReduce[sols /. subs]]

which will give you three real solutions:

{{y -> (1/104)*(-333 - Sqrt[49289]), x -> (1/104)*(179 + Sqrt[49289])},
 {y -> (1/104)*(-333 + Sqrt[49289]), x -> (1/104)*(179 - Sqrt[49289])},
 {y -> 50, x -> 1}}

You want the solution(s) resulting in real L:

L0=Select[L/.subs/.nsol,Im[#]==0&]

Here is your extremal value of L:

{-100 Log[2]+200 Log[50]-102 Log[51]-106 Log[53]}

which is approx.
N[L0,16]
-108.8072743447715

> 
> N[Sols, $MaxPrecision] /.{a->50,b->50,c\[Rule]2,d\[Rule]4,e->100,f\[Rule]2,
>     g->100}

$MaxPrecision is usually set to Infinity and its use is to no avail in
this context.

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


-- 
Peter Pein
Berlin


  • Prev by Date: Re: Reduce bug report
  • Next by Date: Re: Precision/FullSimplify
  • Previous by thread: Re: Precision/FullSimplify
  • Next by thread: Re: Precision/FullSimplify