How do I Regression Fit a SinWave to a dataset?

조회 수: 278 (최근 30일)
Clifford Shelton
Clifford Shelton 2012년 4월 30일
댓글: Devon Cogan 2016년 7월 18일
I have a dataset and I want to best fit a sinewave to the plotted data set. This process I think is called a regression...but all the info I come across is about linear regressions only.
Any help would be most appreciated!
  댓글 수: 1
Arjun Jaitli
Arjun Jaitli 2014년 11월 20일
My question is for Wayne King - When you finally plot the fitted curve (yhat) and the actual data (y), is yhat the error or in other words the least square difference between the actual data and the sine fit?

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

채택된 답변

Wayne King
Wayne King 2012년 4월 30일
You need to know what periods you want to fit. You had another post where you talked about fitting city population for a period of 50 years. You did not say how often the data are sampled, I'll assume yearly. Just substitute your data for y (as a column vector)
t = (1:50)';
X = ones(50,3);
X(:,2) = cos((2*pi)/50*t);
X(:,3) = sin((2*pi)/50*t);
y = 2*cos((2*pi)/50*t-pi/4)+randn(size(t));
y = y(:);
beta = X\y;
yhat = beta(1)+beta(2)*cos((2*pi)/50*t)+beta(3)*sin((2*pi)/50*t);
plot(t,y,'b');
hold on
plot(t,yhat,'r','linewidth',2);
If you have the Statistics Toolbox, you can do the same thing with regress()
If you don't know the periods, it is best to use Fourier analysis.
  댓글 수: 2
Clifford Shelton
Clifford Shelton 2012년 5월 1일
Your answer to my question yesterday has been EXTREMELY helpful! Thanks so much. Some issues though:
Let's say I had a data set of city population statistics and it was yearly data for 501 years. And I wanted to see if there was a 25.7 year cycle to that data.
How would I go about fitting a sin wave with a 25.7 year period to my 501 year data set?
Furthermore, when I try to substitute y with my data and create a column vector as suggested in your posted code:
>> X = ones(501,3);
>> X(:,2) = cos((2*pi)/501*t);
>> X(:,3) = sin((2*pi)/501*t);
>> y = Score(:);
>> y = 2*cos((2*pi)/50*t-pi/4)+randn(size(t));
>> y = y(:);
>> beta = X/y;
??? Error using ==> mldivide
Matrix dimensions must agree.
I get the above error. What am I doing wrong in constructing my column vector with the data I call "Score"?
Thank you SO much!
Devon Cogan
Devon Cogan 2016년 7월 18일
The error you're getting means that your column vector "Score" does not have 501 points of data. If you are trying to fit a curve to a data set of only 25.7 years, isolate the 26 data points you want to analyze and change all the "501"'s to "26"'s.

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

추가 답변 (1개)

Richard Willey
Richard Willey 2012년 5월 1일
Here's some simple code that illustrates how to perform nonlinear regression using the 12a release of Statistics Toolbox.
Note: NonLinearModel.fit requires that you provide starting conditions for the various parameters. (Providing good starting conditions helps to ensure that the optimization solvers converge on a global solution rather than a local solution)
%%Generate some data
X = 2* pi*rand(100,1);
X = sortrows(X);
Y = 9 + 7*sin(2*X + 4*pi) + randn(100,1);
scatter(X,Y)
Generate a fit
% Note that we need to pass three sets of input arguments to NonLinearModel
% # The X and Y data
% # A string describing our model
% # Starting conditions for the optimization solvers
% Generate some good starting conditions for the solvers
scatter(X, Y)
hold on
B0 = mean(Y); % Vertical shift
B1 = (max(Y) - min(Y))/2; % Amplitude
B2 = 2; % Phase (Number of peaks)
B3 = 0; % Phase shift (eyeball the Curve)
myFit = NonLinearModel.fit(X,Y, 'y ~ b0 + b1*sin(b2*x1 + b3)', [B0, B1, B2, B3])
% Note that all the coefficient estimates are very good except for b3 where
% any even integer is equally valid
%%look at the complete set of methods
methods(myFit)
%%Generate a plot
hold on
plot(X, myFit.Fitted)
hold off
%%Generate a fit using an alternative syntax
myFit2 = NonLinearModel.fit(X,Y, @(b,x)(b(1) + b(2)*sin(b(3)*x + b(4))), [B0, B1, B2, B3])

카테고리

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

Community Treasure Hunt

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

Start Hunting!

Translated by