How to efficiently perform Fminsearch fitting

6 views (last 30 days)
Hello
I am looking to make the below code efficient. So that I can update when required. Right now its clunky and I findi it difficult to make sense of the fitting process.
The goal is to fit the data ( Itheory) inot the model equation to display:
  1. The best Ro
  2. plot Imeas vs I theory ( to compare the fit)
  3. SSE
  4. montior the progress of optimization using the syntax
[x,fval,exitflag,output] = fminsearch(fun,x0,options)
Below is the code and also attaching the data set
%% Author: Anand Rathnam
%%
%{
Objective:
To determine the radius of the copper wire by fitting the CT
scan intensity data to the absorption intenisty decay model
%}
%{
Variable/Constants definitions
1. Imeas - Measured Intensity
2. y - size of the measured data. Imeas is a function y
3. k - wave number
4. bet_a - scattering coeff
5. R - Radius of the wire
6. Itheory - Predicted intensity
%}
%% Loading and plotting the Measured data
load('Imeas.txt', 'Imeas') % Loading the intensity data
Imeas=double(Imeas); % updating the data type to double to faciliatte computation
I_size = numel(Imeas); % Finding the number of Imeas array elements
y = 0:I_size-1;
% Plotting the measured data
figure(1)
clf
h1=plot(y,Imeas);
% set(gca,'XTick',x)
xlabel(' y');
ylabel('Imeas (Measured Intensity)');
%% Known values
bet_a = 7.30026262E-10;
k =( 2*pi)/0.00000008;
%Initiall guess
Ro = 0.001; % 1 micrometer as the initial radius guess
%% Calling the functions
% 1. Calling the "Model prediction" function and plotting
for i = 1:numel(Imeas)
Itheory = @(Ro)(Imeas(i)*exp(-Ro*-2*bet_a*k*2*Ro*(1-(y(i)/Ro).^2)));
Ro = 0.001; % Initial guess value
[best_R(i),best_Error(i)] = fminsearch(Itheory,Ro);
I_th(i) = Imeas(i).*exp(-best_R(i)*-2*bet_a*k*2*best_R(i)*(1- (y(i)/best_R(i)).^2));
end
% Finally, Error is:
hold on
h2=plot(y,I_th);
best_R(I_size)
best_Error(I_size)
err_sq = sum((I_th(:)-Imeas(:)).^2) % No need to have another sub-function here
  1 Comment
Anand Ra
Anand Ra on 1 May 2022
Above code was a modification based on some help in this forum. Below is my initial code. I would greatly appreciate if this can be fixed.
%% Author: Anand Rathnam
%%
%{
Objective:
To determine the radius of the copper wire by fitting the CT
scan intensity data to the absorption intenisty decay model
%}
%{
Variable/Constants definitions
1. Imeas - Measured Intensity
2. y - size of the measured data. Imeas is a function y
3. k - wave number
4. bet_a - scattering coeff
5. R - Radius of the wire
6. Itheory - Predicted intensity
%}
clc;
close all;
%% Loading and plotting the Measured data
load('Imeas.txt', 'Imeas') % Loading the intensity data
Imeas=double(Imeas); % updating the data type to double to faciliatte computation
I_size = numel(Imeas); % Finding the number of Imeas array elements
y = 0:I_size-1;
% Plotting the measured data
figure(1)
clf
h1=plot(y,Imeas);
set(gca,'YLim',[0,1.2]);
% set(gca,'XLim',[0,3.5]);
% set(gca,'XTick',x)
xlabel(' y');
ylabel('Imeas ( Measured Intensity');
%% Known values
bet_a = 7.30026262E-10;
k =( 2*pi)/0.00000008;
%Initiall guess
Ro = 0.000001; % 1 micrometer as the initial radius guess
%% Calling the functions
% 1. Calling the "Model prediction" function and plotting
Itheory = Ifunc(y,Ro,Imeas, bet_a, k, I_size);
hold on
h2=plot(y,Itheory);
% legend([h1,h2],{'Data','First Guess'},'Location','NorthWest');
% 2. Calling the "Error function" to calculate the SSE
err_sq = Ierr(Imeas,Itheory)
% 3. Finding the best fitting parameter
fun = @Itheory;
[best_R,best_Error] = fminsearch(fun,Ro)
%% 1: Function for the "Model prediction"
%{
This function:
1. Contains the model that predicts the data
2. Takes the input arguments: initial guess, y, Imeas
3. returns an array of predicted data: Itheory
%}
function Itheory = Ifunc(y,Ro,Imeas, bet_a, k, I_size)
% Itheory = Imeas.*exp(-2*bet_a*k*2*R.*(1- (y/R).^2));
for i=1:I_size
Itheory(i) = Imeas(i).*exp(-Ro*-2*bet_a*k*2*Ro*(1- (y(i)/Ro).^2));
end
return
end
%% 2. Error function
%{
This function:
1. Calculates the difference between the data and the model's prediction
2. Calulates and returns sums of squared error (SSE)
3. The SSE calculated is between our model and the data for our first guess
%}
function err_sq = Ierr(Imeas,Itheory)
% Itheory = Ifunc(y,Ro,Imeas);
err_sq = sum( (Itheory(:)-Imeas(:)).^2);
end

Sign in to comment.

Accepted Answer

Matt J
Matt J on 1 May 2022
There's no need for a for loop in Ifunc
function Itheory = Ifunc(y,Ro,Imeas, bet_a, k, I_size)
% Itheory = Imeas.*exp(-2*bet_a*k*2*R.*(1- (y/R).^2));
Itheory = Imeas.*exp(-Ro*-2*bet_a*k*2*Ro*(1- (y/Ro).^2));
end
  14 Comments
Torsten
Torsten on 8 May 2022
fun = @(z) Ierr( Imeas, Ifunc(z(1), z(2),Io,mu_PMMA,mu_Cu));
[best_z,best_Fval,best_ExitFlag] = fminsearch(fun,[z_Cu,z_PMMA],options)
best_z_Cu = best_z(1);
best_z_PMMA = best_z(2);
instead of
fun = @(z) Ierr( Imeas, Ifunc(z(1), z(2),Io,mu_PMMA,mu_Cu));
[best_z_Cu,best_z_PMMA,best_Error] = fminsearch(fun,[z_Cu,z_PMMA],options)
fminsearch and all other solvers work with vectors for the solution variables (in your case [z_Cu,z_PMMA]) instead of handling the solution variables separately. Thus you always have to extract your scalar variables from a (in your case) 2-dimensional vector.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!