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