Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.}. >>
share|improve this question
    
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

1 Answer 1

up vote 1 down vote accepted
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}}
share|improve this answer
    
Thanks to you and 0x4A4D. It works. –  dreichler Jun 22 '13 at 11:21

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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