Want to use 'fitnlm' to fit the surviving function of a range of distributions to my data
3 views (last 30 days)
Show older comments
I have multiple sets of files of the same type of data (radiation dose and surviving fraction of cells) which I want to fit distributions (e.g. Rayleigh, Exponential, Log logistic, and Pareto). I figured I'd want to use the fitnlm function like from the documentation..
modelfun = @(b,x)b(1) + b(2)*x(:,1).^b(3) + ...
b(4)*x(:,2).^b(5);
beta0 = [-50 500 -1 500 -1];
mdl = fitnlm(X,y,modelfun,beta0)
If I were to try and put my own variables in and involve a for loop to include all the data sets I have
for k = 1:N
cdfRayl{k} = @(sigma(k), dose(k)) 1 - exp(-(dose(k).^2)/(2*sigma(k)^2));
beta0 = [0];
yfitRayl{k} = fitnlm(dose{k},sf{k},cdfRayl{k},beta0);
end
However, when I try and run this it returns the error
"Invalid expression. When calling a function or indexing a variable, use parentheses. Otherwise, check for mismatched delimiters."
This error is referring to the cdfRayl line. I think it is having trouble with the @(sigma(k), dose(k)) part.
Is anyone able to suggest another way to code this? I used fitlm to do the same thing for gamma and poisson distributions for which the code looked like this
for k = 1:length(theFiles)
[lmgammaFit,devGamma,statsGamma] = glmfit(dose{k},sf{k},'Gamma');
yfitGamma{k} = glmval(lmgammaFit,dose{k},'reciprocal');
end
I'm basically trying to do the same thing but with the distributions listed above.
2 Comments
Image Analyst
on 3 Jan 2020
Make it easy for us to help you by attaching a .mat file with your data in it.
Answers (2)
Image Analyst
on 3 Jan 2020
Edited: Image Analyst
on 3 Jan 2020
Here's the code:
% Uses fitnlm() to fit a non-linear model (a Rayleigh curve) through noisy data.
% Requires the Statistics and Machine Learning Toolbox, which is where fitnlm() is contained.
% Initialization steps.
clc; % Clear the command window.
close all; % Close all figures (except those of imtool.)
clear; % Erase all existing variables. Or clearvars if you want.
workspace; % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 20;
s = load('survivaldata.mat')
dose = s.dose;
sf = s.sf;
% Create the X coordinates.
X = dose;
% Y is the survival factor.
Y = sf; % Get a vector.
% Now we have noisy training data that we can send to fitnlm().
% Plot the noisy initial data.
plot(X, Y, 'b*', 'LineWidth', 2, 'MarkerSize', 15);
grid on;
% Convert X and Y into a table, which is the form fitnlm() likes the input data to be in.
tbl = table(X(:), Y(:));
% Define the model as Y = a + exp(-b*x)
% Note how this "x" of modelfun is related to big X and big Y.
% x((:, 1) is actually X and x(:, 2) is actually Y - the first and second columns of the table.
modelfun = @(sigma2,x) (x(:, 1)/sigma2) .* exp(-x(:, 1).^2 / (2 * sigma2));
beta0 = 1; % Guess values to start with. Just make your best guess.
% Now the next line is where the actual model computation is done.
mdl = fitnlm(tbl, modelfun, beta0);
% Now the model creation is done and the coefficients have been determined.
% YAY!!!!
% Extract the coefficient values from the the model object.
% The actual coefficients are in the "Estimate" column of the "Coefficients" table that's part of the mode.
sigma2Estimate = mdl.Coefficients{:, 'Estimate'}
sigma = sqrt(sigma2Estimate)
% Create smoothed/regressed data using the model:
yFitted = (X / sigma^2) .* exp(-X / (2 * sigma^2));
% Now we're done and we can plot the smooth model as a red line going through the noisy blue markers.
hold on;
plot(X, yFitted, 'r-', 'LineWidth', 2);
grid on;
title('Exponential Regression with fitnlm()', 'FontSize', fontSize);
xlabel('X', 'FontSize', fontSize);
ylabel('Y', 'FontSize', fontSize);
legendHandle = legend('Noisy Y', 'Fitted Y', 'Location', 'north');
legendHandle.FontSize = 25;
% Set up figure properties:
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);
% Get rid of tool bar and pulldown menus that are along top of figure.
% set(gcf, 'Toolbar', 'none', 'Menu', 'none');
% Give a name to the title bar.
set(gcf, 'Name', 'Demo by ImageAnalyst', 'NumberTitle', 'Off')
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/259679/image.png)
As you can see, because your data does not go up then down like a Rayleigh curve is supposed to do, the fit is not so good. Are you sure you don't want to fit it to an exponential decay curve?
0 Comments
Image Analyst
on 3 Jan 2020
Edited: Image Analyst
on 3 Jan 2020
For what it's worth, here is the fit to an exponential decay:
% Uses fitnlm() to fit a non-linear model (an exponential decay curve) through noisy data.
% Requires the Statistics and Machine Learning Toolbox, which is where fitnlm() is contained.
% Initialization steps.
clc; % Clear the command window.
close all; % Close all figures (except those of imtool.)
clear; % Erase all existing variables. Or clearvars if you want.
workspace; % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 20;
s = load('survivaldata.mat')
dose = s.dose;
sf = s.sf;
% Create the X coordinates
X = dose;
% Get Y
Y = sf; % Get a vector.
% Now we have noisy training data that we can send to fitnlm().
% Plot the noisy initial data.
plot(X, Y, 'b*', 'LineWidth', 2, 'MarkerSize', 15);
grid on;
% Convert X and Y into a table, which is the form fitnlm() likes the input data to be in.
tbl = table(X(:), Y(:));
% Define the model as Y = a + exp(-b*x)
% Note how this "x" of modelfun is related to big X and big Y.
% x((:, 1) is actually X and x(:, 2) is actually Y - the first and second columns of the table.
modelfun = @(b,x) b(1) * exp(-b(2)*x(:, 1)) + b(3);
beta0 = [11, .5, 6]; % Guess values to start with. Just make your best guess.
% Now the next line is where the actual model computation is done.
mdl = fitnlm(tbl, modelfun, beta0);
% Now the model creation is done and the coefficients have been determined.
% YAY!!!!
% Extract the coefficient values from the the model object.
% The actual coefficients are in the "Estimate" column of the "Coefficients" table that's part of the mode.
coefficients = mdl.Coefficients{:, 'Estimate'}
% Create smoothed/regressed data using the model:
yFitted = coefficients(1) * exp(-coefficients(2)*X) + coefficients(3);
% Now we're done and we can plot the smooth model as a red line going through the noisy blue markers.
hold on;
plot(X, yFitted, 'r-', 'LineWidth', 2);
grid on;
title('Exponential Regression with fitnlm()', 'FontSize', fontSize);
xlabel('X', 'FontSize', fontSize);
ylabel('Y', 'FontSize', fontSize);
legendHandle = legend('Noisy Y', 'Fitted Y', 'Location', 'north');
legendHandle.FontSize = 25;
% Set up figure properties:
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);
% Get rid of tool bar and pulldown menus that are along top of figure.
% set(gcf, 'Toolbar', 'none', 'Menu', 'none');
% Give a name to the title bar.
set(gcf, 'Name', 'Demo by ImageAnalyst', 'NumberTitle', 'Off')
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/259682/image.png)
5 Comments
Image Analyst
on 5 Jan 2020
So it's really after the "for" but before the "fitnlm"? Can you just attach the whole m-file?
See Also
Categories
Find more on Probability Distributions 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!