How do I make a better algorithm to fit my data?

Hey, I'm trying to determine the origin path of two variables by fitting an ellipse and running optimization in order to see if I can get a good estimate of the center.

This is the current code that formats my data in a readable viewform then tries to fit the equation of an ellipse:

clc; close all; clear

% load data

load("JR_FAIS_default.mat")

load("JR_HC_default.mat")

load("JR_FAIS_weak50.mat")

load("JR_HC_weak50.mat")

load("JR_FAIS_weak25.mat")

load("JR_HC_weak25.mat")

load("JR_FAIS_strong25.mat")

load("JR_FAIS_strong50.mat")

load("JR_HC_strong25.mat")

load("JR_HC_strong50.mat")

% combine the groups

full_group_fr_default = [fr_FAIS_default; fr_HC_default];

full_group_fr_w25 = [fr_FAIS_weak25;fr_HC_weak25];

full_group_fr_w50 = [fr_FAIS_weak50;fr_HC_weak50];

full_group_fr_s25 = [fr_FAIS_strong25;fr_HC_strong25];

full_group_fr_s50 = [fr_FAIS_strong50;fr_HC_strong50];

full_group_fx_default = [fx_FAIS_default; fx_HC_default];

full_group_fx_w25 = [fx_FAIS_weak25;fx_HC_weak25];

full_group_fx_w50 = [fx_FAIS_weak50;fx_HC_weak50];

full_group_fx_s25 = [fx_FAIS_strong25;fx_HC_strong25];

full_group_fx_s50 = [fx_FAIS_strong50;fx_HC_strong50];

full_group_fy_default = [fy_FAIS_default; fy_HC_default];

full_group_fy_w25 = [fy_FAIS_weak25;fy_HC_weak25];

full_group_fy_w50 = [fy_FAIS_weak50;fy_HC_weak50];

full_group_fy_s25 = [fy_FAIS_strong25;fy_HC_strong25];

full_group_fy_s50 = [fy_FAIS_strong50;fy_HC_strong50];

full_group_fz_default = [fz_FAIS_default; fz_HC_default];

full_group_fz_w25 = [fz_FAIS_weak25;fz_HC_weak25];

full_group_fz_w50 = [fz_FAIS_weak50;fz_HC_weak50];

full_group_fz_s25 = [fz_FAIS_strong25;fz_HC_strong25];

full_group_fz_s50 = [fz_FAIS_strong50;fz_HC_strong50];

% calculating for weak 50

% calculating the angles + resultant vector

Fres_w50 = sqrt(sum(full_group_fx_w50.^2) + sum(full_group_fy_w50.^2) + sum(full_group_fz_w50.^2));

phi_full_group_w50 = (atan2(full_group_fy_w50,full_group_fx_w50));

alpha_full_group_w50 = (atan2(full_group_fz_w50,sqrt(full_group_fx_w50.^2 + full_group_fy_w50.^2)));

theta_full_group_w50 = atan2(full_group_fz_w50, sqrt(full_group_fx_w50.^2 + full_group_fy_w50.^2));

% calculating for weak 25

% calculating the angles + resultant vector

phi_full_group_w25 = atan2(full_group_fy_w25,full_group_fx_w25);

alpha_full_group_w25 = atan2(full_group_fz_w25,sqrt(full_group_fx_w25.^2 + full_group_fy_w25.^2));

theta_full_group_w25 = atan2(full_group_fz_w25, sqrt(full_group_fx_w25.^2 + full_group_fy_w25.^2));

% calculating for default

% calculating the angles + resultant vector

phi_full_group_default = atan2(full_group_fy_default,full_group_fx_default);

alpha_full_group_default = atan2(full_group_fz_default,sqrt(full_group_fx_default.^2 + full_group_fy_default.^2));

theta_full_group_default = atan2(full_group_fz_default, sqrt(full_group_fx_default.^2 + full_group_fy_default.^2));

% calculating for strong 25

% calculating the angles + resultant vector

phi_full_group_s25 = atan2(full_group_fy_s25,full_group_fx_s25);

alpha_full_group_s25 = atan2(full_group_fz_s25,sqrt(full_group_fx_s25.^2 + full_group_fy_s25.^2));

theta_full_group_s25 = atan2(full_group_fz_s25, sqrt(full_group_fx_s25.^2 + full_group_fy_s25.^2));

% calculating for strong 50

% calculating the angles + resultant vector

phi_full_group_s50 = atan2(full_group_fy_s50,full_group_fx_s50);

