Fitting two curves with shared parameters to two datasets
3 views (last 30 days)
Show older comments
I want to fit the two functions
to a dataset of the form . Here are a dimensionless distance and two angles (degrees) respectively.
I generate the dataset using the above functions with some Gaussian noise on top. The data generation seems to work fine.
Yet, the call "lsqcurvefit(fun, x0, alpha_data_tab, beta_data_tab);" gives me an error which indicates that there are some dimension mismatch. However, I do not see it. Can anyone help me with this issue?
clc
close all
clear all
set(0,'defaulttextinterpreter','latex')
set(0,'defaultAxesTickLabelInterpreter','latex');
set(0,'defaultLegendInterpreter','latex');
set(0, 'DefaultLineLineWidth', 2);
set(0,'defaultAxesFontSize',15)
%True parameters
Af0_true = -0.017767; %Fast initial amplitude
As0_true = 0.017767; %Slow initial amplitude
lambda_f_true = -4.0571e-04; %Fast damping parameter; 1/cm
lambda_s_true = 6.4718e-05; %Slow damping parameter; 1/cm
angular_vel_f_true = 0.0228; %Angular velocity; rad/cm
angular_vel_s_true = 0.0050; %Angular velocity; rad/cm
dx = 0.01; %cm
s_true = 0:dx:1000; %cm
std = 0.1; %deg
spacing = 4;
s_data1 = 275:spacing:300;
s_data2 = 350:spacing:375;
s_data3 = 425:spacing:450;
s_data4 = 500:spacing:525;
s_data = [s_data1, s_data2, s_data3, s_data4];
[s, alpha, beta, s_data, alpha_data, beta_data] = dataMaker(Af0_true, As0_true, lambda_f_true, lambda_s_true, angular_vel_f_true, angular_vel_s_true, s_true, s_data, std);
alpha_data_tab = [s_data(:), alpha_data(:)]; % Concatenate Vectors
beta_data_tab = [s_data(:), beta_data(:)]; % Concatenate Vectors
fun1 = @(v,s_data)v(1)*exp(v(3)*s_data).*cos(v(5)*s_data) + v(2)*exp(v(4)*s_data).*cos(v(6)*s_data);
fun2 = @(v,s_data)v(1)*exp(v(3)*s_data).*sin(v(5)*s_data) + v(2)*exp(v(4)*s_data).*sin(v(6)*s_data);
fun = @(v,s_data)[fun1(v,s_data)', fun2(v,s_data)']; % Composite Function
%initial guess
percent = 0.2; %Percent deviation
x0 = [Af0_true-percent*Af0_true, As0_true-percent*As0_true, lambda_f_true-percent*lambda_f_true, lambda_s_true-percent*lambda_s_true,...
angular_vel_f_true-percent*angular_vel_f_true, angular_vel_s_true-percent*angular_vel_s_true];
fun(x0, s_data)
X = lsqcurvefit(fun, x0, alpha_data_tab, beta_data_tab);
function [s_total, alpha_total, beta_total, s_data, alpha_data, beta_data] = dataMaker(Af0_true, As0_true, lambda_f_true, lambda_s_true, angular_vel_f_true, angular_vel_s_true, s_total, s_data, std_deg)
s_total = s_total;
std_rad = std_deg*(pi/180); %rad
%True parameter vector
v0 = [Af0_true, As0_true, lambda_f_true, lambda_s_true, angular_vel_f_true, angular_vel_s_true];
alpha_total = v0(1)*exp(v0(3)*s_total).*cos(v0(5)*s_total) + v0(2)*exp(v0(4)*s_total).*cos(v0(6)*s_total);
beta_total = v0(1)*exp(v0(3)*s_total).*sin(v0(5)*s_total) + v0(2)*exp(v0(4)*s_total).*sin(v0(6)*s_total);
s_data = s_data;
alpha_data = v0(1)*exp(v0(3)*s_data).*cos(v0(5)*s_data) + v0(2)*exp(v0(4)*s_data).*cos(v0(6)*s_data);
beta_data = v0(1)*exp(v0(3)*s_data).*sin(v0(5)*s_data) + v0(2)*exp(v0(4)*s_data).*sin(v0(6)*s_data);
alpha_noise = std_rad*randn(1,length(s_data));
beta_noise = std_rad*randn(1,length(s_data));
alpha_data = alpha_data + alpha_noise;
beta_data = beta_data + beta_noise;
end
0 Comments
Accepted Answer
Torsten
on 13 Nov 2023
Edited: Torsten
on 13 Nov 2023
clc
close all
clear all
set(0,'defaulttextinterpreter','latex')
set(0,'defaultAxesTickLabelInterpreter','latex');
set(0,'defaultLegendInterpreter','latex');
set(0, 'DefaultLineLineWidth', 2);
set(0,'defaultAxesFontSize',15)
%True parameters
Af0_true = -0.017767; %Fast initial amplitude
As0_true = 0.017767; %Slow initial amplitude
lambda_f_true = -4.0571e-04; %Fast damping parameter; 1/cm
lambda_s_true = 6.4718e-05; %Slow damping parameter; 1/cm
angular_vel_f_true = 0.0228; %Angular velocity; rad/cm
angular_vel_s_true = 0.0050; %Angular velocity; rad/cm
dx = 0.01; %cm
s_true = 0:dx:1000; %cm
std = 0.1; %deg
spacing = 4;
s_data1 = 275:spacing:300;
s_data2 = 350:spacing:375;
s_data3 = 425:spacing:450;
s_data4 = 500:spacing:525;
s_data = [s_data1, s_data2, s_data3, s_data4];
[s, alpha, beta, s_data, alpha_data, beta_data] = dataMaker(Af0_true, As0_true, lambda_f_true, lambda_s_true, angular_vel_f_true, angular_vel_s_true, s_true, s_data, std);
%alpha_data_tab = [s_data(:), alpha_data(:)]; % Concatenate Vectors
%beta_data_tab = [s_data(:), beta_data(:)]; % Concatenate Vectors
xdata = s_data;
ydata = [alpha_data(:),beta_data(:)];
fun1 = @(v,s_data)v(1)*exp(v(3)*s_data).*cos(v(5)*s_data) + v(2)*exp(v(4)*s_data).*cos(v(6)*s_data);
fun2 = @(v,s_data)v(1)*exp(v(3)*s_data).*sin(v(5)*s_data) + v(2)*exp(v(4)*s_data).*sin(v(6)*s_data);
fun = @(v,s_data)[fun1(v,s_data)', fun2(v,s_data)']; % Composite Function
%initial guess
percent = 0.2; %Percent deviation
x0 = [Af0_true-percent*Af0_true, As0_true-percent*As0_true, lambda_f_true-percent*lambda_f_true, lambda_s_true-percent*lambda_s_true,...
angular_vel_f_true-percent*angular_vel_f_true, angular_vel_s_true-percent*angular_vel_s_true];
%X = lsqcurvefit(fun, x0, alpha_data_tab, beta_data_tab);
X = lsqcurvefit(fun, x0, xdata, ydata);
function [s_total, alpha_total, beta_total, s_data, alpha_data, beta_data] = dataMaker(Af0_true, As0_true, lambda_f_true, lambda_s_true, angular_vel_f_true, angular_vel_s_true, s_total, s_data, std_deg)
s_total = s_total;
std_rad = std_deg*(pi/180); %rad
%True parameter vector
v0 = [Af0_true, As0_true, lambda_f_true, lambda_s_true, angular_vel_f_true, angular_vel_s_true];
alpha_total = v0(1)*exp(v0(3)*s_total).*cos(v0(5)*s_total) + v0(2)*exp(v0(4)*s_total).*cos(v0(6)*s_total);
beta_total = v0(1)*exp(v0(3)*s_total).*sin(v0(5)*s_total) + v0(2)*exp(v0(4)*s_total).*sin(v0(6)*s_total);
s_data = s_data;
alpha_data = v0(1)*exp(v0(3)*s_data).*cos(v0(5)*s_data) + v0(2)*exp(v0(4)*s_data).*cos(v0(6)*s_data);
beta_data = v0(1)*exp(v0(3)*s_data).*sin(v0(5)*s_data) + v0(2)*exp(v0(4)*s_data).*sin(v0(6)*s_data);
alpha_noise = std_rad*randn(1,length(s_data));
beta_noise = std_rad*randn(1,length(s_data));
alpha_data = alpha_data + alpha_noise;
beta_data = beta_data + beta_noise;
end
0 Comments
More Answers (0)
See Also
Categories
Find more on Get Started with Curve Fitting Toolbox in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!