Want to use 'fitnlm' to fit the surviving function of a range of distributions to my data

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.

[Image Analyst]

Make it easy for us to help you by attaching a .mat file with your data in it.

[mel1708]

Sure thing, here it is. Sorry I took a while to reply. It was overnight for me here in Australia. This is one of the files of data I have, but the others are the same type of data.

Answer 1

[Image Analyst]

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') 

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?

Answer 2

[Image Analyst]

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') 

[mel1708]

Thank you so much for all of your help! I have incorporated the exponential fit into my code successfully. I have tried doing the same to show the Rayleigh distribution for all of the data sets I have (i.e. performing the fit for each data set imported by looping over k) but I can't seem to get the fit working in the for loop.

This is how I've used your code, but 'mdl' and hence the calculation of sigma doesn't seem to be working in the loop because when it comes to calculating yfitRayl{k} it gives me an error stating incorrect dimensions for matrix multiplication. The mdl and sigma functions are only giving me results for k = 1. Am I using this code correctly inside the for loop?

%==============================
%   RAYLEIGH DISTRIBUTION
%==============================
for k = 1:length(theFiles)
    modelfun = @(sigma2,x) (x/sigma2) .* exp(-x.^2 / (2 * sigma2));  
    beta0 = 1;
    mdl{k} = fitnlm(tbl, modelfun, beta0);

    sigma2Estimate{k} = mdlRayl{k}.Coefficients{:, 'Estimate'}
    sigma{k} = sqrt(sigma2Estimate{k});
    yfitRayl{k} = (dose{k} / sigma{k}^2) * exp(-dose{k} / (2 * sigma{k}^2));
    
    figure(11);
    hold on
    plot(dose{k},sf{k},'o','Color',col(k,:),'Marker','o','LineStyle','none','DisplayName',txt);
    h = errorbar(dose{k},sf{k},errs{k},'Marker','none','LineStyle','none','Color',[col(k,1), col(k,2), col(k,3)],'HandleVisibility','off');
    legend('Location','southoutside');
    xlabel('Dose (Gy)'); ylabel('Surviving Fraction'); grid on; legend show;
    plot(dose{k},yfitRayl{k},'Color',col(k,:),'HandleVisibility','off');
end
    hold off

[Image Analyst]

Attach the one/first file that it fails on.

Also, I don't the code where you're updating tbl according to the new file. If you don't do that then why would mdl{k} be any different on any iteration? It will be the same.

[mel1708]

Here is all of the data I am currently testing the code on. So each 6x1 double in the cell array is a group of data like I attached in the first version of survivaldata.mat.

Oops, here is the code I used to define tbl, which I put at the top of the section so I can just call it for each distribution. This is why it's not in the code block I posted above.

    tbl = table(dose{k},sf{k});
    

So basically I want to perform a whole new Rayleigh fit for each k and plot them all onto one plot. I coded the exponential fit the same way and ended up with quite a nice looking fit so thankyou for suggesting this.

[Image Analyst]

So it's really after the "for" but before the "fitnlm"? Can you just attach the whole m-file?

[mel1708]

The exponential fit you coded above is on lines 191 to 208 and appears to be working fine. The next section from 213 to 231 is the section I can't get to work.