alpha_full_group_s50 = atan2(full_group_fz_s50,sqrt(full_group_fx_s50.^2 + full_group_fy_s50.^2));

theta_full_group_s50 = atan2(full_group_fz_s50, sqrt(full_group_fx_s50.^2 + full_group_fy_s50.^2));

for i = 1:28

x1_s50(i,:) = rad2deg(phi_full_group_s50(i,:))-90;

y1_s50(i,:) = rad2deg(alpha_full_group_s50(i,:));

%%scale/normalize for optimization to not panic

x1_scaled_s50(i,:) = (x1_s50(i,:) - mean(x1_s50(i,:))) / max(abs(x1_s50(i,:) - mean(x1_s50(i,:))));

y1_scaled_s50(i,:) = (y1_s50(i,:) - mean(y1_s50(i,:))) / max(abs(y1_s50(i,:) - mean(y1_s50(i,:))));

major_guess(i,:) = (max(x1_scaled_s50(i,:)) - min(x1_scaled_s50(i,:)))*.5;

minor_guess(i,:) = (max(y1_scaled_s50(i,:)) - min(y1_scaled_s50(i,:)))*.5;

%%initial guess for frame of reference on LSR

initial_guess_s50(i,1:5) = [mean(x1_scaled_s50(i,:)), mean(y1_scaled_s50(i,:)),-0.2,0.3,pi/2]; % Assuming semi-major and semi-minor axes of 0.5 initially, no rotation

end

% Define the objective function (residuals) [s50]

objective_function = @(params, xy) ((((xy(:,1) - params(1)) * cos(params(5)) + (xy(:,2) - params(2)) * sin(params(5)))./params(3)).^2 + ...

(((xy(:,1) - params(1)) * sin(params(5)) - (xy(:,2) - params(2)) * cos(params(5)))./params(4)).^2 - 1);

% Define lower and upper bounds for parameters

lb = [-1, -1, -0, -0, -pi]; % Lower bounds for center_x, center_y, semi_major_axis, semi_minor_axis, rotation_angle

ub = [1, 1, 1, 1, pi]; % Upper bounds for center_x, center_y, semi_major_axis, semi_minor_axis, rotation_angle

options = optimoptions ('lsqcurvefit', 'Algorithm','trust-region-reflective', 'MaxIterations',4000, 'StepTolerance',1e-4);

for k = 1:28

% Use lsqcurvefit for optimization [s50]

optimized_params_s50(k,:) = lsqcurvefit(@(params, xy) objective_function(params, xy), initial_guess_s50(k,:), [x1_scaled_s50(k,:)', y1_scaled_s50(k,:)'], zeros(size(x1_scaled_s50(k,:)')), lb, ub, options);

end

for h = 1:28

% Extract ellipse parameters (center coordinates, semi-major axis, semi-minor axis, and rotation angle)

center_x_s50(h,:) = optimized_params_s50(h,1) * max(abs(x1_s50(h,:) - mean(x1_s50(h,:)))) + mean(x1_s50(h,:));

center_y_s50(h,:) = optimized_params_s50(h,2) * max(abs(y1_s50(h,:) - mean(y1_s50(h,:)))) + mean(y1_s50(h,:));

semi_major_axis_s50(h,:) = optimized_params_s50(h,3) * max(abs(x1_s50(h,:) - mean(x1_s50(h,:))));

semi_minor_axis_s50(h,:) = optimized_params_s50(h,4) * max(abs(y1_s50(h,:) - mean(y1_s50(h,:))));

rotation_angle_s50(h,:) = optimized_params_s50(h,5);

distance_s50(h,:) = sqrt((x1_s50(h,:)-center_x_s50(h,:)).^2 + (y1_s50(h,:)-center_y_s50(h,:)).^2);

end

for i = 1:28

theta = linspace(0,2*pi,101);

x(i,:) = center_x_s50(i,:) + semi_major_axis_s50(i,:) * cos(theta) - semi_minor_axis_s50(i,:) * sin(theta) * sin(rotation_angle_s50(i,:));

y(i,:) = center_y_s50(i,:) + semi_major_axis_s50(i,:) * cos(theta) + semi_minor_axis_s50(i,:) * sin(theta) * cos(rotation_angle_s50(i,:));

end

However, when I plot one against another, there is not a great estimate that is formed for example:

How should I go about optimizing this so that at least the rotation and fit is better?

[Mathieu NOE]

hello

you need to provide your data (mat files) if you want someone to help you

there are several submissions available on the Fex about fitting an elipse :

File Exchange - MATLAB Central (mathworks.com)

[Mathieu NOE]

hello again

problem solved ?