필터 지우기
필터 지우기

How to fit a curve using "power" fitting or "custom fitting"?

조회 수: 8 (최근 30일)
ishita agrawal
ishita agrawal 2017년 9월 16일
댓글: ishita agrawal 2017년 9월 22일
I have data which I need to fit using following equation:
y= f(x)= a*x^b+c.
Code:
r=scatter(npp7,lk_2k1);
r.MarkerEdgeColor = 'r';
r.MarkerFaceColor = [0.9 0.9 0.9];
hold on
% Power fit - %y=f(x)=a*x^b+c
a=28.54;
b=0.4634
c=-3.289;
x = data(:);
y = a*x^b+c;
%f = fit(x,y);
p = plot(x,y)
p(1).LineWidth = 2;
c = p.Color;
p.Color = 'r';
It shows: Error using ^ One argument must be a square matrix and the other must be a scalar. Use POWER (.^) for elementwise power.
But if I use(.^), it shows multiple fit lines as shown in attached figure. I want just one fit line for same equation.

채택된 답변

Image Analyst
Image Analyst 2017년 9월 16일
I attach a demo to do an exponential fit. Adapt as needed.
  댓글 수: 3
Image Analyst
Image Analyst 2017년 9월 16일
편집: Image Analyst 2017년 9월 16일
Below is the full demo:
% Uses fitnlm() to fit a non-linear model (a power law 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;
% Create the X coordinates: 30 points from 0.01 to 20, inclusive.
X = linspace(0.01, 20, 30);
% Define function that the X values obey.
a = 10 % Arbitrary sample values I picked.
b = 0.4
c = 2
Y = a * X .^ b + c; % Get a vector. No noise in this Y yet.
% Add noise to Y.
Y = Y + 0.8 * randn(1, length(Y));
% 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 * (x .^ b) + c
% 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) * x(:, 1) .^ + b(2) + b(3);
beta0 = [10, .4, 2]; % 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) * X .^ coefficients(2) + 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('Power Law Regression with fitnlm()', 'FontSize', fontSize);
xlabel('X', 'FontSize', fontSize);
ylabel('Y', 'FontSize', fontSize);
legendHandle = legend('Noisy Y', 'Fitted Y', 'Location', 'north');
legendHandle.FontSize = 25;
message = sprintf('Coefficients for Y = a * X ^ b + c:\n a = %8.5f\n b = %8.5f\n c = %8.5f',...
coefficients(1), coefficients(2), coefficients(3));
text(8, 15, message, 'FontSize', 23, 'Color', 'r', 'FontWeight', 'bold', 'Interpreter', 'none');
% Set up figure properties:
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0, 0.04, 1, 0.96]);
% 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')
If you need more help, attach npp7 and lk_2k1.
ishita agrawal
ishita agrawal 2017년 9월 22일
This is good enough. Thank you so much.

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

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Linear and Nonlinear Regression에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by