1
$\begingroup$

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.}. >>
$\endgroup$
8
  • $\begingroup$ Search the website before asking your question. Here is an example mathematica.stackexchange.com/questions/10089/… . Hope it helps! $\endgroup$ Commented Jun 19, 2013 at 20:58
  • $\begingroup$ I have modified the question, now it shouldn't be a duplicate anymore $\endgroup$
    – dreichler
    Commented Jun 20, 2013 at 9:34
  • $\begingroup$ 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? $\endgroup$ Commented Jun 20, 2013 at 9:35
  • $\begingroup$ I tried a graphical fit with a manipulate i get roughly a=-0.04, b=-0.97 and d=-0.17 $\endgroup$
    – dreichler
    Commented Jun 20, 2013 at 9:52
  • $\begingroup$ 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}]. $\endgroup$ Commented Jun 20, 2013 at 9:58

1 Answer 1

3
$\begingroup$
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}} = *) [email protected],  
 (* {d,v} = *) {Mean@#,Variance@#}&@Total[(u.Transpose@pts)^2]}]

{{{0.0429303,1.00828},{1.00828,46.9747}},{0.187624,0.0051474}}
$\endgroup$
1
  • $\begingroup$ Thanks to you and 0x4A4D. It works. $\endgroup$
    – dreichler
    Commented Jun 22, 2013 at 11:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.