0
$\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$
  • $\begingroup$ Search the website before asking your question. Here is an example mathematica.stackexchange.com/questions/10089/… . Hope it helps! $\endgroup$ – PlatoManiac Jun 19 '13 at 20:58
  • $\begingroup$ I have modified the question, now it shouldn't be a duplicate anymore $\endgroup$ – dreichler Jun 20 '13 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$ – J. M. will be back soon Jun 20 '13 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 Jun 20 '13 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$ – J. M. will be back soon Jun 20 '13 at 9:58
1
$\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}} = *) Transpose@u.u,  
 (* {d,v} = *) {Mean@#,Variance@#}&@Total[(u.Transpose@pts)^2]}]

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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