周期的な成分の有意性検定

この例では、フィッシャー g 統計量を使用してホワイト ノイズを伴う正弦波成分の有意性を評価する方法を示します。フィッシャー g 統計量は、周波数範囲 (0, Fs/2) の 1/2 にわたるピリオドグラムのすべての値に対する、最大ピリオドグラム値の比です。g 統計量と正確な分布の詳細については、参考文献に記載されています。

平均 0、分散 1 のホワイト ガウス ノイズを伴う 100 Hz 正弦波から構成される信号を作成します。正弦波の振幅は 0.25 です。サンプル レートは 1 kHz です。再現性のある結果を得るために、乱数発生器を既定の状態に設定します。

rng default

Fs = 1e3;
t = 0:1/Fs:1-1/Fs;
x = 0.25*cos(2*pi*100*t)+randn(size(t));

periodogram を使用して信号のピリオドグラムを求めます。0 とナイキスト周波数 (Fs/2) は除外します。ピリオドグラムをプロットします。

[Pxx,F] = periodogram(x,rectwin(length(x)),length(x),Fs);
Pxx = Pxx(2:length(x)/2);

periodogram(x,rectwin(length(x)),length(x),Fs)

ピリオドグラムの最大値を検索します。フィッシャーの g 統計量は、全ピリオドグラムの合計に対する最大ピリオドグラム値の比です。

[maxval,index] = max(Pxx);
fisher_g = Pxx(index)/sum(Pxx)

fisher_g = 0.0381

ピリオドグラムの最大値は 100 Hz で発生しますが、それはピリオドグラムの最大値のインデックスに対応する周波数を検索して確認できます。

F = F(2:end-1);
F(index)

ans = 100

参考文献で詳しく述べられている分布結果を使用して、フィッシャー g 統計量の有意水準 pval を決定します。以下の MATLAB® コードは、[2] の方程式 (6) を実装しています。二項係数を計算するときは、ガンマ関数の対数を使用してオーバーフローを回避します。

N = length(Pxx);
nn = 1:floor(1/fisher_g);

I = (-1).^(nn-1).*exp(gammaln(N+1)-gammaln(nn+1)-gammaln(N-nn+1)).*(1-nn*fisher_g).^(N-1);

pval = sum(I)

pval = 2.0163e-06

p 値は 0.00001 よりも小さく、このことは 100 Hz における周期的な成分が有意であることを示しています。フィッシャー g 統計量の解釈は、他の周期性の存在によって複雑になります。複数の周期性がある場合の変更については [1] を参照してください。

参考文献

[1] Percival, Donald B. and Andrew T. Walden. Spectral Analysis for Physical Applications. Cambridge, UK: Cambridge University Press, 1993.

[2] Wichert, Sofia, Konstantinos Fokianos, and Korbinian Strimmer. "Identifying Periodically Expressed Transcripts in Microarray Time Series Data." Bioinformatics. Vol. 20, 2004, pp. 5-20.

参考

nchoosek | periodogram