Perform fminsearch to find the best value for the variable in a model

Hello

I trying to fit data set into a model to predict the best value of a variable using fminsearch.

I am getting a couple of errors and I am unable to trouble. Kindly requesting for assistance.

Below is my code, followed by errors generated when I ran the code.

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

hold on

h2=plot(y,Itheory,'r-');

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,options)

%% 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,R,Imeas)

Itheory = Imeas.*exp(-2*bet_a*k*2*R*(1- (y/R)^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

Errors generated:

Unrecognized function or variable 'bet_a'.

Error in Model_Fitting>Ifunc (line 104)

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

Error in Model_Fitting (line 73)

Itheory = Ifunc(y,Ro,Imeas);

Answer 1

[Sulaymon Eshkabilov]

Note the following two points in your code that must be identical and input variable names must match as defined above:

Itheory = Ifunc(y,Ro,Imeas, bet_a, k, I_size);

function Itheory = Ifunc(y,R,Imeas, bet_a, k, I_size)  
...
    

You have initiated: Ro = 0.000001;

But calling in your fcn by: R that is not defined.

Moreover, you have obtained numerical values of Itheory and later calling it again as a function. You'd need to work here a bit differently. In this case, it is better to employ the function handle approach.

This small code takes over the two sub-functions.

...
for ii = 1:numel(Imeas)
Itheory = @(Ro)(Imeas(ii)*exp(-Ro*-2*bet_a*k*2*Ro*(1-(y(ii)/Ro)^2));
Ro = 0.000001;       % Initial guess value
[best_R(ii),best_Error(ii)] = fminsearch(Itheory,Ro);
I_th(ii) = (Imeas(ii)*exp(-best_R(ii)*-2*bet_a*k*2*best_R(ii)*(1-(y(ii)/best_R(ii)))^2));
end
% Finally, Error is:
err_sq = sum((I_th(:)-Imeas(:)).^2);  % No need to have another sub-function here 

[Anand Ra]

@Sulaymon EshkabilovThanks, I fixed it. I am trying to find the best value for Ro.

But still getting errors.

I have also attached the data set for your reference if you want to run the code and verify.

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

%}

%% Function definitions

%{

1.

%}

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

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

[Anand Ra]

Error produced when running the above code:

Unrecognized function or variable 'Itheory'.

Error in fminsearch (line 201)

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

Error in Model_Fitting (line 89)

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

>>

Answer 2

[Sulaymon Eshkabilov]

There are a couple potential of errs in your code. (1) Ifunc fcn has to have two more input variables, viz. bet_a, k and thus, while calling it these variables to be included, e.g.:

...
% 1. Calling the "Model prediction" function and plotting
Itheory = Ifunc(y,R,Imeas, bet_a, k);  % Note: the input variables
hold on
h2=plot(y,Itheory,'r-');
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,options)
 
function Itheory = Ifunc(y,R,Imeas, bet_a, k)   % Note: the input variables
Itheory = Imeas.*exp(-2*bet_a*k*2*R*(1- (y/R)^2));
end
...

[Anand Ra]

@Sulaymon Eshkabilov

Thanks. I am still facing errors. Also I made some changes to the code.

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

%}

%% Function definitions

%{

1.

%}

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,R,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(-R*-2*bet_a*k*2*R*(1- (y(i)/R).^2));

end

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 produced:

Unrecognized function or variable 'Itheory'.

Error in fminsearch (line 201)

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

Error in Model_Fitting (line 89)

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

Answer 3

[Sulaymon Eshkabilov]

One pranthesis was missing:

...
for ii = 1:numel(Imeas)
Itheory = @(Ro)(Imeas(ii)*exp(-Ro*-2*bet_a*k*2*Ro*(1-(y(ii)/Ro)^2)));  %  ) missing
Ro = 0.000001;       % Initial guess value
[best_R(ii),best_Error(ii)] = fminsearch(Itheory,Ro);
I_th(ii) = (Imeas(ii)*exp(-best_R(ii)*-2*bet_a*k*2*best_R(ii)*(1-(y(ii)/best_R(ii)))^2));
end

It is working correctly and without any err. Verify your formulations for Itheory.

[Anand Ra]

@Sulaymon EshkabilovThank you. Function handle approach is always confusing to me.

But yes, it works. Thanks

Btw what is the best Ro? and what is the error? the last values?

Also, I am looking to monitor the iteration

1. Using the syntax

[x,fval,exitflag,output] = fminsearch(fun,x0,options)

2. And also plot the measured vs theoretical to find the fit.

3. find the best Ro and error

Yes, I theory model- function is weird, might need retrospection. ( Will look into it, thanks for letting me know)

[Sulaymon Eshkabilov]

Yes, the last value is the found best Ro value and min err.

So, you'd accept the ans., right :)

All the best.

[Anand Ra]

@Sulaymon Eshkabilov Thanks will do. Appreciate the service.

As I mentioned, I am looking to display

1. Monitor optimization using the syntax

[x,fval,exitflag,output] = fminsearch(fun,x0,options)

2. And also plot the measured vs theoretical to find the fit.

3. find the best Ro and error

The fminsearch inside the forloop makes it look clunky for me to modify anything.