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

7 views (last 30 days)
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);

Accepted Answer

Sulaymon Eshkabilov
Sulaymon Eshkabilov on 1 May 2022
Edited: Sulaymon Eshkabilov on 1 May 2022
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
  2 Comments
Anand Ra
Anand Ra on 1 May 2022
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)
>>

Sign in to comment.

More Answers (2)

Sulaymon Eshkabilov
Sulaymon Eshkabilov on 30 Apr 2022
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
...
  1 Comment
Anand Ra
Anand Ra on 1 May 2022
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)

Sign in to comment.


Sulaymon Eshkabilov
Sulaymon Eshkabilov on 1 May 2022
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.
  3 Comments
Anand Ra
Anand Ra on 1 May 2022
@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.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!