Want to use 'fitnlm' to fit the surviving function of a range of distributions to my data
조회 수: 3 (최근 30일)
이전 댓글 표시
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
Image Analyst
2020년 1월 3일
Make it easy for us to help you by attaching a .mat file with your data in it.
답변 (2개)
Image Analyst
2020년 1월 3일
편집: Image Analyst
2020년 1월 3일
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?
댓글 수: 0
Image Analyst
2020년 1월 3일
편집: Image Analyst
2020년 1월 3일
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')
댓글 수: 5
Image Analyst
2020년 1월 5일
So it's really after the "for" but before the "fitnlm"? Can you just attach the whole m-file?
참고 항목
카테고리
Help Center 및 File Exchange에서 Probability Distributions에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!