Re: Precision/FullSimplify

*To*: mathgroup at smc.vnet.net*Subject*: [mg56285] Re: Precision/FullSimplify*From*: dh <dh at metrohm.ch>*Date*: Thu, 21 Apr 2005 05:36:03 -0400 (EDT)*References*: <d458t7$i0o$1@smc.vnet.net>*Sender*: owner-wri-mathgroup at wolfram.com

Hi Leigh, by using N[..] you are using approximate numbers. Whats wrong if you get an imaginary part that has to be considered zero if one takes into account the precision of the result. Here is a simple example that may clarify the problem. Sqrt[-10^2 3.] - 10 Sqrt[-3.] gives 0. + 3.552713678800501*^-15*I this is zero with an accuracy of approx 15 digits after the decimal separator. Considering machine precision of 16, this result seems quit reasonable. Sincerely, Daniel 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 > > 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) >