Steps for estimating alpha in 1/f^alpha noise using OLS

2 次查看（过去 30 天）
Ulrik William Nash
Ulrik William Nash 2017-5-28
回答： Nachiket Katakkar 2017-6-1
I have a signal (the vector "pink"), which I know is pink noise, 1/f^alfa, where alpha = 1. But let's assume alpha is unknown, and I need to estimate alpha. In other words, Let's assume I just have the signal.
I use the following procedure to get into a position to estimate alpha using OLS:
F = abs(fft(pink)); x = log(1:(length(pink)/2)); y = log(F(1:(length(pink)/2)));
Now I run OLS regression:
[r,m,b] = regression(x,y)
where I get
r = -0.5948, m = -0.4915, b = 6.8223.
I would have thought my estimate of alpha was m, but m is not equal to 1. So, what am I doing wrong? Is there an additional step I am forgetting, or is my procedure plain wrong?
  3 个评论
显示 1更早的评论
Ulrik William Nash
Ulrik William Nash 2017-5-28
编辑：John D'Errico 2017-5-28
Hi John,
I didn't expect that twist.
Here is the code for the particular function:
function [r,m,b] = regression(targets,outputs,flag)
%REGRESSION Linear regression.
%
% <a href="matlab:doc regression">regression</a> calculates the linear regression between each element
% of the network response and the corresponding target.
%
% [R,M,B] = <a href="matlab:doc regression">regression</a>(T,Y) takes cell array or matrix targets T and
% output Y, each with total matrix rows of N, and returns the linear
% regression for each of the N rows: the regression values R, slopes M,
% and y-intercepts B.
%
% <a href="matlab:doc regression">regression</a>(T,Y,'one') returns scalar R, M and B values across all
% rows of targets and outputs.
%
% Here a feedforward network is trained and regression performed on its
% targets and outputs.
%
% [x,t] = <a href="matlab:doc simplefit_dataset">simplefit_dataset</a>;
% net = <a href="matlab:doc feedforwardnet">feedforwardnet</a>(10);
% net = <a href="matlab:doc train">train</a>(net,x,t);
% y = net(x);
% [r,m,b] = <a href="matlab:doc regression">regression</a>(t,y)
%
% See also PLOTREGRESSION
% Copyright 2010 The MathWorks, Inc.
if nargin < 2, error(message('nnet:Args:NotEnough')); end
if iscell(targets), targets = cell2mat(targets); end
if iscell(outputs), outputs = cell2mat(outputs); end
if all(size(targets) ~= size(outputs))
error(message('nnet:NNData:TYMismatch'))
end
if (nargin >= 3) && ischar(flag) && strcmp(flag,'one')
targets = targets(:)';
outputs = outputs(:)';
end
[N,Q] = size(targets);
m = zeros(N,1);
b = zeros(N,1);
r = zeros(N,1);
for i=1:N
t = targets(i,:);
y = outputs(i,:);
ignore = find(isnan(t) | isnan(y));
t(ignore) = [];
y(ignore) = [];
Quse = Q - length(ignore);
h = [t' ones(size(t'))];
yt = y';
rankStatus = warning('off','MATLAB:rankDeficientMatrix');
rankRestore = onCleanup(@() warning(rankStatus));
theta = h\yt;
m(i) = theta(1);
b(i) = theta(2);
yn = y - mean(y);
tn = t - mean(t);
sty = std(yn);
stt = std(tn);
r(i) = yn*tn'/(Quse - 1);
if (sty~=0)&&(stt~=0)
r(i) = r(i)/(sty*stt);
end
end
Ulrik William Nash
Ulrik William Nash 2017-5-28
I just double-checked that the "regression" and "regress" functions compute the same output, which they do, except the first also calculates r.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Nachiket Katakkar
The regression function is part of the Neural Network Toolbox:
http://www.mathworks.com/help/nnet/ref/regression.html
The coefficient "m" is the slope of the curve between the inputs "x" and "y". The example on the documentation page shows a value of m equal to 1, however, this would not necessarily be the case for any regression problem.
Consider this example with random data:
>> pink = rand(1,100);
F = abs(fft(pink));
x = log(1:(length(pink)/2));
y = log(F(1:(length(pink)/2)));
[r,m,b] = regression(x,y)
r =
-0.3448
m =
-0.3469
b =
1.7326
You can visualize the relationship between "x" and 'y" using:
>> plotregression(x,y)
Also notice, that in your case the value of "b", the offset of regression is also quite high. The regression function seems to be working as expected, so I would recommend confirming that the algorithm you are using is correct.

类别

在 Help Center 和 File Exchange 中查找有关 Linear Regression 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by