How to efficiently perform Fminsearch fitting

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

[Anand Ra]

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

Answer 1

[Matt J]

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

[Anand Ra]

@Matt J

Thank you for the response. However, I am getting a matrix output when I update the function per your say. Expecting an array output from model function. ( Itheory=array).

Also attaching the data set for reference.

Below is the updated code and also the error produced.

%% 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));

Itheory = Imeas.*exp(-Ro*-2*bet_a*k*2*Ro*(1- (y/Ro).^2));

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

Error

Arrays have incompatible sizes for this operation.

Error in Model_Fitting>Ierr (line 117)

err_sq = sum( (Itheory(:)-Imeas(:)).^2);

Error in Model_Fitting (line 75)

err_sq = Ierr(Imeas,Itheory)

[Matt J]

Make sure Imeas and y have the same dimensions, e.g., are both column vectors rather than row vectors.

[Anand Ra]

@Matt J. I changed both to column vectors. That did fix one of the problem . But still getting he following error.

Unrecognized function or variable 'Itheory'.

Error in fminsearch (line 201)

fv(:,1) = funfcn(x,varargin{:});

Error in Model_Fitting (line 82)

[best_R,best_Error] = fminsearch(fun,Ro)

%% 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));

Itheory = Imeas.*exp(-Ro*-2*bet_a*k*2*Ro*(1- (y/Ro).^2));

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

[Matt J]

fun = @(R) Ifunc(y,R,Imeas, bet_a, k, I_size) ;
[best_R,best_Error] = fminsearch(fun,Ro)

[Anand Ra]

@Matt J Got it and fixed it. But its still throwing me the following error:

Unable to perform assignment because the size of the

left side is 1-by-1 and the size of the right side is

1452-by-1.

Error in fminsearch (line 201)

fv(:,1) = funfcn(x,varargin{:});

Error in Model_Fitting (line 84)

[best_R,best_Error] = fminsearch(fun,Ro);

Did verify the variable sizes for consistency. Not sure whats going on.

%% 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;

clearvars;

%% 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 = @(R) Ifunc(y,Ro,Imeas, bet_a, k, I_size) ;

[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));

Itheory = Imeas.*exp(-Ro*-2*bet_a*k*2*Ro*(1- (y/Ro).^2));

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

[Matt J]

Whoops. That should have been,

fun = @(R) Ierr( Imeas,  Ifunc(y,R,Imeas, bet_a, k, I_size)) ;

[Anand Ra]

@Matt. You've been incredible. Thanks a ton!

[Anand Ra]

@Matt J

I am trying to determine two varibales through data fitting by modifying the above code that you helped me fix. And I am unable to fix the error. Below is the code and attached is the data set.

%% Author: Anand Rathnam

% Model Fitting using fminsearch

%%

%{

Objective:

To determine the thickness of copper wire and its coating(PMMA) by curve fitting.

The CT scan decaying intensity data is fitted to the attenuation model to determine the thickness.

%}

%{

Variable/Constants definitions

1. Imeas - Measured CT Intensity (Normalized)

2. y - Number of data points

3. k - wave number

4. beta - refractive index

5. R - Radius of the wire

6. Itheory - Predicted intensity

7. z_Cu = thickness of copper wire

8. z_PMMA = thickness of coating

%}

clc;

close all;

clearvars;

%% Loading and plotting the Measured data

load('Imeas.txt', 'Imeas') % Loading the measured 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,'r-o');

set(gca,'YLim',[0,1.2]);

% set(gca,'XLim',[0,3.5]);

% set(gca,'XTick',x)

grid on

grid minor

xlabel('y');

ylabel('Imeas(Measured Intensity)');

%Initial guess

z_Cu = 1800; % microns

z_PMMA =40; % microns

Io = 1;

%%Absorption coefficient, mu

lamda = 7.1255e-5; %micron

k = (2*pi)/lamda; % (micron)^-1

beta_Cu = 2.590970004E-07; % Refractive index From https://henke.lbl.gov/optical_constants/

beta_PMMA = 5.03136866E-10; % Refractive index From https://henke.lbl.gov/optical_constants/

