fminsearch for parameter estimation
Show older comments
Hi, folks!
I have been problem in using fminsearch for parameter estimation. The problem is to retrieve parameter (such as, oxygen saturation and hemoglobin concentration) from reflectance data.
My data consists of wavelength (x-axis) and reflectance intensity (y-axis). So far, I am not convinced by result obtained because it is similar with the initial guess. Here is the code, please help me!
data = importdata('R-data.txt');
xdata = data(:,1);
ydata = data(:,2);
data_mua = importdata('mua_all.txt');
mua_hbo = data_mua(:,2);
mua_hb = data_mua(:,3);
mua_water = data_mua(:,4);
mua_lipid = data_mua(:,5);
T = [0.001 0.8];
a = 2.8;
c_water = 0.55;
c_lipid = 0.75;
g = 0.8;
k = 1;
rc = 0.04;
rmus = 1.5;
mus = rmus./(1-g);
zo = 1./rmus;
r1 = sqrt(rc.^2+zo.^2);
mua_blood = T(2).*mua_hbo + (1-T(2)).*mua_hb;
mua = T(1).*mua_blood + c_water.*mua_water + c_lipid.*mua_lipid;
mueff = sqrt(3.*mua.*(mua+rmus));
D = 1./3.*(mua+rmus);
r2 = sqrt(zo.^2.*(1+(4.*a)./3).^2+rc.^2);
zb = 2.*k.*D;
fun = @(T,xdata)(rmus./(4.*pi*(mus+mua))).*(zo.*(mueff+r1.^-1).*(exp(-mueff.*r1)./r1.^2)+...
(zo+2.*zb).*(mueff+r2.^-1).*(exp(-mueff.*r2)./r2.^2));
err = @(T) sum((fun(xdata,T)-ydata).^2);
param = zeros(57,1);
param = fminsearch(err,T);
1 Comment
John D'Errico
on 15 Sep 2020
Edited: John D'Errico
on 15 Sep 2020
param is a vector of length 57. Using fminsearch to search a 57 dimension parameter space is a silly thing to do. The code is not designed to work well above roughly 6-8 dimensions. fminsearch is a reasonable solver for small problems. 57 dimensions is NOT small.
On the other hand, T seems to be a vector of length 2. At least that is how T is initially initialized. So your code is almost impossible to understand because you are doing strange things. What are you trying to solve here?
Other examples of why your code is impossible to run might be in these lines:
fun = @(T,xdata)(rmus./(4.*pi*(mus+mua))).*(zo.*(mueff+r1.^-1).*(exp(-mueff.*r1)./r1.^2)+...
(zo+2.*zb).*(mueff+r2.^-1).*(exp(-mueff.*r2)./r2.^2));
err = @(T) sum((fun(xdata,T)-ydata).^2);
so fun is written as a function of arguments T and xdata, in that order.
But how do you use fun? You pass in xdata and T, in the opposite order. MATLAB does not understand the arguments were reversed.
Honestly, I think you just need to learn to code more carefully. read your own code. Look at what it does. Verify every line, to understand what it does and make sure it is doing what you think it is doing. But who am I to know?
Answers (0)
Categories
Find more on MATLAB 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!