Matlab:尝试通过直方图盒计数从时间序列估计多重分形谱
Matlab: trying to estimate multifractal spectrum from time series by histogram box-counting
我正在使用耶鲁分形页面中的方法:
http://classes.yale.edu/fractals/MultiFractals/Moments/TSMoments/TSMoments.html
这组讲座幻灯片(幻灯片 32)也对此进行了阐述:
http://multiscale.emsl.pnl.gov/docs/multifractal.pdf
这个想法是你得到一个数据集,并通过许多直方图检查它,增加条的数量,即分辨率。一旦分辨率足够高,一些条形就会取零值。在每个阶段,我们取落入每个直方图箱(忽略任何零值箱)的结果数量,将其除以数据集的总大小,并将其提高到幂 q。此操作给出了给定时刻和 bin 大小的 'partition function'。引用上面链接的教程:“提供了对非均匀性的选择性表征
测量,正 q 突出最密集的区域,负 q 突出最平滑的区域。
所以我在 Matlab 中使用直方图函数,遍历 bin 大小,对所有非零 bin 内容求和,等等。但是我的分区函数输出数组只是一堆 1。我看不出出了什么问题,其他人可以吗?
英特尔、思科、苹果和其他公司的数据可在同一耶鲁网站上获得:yale.edu/fractals/MultiFractals/Finf(a)/Finf(a).html
N.B。 intel
指的是我原先用作数据集的intel股价
lower = 1; %set lowest level of histogram resolution (bin size)
upper = 300; %set highest level of histogram resolution (bin size)
qlow = -20; %set lowest moment exponent
qhigh = 20; %set highet moment exponent
qstep = 0.25; %set step size between moment exponents
qn= ((qhigh-qlow)/qstep) + 1; %calculates number of steps given qlow, qhigh, qstep
qvalues= linspace(qlow, qhigh, qn); %creates a vector of q values given above parameters
m = min(intel); %find the maximum of the dataset
M = max(intel); %find the minimum of the dataset
for Q = 1:length(qvalues) %loop over moment exponents q
for k = lower:upper %loop over bin sizes
counts = hist(intel, k); %unpack all k histogram height values into 'counts'
counts(counts==0) = []; %delete all zero values in ''counts
Zq = counts ./ length(intel);
Zq = Zq .^ qvalues(Q);
Zq = sum(Zq);
partitions(k-(lower-1), Q) = Zq; %store Zq in the kth row and the Qth column of 'partitions'
end
end
您的代码似乎总体上没有错误,但我做了一些更改,因为您在循环中执行了不必要的重复(我将外循环移到了内部,"vectorized" 因为所有时刻的计算都可以同时执行给定直方图。此外,它正在构建花费时间最长的直方图)。
intel = cumsum(randn(64,1)); % <-- mock random walk
Ni =length(intel);
%figure, plot(intel)
lower = 1; %set lowest level of histogram resolution (bin size)
upper = 300; %set highest level of histogram resolution (bin size)
qlow = -20; %set lowest moment exponent
qhigh = 20; %set highet moment exponent
qstep = 0.25; %set step size between moment exponents
qn= ((qhigh-qlow)/qstep) + 1; %calculates number of steps given qlow, qhigh, qstep
qvalues= linspace(qlow, qhigh, qn); %creates a vector of q values given above parameters
m = min(intel); %find the maximum of the dataset
M = max(intel); %find the minimum of the dataset
partitions = zeros(upper-lower+1,length(qvalues));
for k = lower:upper %loop over bin sizes
% (1) Select a bin size r and partition [m,M] into intervals of size r:
% [m, m+r), [m+r, m+2r), ..., [m+kr, M], where m+kr < M <= m+(k+1)r.
% Call these bins B0, ..., Bk.
edges = linspace(m,M,k+1);
edges(end)=Inf;
% (2) For each j, 0 <= j <= k, count the number of xi that lie in bin Bj. Call this number nj. Ignore all nj that equal 0 after all the xi have been counted..
counts = histc(intel, edges); %unpack all k histogram height values into 'counts'
counts(counts==0) = []; %delete all zero values in ''counts
% (3) Now compute the qth moment, Mrq = (n0/N)q + ... + (nk/N)q, where the sum is over all nonzero ni.
% Zq = counts/Ni;
partitions(k, :) = sum( (counts/Ni) .^ qvalues); %store Zq in the kth row and the Qth column of 'partitions'
end
figure, hold on
loglog(1./[1:k]', partitions(:,1),'g.-')
loglog(1./[1:k]', partitions(:,80),'b.-')
loglog(1./[1:k]', partitions(:,160),'r.-')
% (4) Perform linear regressions here to get alpha(r) ....
我正在使用耶鲁分形页面中的方法:
http://classes.yale.edu/fractals/MultiFractals/Moments/TSMoments/TSMoments.html
这组讲座幻灯片(幻灯片 32)也对此进行了阐述:
http://multiscale.emsl.pnl.gov/docs/multifractal.pdf
这个想法是你得到一个数据集,并通过许多直方图检查它,增加条的数量,即分辨率。一旦分辨率足够高,一些条形就会取零值。在每个阶段,我们取落入每个直方图箱(忽略任何零值箱)的结果数量,将其除以数据集的总大小,并将其提高到幂 q。此操作给出了给定时刻和 bin 大小的 'partition function'。引用上面链接的教程:“提供了对非均匀性的选择性表征 测量,正 q 突出最密集的区域,负 q 突出最平滑的区域。
所以我在 Matlab 中使用直方图函数,遍历 bin 大小,对所有非零 bin 内容求和,等等。但是我的分区函数输出数组只是一堆 1。我看不出出了什么问题,其他人可以吗?
英特尔、思科、苹果和其他公司的数据可在同一耶鲁网站上获得:yale.edu/fractals/MultiFractals/Finf(a)/Finf(a).html
N.B。 intel
指的是我原先用作数据集的intel股价
lower = 1; %set lowest level of histogram resolution (bin size)
upper = 300; %set highest level of histogram resolution (bin size)
qlow = -20; %set lowest moment exponent
qhigh = 20; %set highet moment exponent
qstep = 0.25; %set step size between moment exponents
qn= ((qhigh-qlow)/qstep) + 1; %calculates number of steps given qlow, qhigh, qstep
qvalues= linspace(qlow, qhigh, qn); %creates a vector of q values given above parameters
m = min(intel); %find the maximum of the dataset
M = max(intel); %find the minimum of the dataset
for Q = 1:length(qvalues) %loop over moment exponents q
for k = lower:upper %loop over bin sizes
counts = hist(intel, k); %unpack all k histogram height values into 'counts'
counts(counts==0) = []; %delete all zero values in ''counts
Zq = counts ./ length(intel);
Zq = Zq .^ qvalues(Q);
Zq = sum(Zq);
partitions(k-(lower-1), Q) = Zq; %store Zq in the kth row and the Qth column of 'partitions'
end
end
您的代码似乎总体上没有错误,但我做了一些更改,因为您在循环中执行了不必要的重复(我将外循环移到了内部,"vectorized" 因为所有时刻的计算都可以同时执行给定直方图。此外,它正在构建花费时间最长的直方图)。
intel = cumsum(randn(64,1)); % <-- mock random walk
Ni =length(intel);
%figure, plot(intel)
lower = 1; %set lowest level of histogram resolution (bin size)
upper = 300; %set highest level of histogram resolution (bin size)
qlow = -20; %set lowest moment exponent
qhigh = 20; %set highet moment exponent
qstep = 0.25; %set step size between moment exponents
qn= ((qhigh-qlow)/qstep) + 1; %calculates number of steps given qlow, qhigh, qstep
qvalues= linspace(qlow, qhigh, qn); %creates a vector of q values given above parameters
m = min(intel); %find the maximum of the dataset
M = max(intel); %find the minimum of the dataset
partitions = zeros(upper-lower+1,length(qvalues));
for k = lower:upper %loop over bin sizes
% (1) Select a bin size r and partition [m,M] into intervals of size r:
% [m, m+r), [m+r, m+2r), ..., [m+kr, M], where m+kr < M <= m+(k+1)r.
% Call these bins B0, ..., Bk.
edges = linspace(m,M,k+1);
edges(end)=Inf;
% (2) For each j, 0 <= j <= k, count the number of xi that lie in bin Bj. Call this number nj. Ignore all nj that equal 0 after all the xi have been counted..
counts = histc(intel, edges); %unpack all k histogram height values into 'counts'
counts(counts==0) = []; %delete all zero values in ''counts
% (3) Now compute the qth moment, Mrq = (n0/N)q + ... + (nk/N)q, where the sum is over all nonzero ni.
% Zq = counts/Ni;
partitions(k, :) = sum( (counts/Ni) .^ qvalues); %store Zq in the kth row and the Qth column of 'partitions'
end
figure, hold on
loglog(1./[1:k]', partitions(:,1),'g.-')
loglog(1./[1:k]', partitions(:,80),'b.-')
loglog(1./[1:k]', partitions(:,160),'r.-')
% (4) Perform linear regressions here to get alpha(r) ....