% Absorption coefficient, mu = 2*beta*k

mu_Cu = beta_Cu*2*k; %(micron)^-1

mu_PMMA = beta_PMMA*2*k; %(micron)^-1

%% Calling the functions

% 1. Calling the "Model prediction" function and plotting

Itheory = Ifunc(z_Cu, z_PMMA,Io,mu_PMMA,mu_Cu);

hold on

title('Model Fitting')

grid on

grid minor

h2=plot(y,Itheory);

options = optimset('Display','iter','PlotFcns',@optimplotfval);

legend([h1,h2],{'Measured Data','Predicted Data'},'Location','NorthWest');

% 2. Calling the "Error function" to calculate the SSE

err_sq = Ierr(Imeas,Itheory');

% 3. Finding the best fitting parameter

fun = @(z_Cu, z_PMMA) Ierr( Imeas, Ifunc(z_Cu, z_PMMA,Io,mu_PMMA,mu_Cu));

%[best_R,best_Error] = fminsearch(fun,Ro)

[best_z_Cu,best_z_PMMA,best_Error] = fminsearch(fun,z_Cu, z_PMMA,options)

% hold on

% h3 = plot(Ro,err_sq,'ro','MarkerFaceColor','r');

%% 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(z_Cu, z_PMMA,Io,mu_PMMA,mu_Cu )

start = 0;

End =1451;

points = 1451; % for micron steps

steps = (z_Cu+z_PMMA)/(points-1);

for i=1:points

Itheory(i) = Io*exp(-1*(mu_PMMA+mu_Cu)*(steps*(i-1)));

end

end

%% 2. Error function -Writing it as a function to allow

%updating the code if we want to monitor the progress of minimization

%{

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

The error is:

Error using fminsearch (line 105)

Argument 3 must be an options structure.

Error in Model_Fitting (line 107)

[best_z_Cu,best_z_PMMA,best_Error] = fminsearch(fun,z_Cu, z_PMMA,options)

[Torsten]

fun = @(z) Ierr( Imeas, Ifunc(z(1), z(2),Io,mu_PMMA,mu_Cu));
%[best_R,best_Error] = fminsearch(fun,Ro)
[best_z_Cu,best_z_PMMA,best_Error] = fminsearch(fun,[z_Cu,z_PMMA],options)

instead of

fun = @(z_Cu, z_PMMA) Ierr( Imeas, Ifunc(z_Cu, z_PMMA,Io,mu_PMMA,mu_Cu));
%[best_R,best_Error] = fminsearch(fun,Ro)
[best_z_Cu,best_z_PMMA,best_Error] = fminsearch(fun,z_Cu, z_PMMA,options)

[Anand Ra]

@Torsten Thank you for the response. I am still getting an error after updating the code per your suggestoom. Below is the error:

Unable to perform assignment because the size of

the left side is 1-by-1 and the size of the right

side is 1-by-1451.

Error in fminsearch (line 201)

fv(:,1) = funfcn(x,varargin{:});

Error in Fitting_two_thickness (line 77)

[best_z_Cu,best_z_PMMA,best_Error] =

fminsearch(fun,[z_Cu,z_PMMA],options)

>>

[Torsten]

Check the sizes of Imeas and Itheory in Ierr:

function err_sq = Ierr(Imeas,Itheory)
% Itheory = Ifunc(y,Ro,Imeas);
size(Imeas)
size(Itheory)
err_sq = sum( (Itheory(:)'-Imeas(:)).^2);
end

My guess is that

function err_sq = Ierr(Imeas,Itheory)
  err_sq = sum( (Itheory(:)-Imeas(:)).^2);
end

is correct.

[Anand Ra]

@Torsten Great, thank you so much for thr fix!

But trying to make sense of the output.

Getting a couple of values for one variable. Is it presenting a range?

Output:

best_z_Cu =

-15.2008 47.4062

best_z_PMMA =

262.5226

best_Error =

1

>> best_z_Cu

[Torsten]

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.

[Anand Ra]

@TorstenThanks a lot!!!