# Fit implicit function

Can anyone help me with fitting an implicit function to data?

The function is given by:

a*x^2 + 2*b*x*y + c*y^2 = e


under the condition

a*c-b^2=1


something like a least square fit with $\chi^2$ value would be really nice. Unfortunately I have no idea how to do a fit with an only implicitly given function.

I tried the solutions according to this thread How to use FindFit to fit an implicit function?, but i got the following errors:

pts = {{0.943851, 0.0255852}, {-2.23615, 0.0255852}, {2.12656, -0.0845208}, {2.32489, -0.0354676}, {-0.284106, -0.0354676}, {-2.64798, 0.0415494},{1.41402,0.0415494}, {-0.114051, 0.0669646}, {-2.99405,0.0669646}, {-3.92389,0.0825293}, {-0.273894,0.0825293}, {-2.66303, 0.0877706}, {-1.76303,0.0877706}}

fitfunc[a_, b_, d_, x_] := y /. Solve[a*x^2 + 2*b*x*y + (1 + b^2)/a*y^2 == d, {y},InverseFunctions -> True]
FindFit[pts, fitfunc[a, b, d, x], {a, b, d}, x]
FindFit::eqineq: Constraints in {(-a b x+Sqrt[a d+a b^2 d-a^2 x^2])/(1+b^2)} are not all equality or inequality constraints. With the exception of integer domain constraints for linear programming, domain constraints or constraints with Unequal (!=) are not supported. >>


I also tried:

fitfunc2[a_?NumericQ, b_?NumericQ, d_?NumericQ, , x_?NumericQ] :=  y /. FindRoot[a*x^2 + 2*b*x*y + (1 + b^2)/a*y^2 == d, {y, 1.}]
FindFit[pts, fitfunc2[a, b, d, x], {a, b, d}, x]
FindFit::nrlnum: The function value {-0.0255852+fitfunc2[1.,1.,1.,0.943851],-0.0255852+fitfunc2[1.,1.,1.,-2.23615],0.0845208\[VeryThinSpace]+fitfunc2[1.,1.,1.,2.12656],0.0354676\[VeryThinSpace]+fitfunc2[1.,1.,1.,2.32489],0.0354676\[VeryThinSpace]+fitfunc2[1.,1.,1.,-0.284106],-0.0415494+fitfunc2[1.,1.,1.,-2.64798],-0.0415494+fitfunc2[1.,1.,1.,1.41402],-0.0669646+fitfunc2[1.,1.,1.,-0.114051],-0.0669646+fitfunc2[1.,1.,1.,-2.99405],-0.0825293+fitfunc2[1.,1.,1.,-3.92389],-0.0825293+fitfunc2[1.,1.,1.,-0.273894],-0.0877706+fitfunc2[1.,1.,1.,-2.66303],-0.0877706+fitfunc2[1.,1.,1.,-1.76303]}
is not a list of real numbers with dimensions {13} at {a,b,d} = {1.,1.,1.}. >>

-
Search the website before asking your question. Here is an example mathematica.stackexchange.com/questions/10089/… . Hope it helps! – PlatoManiac Jun 19 '13 at 20:58
I have modified the question, now it shouldn't be a duplicate anymore – dreichler Jun 20 '13 at 9:34
Okay, but you still need to include pts somewhere, as the rest of us don't have it... on another note, you wouldn't happen to have initial guesses for your parameters, no? – J. M. Jun 20 '13 at 9:35
I tried a graphical fit with a manipulate i get roughly a=-0.04, b=-0.97 and d=-0.17 – dreichler Jun 20 '13 at 9:52
If you're fitting ellipses to data, you might be interested to know that there are special methods like this... but to start you off: NArgMin[{Norm[Function[{x, y}, a x^2 + 2 b x y + c y^2 - e] @@@ pts], b^2 - a c == -1}, {a, b, c, e}]. – J. M. Jun 20 '13 at 9:58

Clear[a,b,c];
NMinimize[{Variance[#.{{a,b},{b,c}}.#& /@ pts], Det@{{a,b},{b,c}} == 1}, {a,b,c}]

(* {0.0051474, {a -> -0.0429303, b -> -1.00828, c -> -46.9747}} *)


You didn't say so, but I suspect you want a > 0 and c > 0. That just complements {a,b,c}:

Clear[a,b,c];
NMinimize[{Variance[#.{{a,b},{b,c}}.#& /@ pts],
Det@{{a,b},{b,c}} == 1 && a > 0 && c > 0}, {a,b,c}]

(* {0.0051474, {a -> 0.0429309, b -> 1.00814, c -> 46.9672}} *)


This handles the constraints implicitly and is much faster. It solves for a rotation angle and for the logs of the two scale factors, making one the negative of the other.

Clear[s,t]; NMinimize[Variance@Total[
({{Cos@t,Sin@t}*E^s,{-Sin@t,Cos@t}*E^-s}.Transpose@pts)^2], {s,t}]

(* {0.0051474, {s -> -1.92503, t -> -0.0214708}} *)

With[{u = {{Cos@t,Sin@t}*E^s,{-Sin@t,Cos@t}*E^-s}/.%[[2]]},
{(* {{a,b},{b,c}} = *) Transpose@u.u,
(* {d,v} = *) {Mean@#,Variance@#}&@Total[(u.Transpose@pts)^2]}]

{{{0.0429303,1.00828},{1.00828,46.9747}},{0.187624,0.0051474}}

-
Thanks to you and 0x4A4D. It works. – dreichler Jun 22 '13 at 11:21