필터 지우기
필터 지우기

How to use fitnlm

조회 수: 64 (최근 30일)
Peter Derivolkov
Peter Derivolkov 2020년 4월 26일
댓글: Image Analyst 2020년 5월 5일
Hello I am trying to fit acurve to my data but I don't know where to begin or what to do? Any help will be appretiated!
This is the graph and code of my data. I have attached the excel file of my data.
A=xlsread('Plague_data.xlsx');
t=A(:,12);
c=A(:,13);
figure (1)
Number.of.cases=plot(t,c,'.','color','b')

채택된 답변

Image Analyst
Image Analyst 2020년 4월 27일
편집: Image Analyst 2020년 4월 27일
Try this:
% Fits a Gaussian to the USA daily deaths for CoVid-19, which includes data for only the left portion of the Gaussian, not the full Gaussian.
%--------------------------------------------------------------------------------------------------------------------------------------------------------
% CLEAN UP - INITIALIZATION STEPS
clc; % Clear the command window.
fprintf('Beginning to run %s.m.\n', mfilename);
close all; % Close all figures (except those of imtool.)
clearvars; % Erase all existing variables.
workspace; % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 20;
%--------------------------------------------------------------------------------------------------------------------------------------------------------
% TRAINING DATA PREPARATION
% Read in data from workbook.
data = readmatrix('Fitnlm_Gaussian_CoVid19.xlsx');
% Get rid of nans
validRows = ~isnan(data(:, 2));
% Get rid of zero counts
validRows = validRows & (data(:, 2) > 0);
data = data(validRows, :)
X = data(:, 1);
Y = data(:, 2);
hFig = figure;
plot(X, Y, 'b.-', 'LineWidth', 2, 'MarkerSize', 30);
grid on;
xlabel('Date', 'FontSize', fontSize);
ylabel('Daily Deaths', 'FontSize', fontSize);
title('CoVid-19 Daily Deaths', 'FontSize', fontSize);
% Make x axis be dates.
datetick('x', 'mmm dd, yyyy', 'keepticks');
ax = gca;
ax.XTickLabelRotation = -45;
% Convert X and Y into a table, which is the form fitnlm() likes the input data to be in.
tbl = table(X(:), Y(:));
%--------------------------------------------------------------------------------------------------------------------------------------------------------
% MODEL DEFINITION
% Define the model as Y = a + b * exp(-(x - c)^2 / d)
% Create an anonymous function for it.
% Note how this "x" of ModelFunction 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.
ModelFunction = @(b, x) b(1) + b(2) * exp(-(x(:, 1) - b(3)).^2/b(4));
%--------------------------------------------------------------------------------------------------------------------------------------------------------
% MODEL CREATION : ESTIMATION OF PARAMETERS
% Guess starting model values to start with. Just make your best guess.
% These are just starting points for the coefficients and will be adjusted during the fit to produce the real coefficients.
beta0 = [0, max(X), mean(X), var(X)]; % A guess at what the coefficients will be.
% Now the next line is where the actual model computation is done.
mdl = fitnlm(tbl, ModelFunction, 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'}
%--------------------------------------------------------------------------------------------------------------------------------------------------------
% MODEL PLOTTING : PLOTTING FITTED/ESTIMATED VALUES OVER THE WHOLE RANGE OF X.
% Let's do a fit, but let's get more points on the fit, beyond just the widely spaced training points,
% so that we'll get a much smoother curve.
xFitted = linspace(min(X), max(X), 1920); % Let's use 1920 points, which will fit across an HDTV screen about one sample per pixel.
% Create smoothed/regressed data using the model:
yFitted = ModelFunction(coefficients, xFitted(:));
% yFitted = coefficients(1) + coefficients(2) * exp(-(xFitted - coefficients(3)).^2 / coefficients(4));
% Now we know that since it's a Gaussian there should not be any negative values. So clip to 0
yFitted = max(0, yFitted);
% Now we're done and we can plot the smooth model as a red line going through the noisy blue markers.
hold on;
plot(xFitted, yFitted, 'r-', 'LineWidth', 2);
grid on;
title('CoVid-19 Daily Deaths in the USA -- Exponential Regression with fitnlm()', 'FontSize', fontSize);
% Put a line at the peak.
darkGreen = [0, 0.5, 0];
xline(coefficients(3), 'Color', darkGreen, 'LineWidth', 3);
% Put up text for the green line.
str = sprintf(' Peak at %s', datestr(coefficients(3)));
text(coefficients(3), 100, str, 'Color', darkGreen, 'FontSize', 14, 'FontWeight', 'bold');
legendHandle = legend('Actual Y', 'Fitted Y', 'Location', 'northwest');
legendHandle.FontSize = 25;
hFig.WindowState = 'maximized';
fprintf('Done running %s.m.\n', mfilename);
  댓글 수: 2
Peter Derivolkov
Peter Derivolkov 2020년 4월 27일
Thank you so much this is perfect!
Image Analyst
Image Analyst 2020년 5월 5일
Attached is a new version where it asks you if you want to allow an offset on the Gaussian or not. Plus it uses a bar chart and has more recent data.

댓글을 달려면 로그인하십시오.

추가 답변 (1개)

Keegan Carvalho
Keegan Carvalho 2020년 4월 27일
Firstly, you need to know what fitting curve you are looking for; linear, quadratic, etc. I've gone with polyfit which you can later use to evaluate on the variables T and C.
A=xlsread('Plague_data.xlsx');
t=A(:,12);
c=A(:,13);
C = isnan(c);
T = isnan(t); % polyfit cannot operate if variables have NaN values
fit1 = polyfit(t(~T),c(~C),1); % polyfit can be used and NaN values are omitted
You can start of with the above and later evaluate using polyval. Read about polyfit, fit.
Hope this helps!
  댓글 수: 2
Image Analyst
Image Analyst 2020년 4월 27일
If it's a typical pandemic curve, it's probably like the left portion of a Gaussian.
Peter Derivolkov
Peter Derivolkov 2020년 4월 27일
Thank you so much I appreciate it!

댓글을 달려면 로그인하십시오.

카테고리

Help CenterFile Exchange에서 Biological and Health Sciences에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by