What does it mean to "fit" it to you? What will you do with that fit? The data is not a function, just a set of points, that happen to represent a closed curve in the plane.
A simple solution might be to convert to polar coordinates.
Now, fit the function r(theta). Be careful, since r is a periodic function. It should have the property that r(0) == r(2*pi). So a polynomial fit would be a bad idea. Instead, perhaps fit this as a sum.
r(theta) = a0 + a1*sin(theta) + b1*cos(theta) + a2*sin(2*theta) + b2*cos(2*theta) + ...
Essentially, I am describing what is essentially a Fourier series there (that is, not an fft, but a Fourier seies.)
As it is constructed, the function ALWAYS has the property of being periodic. And, while you could, in theory, trivially use the curve fitting toolbox for this task, there is no reason to do so. All of the parameters can be estimated using one call to backslash.
x = x(:);
y = y(:);
mux = mean(x);
muy = mean(y);
[Theta,R] = cart2pol(x- mux,y-muy);
[Theta,ind] = sort(Theta);
R = R(ind);
n = numel(x);
M = 6;
A = [ones(n,1),sin(Theta*(1:M)),cos(Theta*(1:M))];
coeffs = A\R;
Rhat = A*coeffs;
Larger values of M will yield a better approximation. A perfect circle would have m == 0. Because your curve has an oval shape, you need to have a few extra terms, but not many. M == 6 will require the estimation of 13 coefficients for the series.
I lack your data, so I generated some random crap data as an example, using ginput.
Nothing exciting. Vaguely ovoid. I'm actually surprised it was that good.
R as a function of Theta seems reasonable.
And now, we can reconstruct the original curve.
[xhat,yhat] = pol2cart(Theta,Rhat);
plot(x,y,'bo',xhat + mux,yhat + muy,'-xr')
Surprisingly good for a simple, moderately low order Fourier series approximation, given the data came from my shaky eye/mouse coordinated hand.