長時間相関（long-range correlation）や 1/f スケーリングの解析では、トレンドを含む非定常時系列から「本来のゆらぎ」を取り出すために、Detrended Fluctuation Analysis（DFA） や Detrended Moving Average Algorithm (DMA) のような、トレンドからの偏差の RMS（root-mean-square）をスケールごとに評価する方法が広く用いられます。
そのような方法の中で、Savitzky-Golayトレンド除去フィルタを用いたDMA は、単純移動平均を使う原型 DMA に対して性能が良いことが知られており、さらに 高次 DMA（higher-order DMA） を用いることで、より高次の多項式トレンドまで除去できるようになります。とはいえ、これらを素朴に実装すると計算量が非常に大きくなり、データが長い場合には計算時間が増大して実用性が低下する、という問題があります。
そのような問題を解決するため、私たちは 高次 DMA を高速に計算するアルゴリズムを開発しました。本記事では、その高速アルゴリズムを R で実装してみます (Pythonのコードも掲載しました)。以下では、その実装の本体である関数 fastDMA と、推定した からスケーリング指数を求めて図示する関数
estimate_scaling_fastDMA の使い方を解説します。この２つの関数は、この記事の最後の節「[重要] DMAに用いる関数」に掲載してありますので、最初に、そのRスクリプトを実行するのを忘れないでください。
とはいえ、時系列の長さが10万点を超えるような場合は、RやPythonで、DMAを実行するのはお勧めしません。この際、C言語プログラムの使い方、少なくともコンパイルする方法を勉強してみてください。その方が、計算が１万倍速いです。
1. DMA考え方
そもそものDMAの考え方は、別の記事を参考にしてください。
長時間相関解析では，DFAとDMAのどっちが良い？ - ケィオスの時系列解析メモランダム
0次DMA（DMA0）
最も基本の DMA で、トレンド推定を単純移動平均で行い、積分時系列から推定されたトレンドを除きます。
2次以上の高次DMA
高次 DMA では、窓内で 多項式（moving average polynomial） を最小二乗で当てはめることで、トレンドを推定します。この多項式フィットでトレンドを推定する方法は、Savitzy-Golayフィルタと呼ばれるものと同じものです。
2. なぜ高速になるのか
高次 DMA をそのまま実装すると、各スケールで各窓ごとに多項式当てはめ（＝多くの和の計算）を繰り返すため、計算量が急増します。そこで、和の再帰更新（recurrence formulas）を使って計算量を大幅に低減し計算を高速化します。
和の再帰更新（recurrence formulas）とは、トレンド推定に必要な部分区間ごとの和を、部分区間を 1 点ずらしたときに、
「前の和」＋「入ってきた点」−「出ていった点」
という形で更新することです。高速アルゴリズムでは、この再帰式を時間方向とスケール方向の両方に展開しています。たとえば、DMA0で用いる単純な移動平均では、窓幅を とし、時系列を
と書くと、位置
における窓内平均は
で定義されます。ここで、窓を 1 点だけ右にずらした における和を考えると、
となりますが、この和は新たにすべてを足し直さなくても、次のように書き換えられます。
DMA0 の高速化は、まさにこの 「前の和」＋「入ってきた点」−「出ていった点」という再帰更新を用いて、窓をスライドさせながら効率よくトレンドを計算している点にあります。
高次 DMA（DMA2, DMA4）では、この考え方を拡張し、単純な和だけでなく、その他のモーメントについても同様の再帰更新式を導入します。
ただし、このような再帰更新を延々と繰り返すと計算結果に丸め誤差が蓄積するため、一定間隔で refresh interval を置き、そこで sum() による再計算を挟んで誤差を抑えます。
3. とりあえず使ってみる
細かい理屈はすっ飛ばし、動かしてみましょう。fastDMA は各スケールにおけるゆらぎ関数 （正確には
とそこから得られる
）を計算し、
estimate_scaling_fastDMA はその結果から 指定したスケール範囲 ( におけるスケーリング指数（log–log 直線の傾き）を推定して図示します。
以下のスクリプトは、最小構成の実行例です。
# 例: 長時間相関を持つデータ（外部パッケージを使う場合）
# ここでは時系列データを x に格納
# 自分のデータを分析する場合は、以下の x を自分のデータに入れ替えてください
N <- 10000
x <- rnorm(N)
# DMA2の実行例：スケーリング指数推定とプロット
res_fit <- estimate_scaling_fastDMA(
x, s1 = 15, s2 = length(x)/20,
order = 2,
scale_correction = TRUE
)
# 推定されたスケーリング指数
res_fit$slope
このスクリプトを実行すれば、以下のようなプロットがえられます。
■ estimate_scaling_fastDMA() の使い方
上のスクリプトで実行した
res_fit <- estimate_scaling_fastDMA( x, s1 = 15, s2 = length(x)/20, order = 2, scale_correction = TRUE )
の部分では、内部で fastDMA() を呼び出し、次の処理をまとめて行います。
を計算
–
をプロット
- 指定範囲
に入る点だけを使って直線回帰
- 回帰直線（破線）と、範囲
- 傾き（スケーリング指数）、
、使用点数、実行時間を図中に表示
ここで指定しているパラメタは次のとおりです。
s1 = 15小スケール側の下限です。小さすぎるスケールは離散化やフィルタ形状の影響が出やすいので、ある程度大きい値から始めます。s2 = length(x)/20大スケール側の上限です。スケールが大きすぎると評価点数が減り、不安定になりやすいので「データ長の 1/20 程度」に抑えています（この係数は目的に応じて調整してください）。scale_correction = TRUEorder に応じて横軸スケールに補正係数を入れるオプションです。 ただし **(s_1, s_2) 自体は補正しません。
推定された傾きを数値として取り出したい場合は
res_fit$slope
を参照してください。あわせて res_fit$r2（決定係数）や res_fit$elapsed_sec（実行時間）も取得できます。
参考文献
- Tsujimoto, Y. Miki, S. Shimatani, and K. Kiyono, “Fast algorithm for scaling analysis with higher-order detrending moving average method,” Physical Review E 93, 053304 (2016).
[重要] DMAに用いる関数
以上の記事に説明にあるコマンドを使うために、R起動後に以下のスクリプトを一度実行してください。
# ------------------------------------------------------------
# fastDMA: Unified Centered DMA (order = 0 / 2 / 4)
#
# 目的:
# Centered DMA（Detrending Moving Average）により、スケールごとの
# ゆらぎ関数 F2（平方平均偏差）を高速アルゴリズムで計算します。
#
# 特徴:
# - order は 0, 2, 4 のみ（局所近似の次数に対応）
# ------------------------------------------------------------
# ---------- scale generator (odd scales) ----------
# minbox から maxbox までを幾何級数的に増やしてスケール列を作る。
# Centered DMA では窓長が奇数であることが自然なので、必ず奇数に丸める。
fastDMA_rscale <- function(minbox, maxbox, boxratio) {
# 予定されるスケール数（log 比で概算）
rslen <- floor(log(maxbox / minbox, base = 10) / log(boxratio, base = 10) + 1)
if (rslen < 1) return(as.integer(minbox))
rs <- integer(rslen)
rs[1] <- as.integer(minbox)
idx <- 2L
ir <- 1L
# boxratio でスケールを増加させつつ、単調増加になるものだけ採用する
while (idx <= rslen && rs[idx - 1L] < maxbox) {
rw <- floor(minbox * (boxratio^ir) + 0.5) # 四捨五入
if (rw > rs[idx - 1L]) {
# 奇数化（(rw %/% 2)*2 で偶数、+1 で奇数）
rs[idx] <- (rw %/% 2L) * 2L + 1L
idx <- idx + 1L
}
ir <- ir + 1L
}
# 末尾の調整（maxbox を超えていたら除外）
idx <- idx - 1L
if (rs[idx] > maxbox) idx <- idx - 1L
rs[seq_len(idx)]
}
# ---------- order 0: one-scale kernel ----------
# order 0（局所平均）で、単一 scale の F2 を計算する。
#
# 考え方:
# scale の移動窓で局所平均 a0 を計算し、窓中心点 y[ic] との差 d を取り、
# d^2 を全窓位置で合計して f2 とする。
#
# 注意:
# - Centered DMA では、中心が定義できる点だけを評価する。
# - 評価点数は (n - scale + 1) = itemp + 1。
#
# 高速化:
# - y10 = 窓内の総和をローリング更新（+新規 - 退出）
# - i_refresh 間隔で sum() により再計算し、浮動小数の累積誤差を抑える
fastDMA0_F2_one <- function(y, scale, i_refresh) {
n <- length(y)
# scale = 2k+1（奇数）とすると、窓中心は k+1 番目
k <- (scale - 1L) %/% 2L
ic <- k + 1L
# 窓を 1 つ進められる回数（有効な中心が取れる範囲）
itemp <- n - scale
# 初期窓 [1:scale] の総和
y10 <- sum(y[1L:scale])
c1 <- 1 / scale
# 初期中心での平均との差
a0 <- c1 * y10
d <- y[ic] - a0
f2 <- d * d
if (itemp > 0L) {
# リフレッシュをしない（または全区間より長い）場合は単純ローリング
if (i_refresh <= 0L || i_refresh > itemp) {
for (i in 1L:itemp) {
y10 <- y10 + y[i + scale] - y[i]
a0 <- c1 * y10
d <- y[i + ic] - a0
f2 <- f2 + d * d
}
} else {
# i_refresh ごとに sum() を取り直して数値誤差を抑制
i <- 1L
while (i <= itemp) {
next_ref <- ((i %/% i_refresh) + 1L) * i_refresh
if (next_ref > itemp) next_ref <- itemp
# next_ref 手前まではローリング更新
if (i <= (next_ref - 1L)) {
for (j in i:(next_ref - 1L)) {
y10 <- y10 + y[j + scale] - y[j]
a0 <- c1 * y10
d <- y[j + ic] - a0
f2 <- f2 + d * d
}
}
# next_ref の位置で窓総和を再計算
if (next_ref <= itemp) {
start <- next_ref + 1L
end <- next_ref + scale
y10 <- sum(y[start:end])
a0 <- c1 * y10
d <- y[next_ref + ic] - a0
f2 <- f2 + d * d
}
i <- next_ref + 1L
}
}
}
# 実際に評価した中心点数 = itemp + 1 = (n - scale + 1)
f2 / as.numeric(itemp + 1L)
}
# ---------- order 2 coefficients ----------
# order 2（局所 2 次）で用いる係数 c1, c2 を返す。
#
# 実装上は、窓内の以下の和を使って局所トレンド a0 を作る：
# y10 = Σ y
# y1ii = Σ y * d^2 （d は中心からの相対インデックス）
# a0 = c1*y10 + c2*y1ii
fastDMA2_cmat <- function(k) {
k <- as.numeric(k)
k2 <- k * k
den <- 8 * k2 * k + 12 * k2 - 2 * k - 3
c1 <- (9 * k2 + 9 * k - 3) / den
c2 <- -15 / den
c(c1 = c1, c2 = c2)
}
# ---------- order 4 coefficients ----------
# order 4（局所 4 次）で用いる係数 c1, c2, c3 を返す。
# 同様に、窓内の多項式モーメント（d^0, d^2, d^4）を使い、
# a0 = c1*y10 + c2*y1ii + c3*y1iv
# を構成する。
fastDMA4_cmat <- function(k) {
k <- as.numeric(k)
k2 <- k * k
k3 <- k2 * k
k4 <- k3 * k
den <- 180 + 72 * k - 800 * k2 - 320 * k3 + 320 * k4 + 128 * k * k4
c1 <- (180 - 750 * k - 525 * k2 + 450 * k3 + 225 * k4) / den
c2 <- (1575 - 1050 * k - 1050 * k2) / den
c3 <- 945 / den
c(c1 = c1, c2 = c2, c3 = c3)
}
# ---------- order 2: all scales, cache-based, in-place ----------
# order 2 の F2 を、複数スケールまとめて計算する。
#
# 高速化の要点:
# 1) あるスケールで用いた「窓内の和（モーメント）」を i_refresh ごとにキャッシュする
# 2) 次のスケールでは、前スケールのキャッシュから初期値を更新して再利用する
# （prev_scale -> scale に変わっても中心基準のモーメント変換で流用できる）
#
# 分母:
# - 実際に評価した中心点数は (n - scale + 1) = itemp + 1
fastDMA2_F2_all <- function(y, scales, i_refresh) {
n <- length(y)
nr <- length(scales)
i_refresh <- as.integer(i_refresh)
if (i_refresh < 1L) i_refresh <- 1L
if (i_refresh > n) i_refresh <- n
# キャッシュは「リフレッシュ地点」ごとに保存する
cache_len <- as.integer(n %/% i_refresh + 1L)
loc_sum0 <- numeric(cache_len) # Σ y
loc_sum1 <- numeric(cache_len) # Σ y*d
loc_sum2 <- numeric(cache_len) # Σ y*d^2
F2 <- numeric(nr)
prev_scale <- 0L
for (k2 in seq_len(nr)) {
scale <- scales[k2]
k <- (scale - 1L) %/% 2L
ic <- k + 1L
itemp <- n - scale
# このスケールでの係数
cc <- fastDMA2_cmat(k)
c1 <- cc[["c1"]]
c2 <- cc[["c2"]]
# iloc は「何番目のキャッシュ地点か」
iloc <- 1L
# 初期窓（j=0 相当）でのモーメント（Σ y, Σ y*d, Σ y*d^2）を準備
if (prev_scale == 0L) {
i0 <- 1L
y10 <- 0.0; y1i <- 0.0; y1ii <- 0.0
nn1 <- 0.0; nn2 <- 0.0
} else {
# 以前のスケール（prev_scale）のキャッシュから、現在の中心へ座標変換して流用
kb <- (prev_scale - 1L) %/% 2L
nn1 <- as.numeric(kb - k); nn2 <- nn1 * nn1
i0 <- prev_scale + 1L
y10 <- loc_sum0[iloc]
y1i <- loc_sum1[iloc] + y10 * nn1
y1ii <- loc_sum2[iloc] + 2 * loc_sum1[iloc] * nn1 + y10 * nn2
}
# i0..scale の分だけ、未加算分を足し込んで「今の scale の窓」に合わせる
if (i0 <= scale) {
idx <- i0:scale
di <- as.numeric(idx - k - 1L) # 中心からの相対位置
yi <- y[idx]
y10 <- y10 + sum(yi)
y1i <- y1i + sum(yi * di)
y1ii <- y1ii + sum(yi * di * di)
}
# 初期窓のモーメントをキャッシュ
loc_sum0[iloc] <- y10
loc_sum1[iloc] <- y1i
loc_sum2[iloc] <- y1ii
# 初期位置のトレンド値と偏差
a0 <- c1 * y10 + c2 * y1ii
d <- y[ic] - a0
f2 <- d * d
if (itemp > 0L) {
# 更新式で頻出する定数
k1 <- as.numeric(k + 1L)
kk <- as.numeric(k)
i <- 1L
while (i <= itemp) {
next_ref <- ((i %/% i_refresh) + 1L) * i_refresh
if (next_ref > itemp) next_ref <- itemp
# リフレッシュ直前までは、モーメントを再帰更新して高速に前進
if (i <= (next_ref - 1L)) {
for (j in i:(next_ref - 1L)) {
y00 <- y10
y0i <- y1i
y0ii <- y1ii
# Σ y をローリング更新
y10 <- y00 + y[j + scale] - y[j]
# Σ y*d と Σ y*d^2 を更新（窓を 1 つずらしたときの再帰式）
temp1 <- y[j] * k1
temp2 <- y[j + scale] * kk
y1i <- y0i - y00 + temp1 + temp2
y1ii <- y0ii - 2.0 * y0i + y00 - temp1 * k1 + temp2 * kk
a0 <- c1 * y10 + c2 * y1ii
d <- y[j + ic] - a0
f2 <- f2 + d * d
}
}
# next_ref の位置でモーメントを再構成（キャッシュ更新点）
if (next_ref <= itemp) {
iloc <- iloc + 1L
if (prev_scale == 0L) {
i0 <- 1L
y10 <- 0.0; y1i <- 0.0; y1ii <- 0.0
} else {
i0 <- prev_scale + 1L
y10 <- loc_sum0[iloc]
y1i <- loc_sum1[iloc] + y10 * nn1
y1ii <- loc_sum2[iloc] + 2 * loc_sum1[iloc] * nn1 + y10 * nn2
}
if (i0 <= scale) {
jj <- i0:scale
idx_y <- next_ref + jj
di <- as.numeric(jj - k - 1L)
yi <- y[idx_y]
y10 <- y10 + sum(yi)
y1i <- y1i + sum(yi * di)
y1ii <- y1ii + sum(yi * di * di)
}
loc_sum0[iloc] <- y10
loc_sum1[iloc] <- y1i
loc_sum2[iloc] <- y1ii
a0 <- c1 * y10 + c2 * y1ii
d <- y[next_ref + ic] - a0
f2 <- f2 + d * d
}
i <- next_ref + 1L
}
}
# 実際に評価した中心点数 = itemp + 1 = (n - scale + 1)
F2[k2] <- f2 / as.numeric(itemp + 1L)
prev_scale <- scale
}
F2
}
# ---------- order 4: all scales, cache-based, in-place ----------
# order 4 の F2 を、複数スケールまとめて計算する（order 2 の拡張）。
#
# 追加で扱うモーメント:
# y1iii = Σ y*d^3
# y1iv = Σ y*d^4
#
# a0 は (d^0, d^2, d^4) の組で構成し、中心点との差の二乗を合計する。
#
# 分母:
# - 実際に評価した中心点数は (n - scale + 1) = itemp + 1
fastDMA4_F2_all <- function(y, scales, i_refresh) {
n <- length(y)
nr <- length(scales)
i_refresh <- as.integer(i_refresh)
if (i_refresh < 1L) i_refresh <- 1L
if (i_refresh > n) i_refresh <- n
cache_len <- as.integer(n %/% i_refresh + 1L)
loc_sum0 <- numeric(cache_len) # Σ y
loc_sum1 <- numeric(cache_len) # Σ y*d
loc_sum2 <- numeric(cache_len) # Σ y*d^2
loc_sum3 <- numeric(cache_len) # Σ y*d^3
loc_sum4 <- numeric(cache_len) # Σ y*d^4
F2 <- numeric(nr)
prev_scale <- 0L
for (k2 in seq_len(nr)) {
scale <- scales[k2]
k <- (scale - 1L) %/% 2L
ic <- k + 1L
itemp <- n - scale
cc <- fastDMA4_cmat(k)
c1 <- cc[["c1"]]; c2 <- cc[["c2"]]; c3 <- cc[["c3"]]
iloc <- 1L
# 初期窓のモーメント準備（prev_scale があればキャッシュから変換して流用）
if (prev_scale == 0L) {
i0 <- 1L
y10 <- 0.0; y1i <- 0.0; y1ii <- 0.0; y1iii <- 0.0; y1iv <- 0.0
nn1 <- 0.0; nn2 <- 0.0; nn3 <- 0.0; nn4 <- 0.0
} else {
kb <- (prev_scale - 1L) %/% 2L
nn1 <- as.numeric(kb - k)
nn2 <- nn1 * nn1
nn3 <- nn2 * nn1
nn4 <- nn2 * nn2
i0 <- prev_scale + 1L
y10 <- loc_sum0[iloc]
y1i <- loc_sum1[iloc] + y10 * nn1
y1ii <- loc_sum2[iloc] + 2 * loc_sum1[iloc] * nn1 + y10 * nn2
y1iii <- loc_sum3[iloc] + 3 * loc_sum2[iloc] * nn1 + 3 * loc_sum1[iloc] * nn2 + y10 * nn3
y1iv <- loc_sum4[iloc] + 4 * loc_sum3[iloc] * nn1 + 6 * loc_sum2[iloc] * nn2 +
4 * loc_sum1[iloc] * nn3 + y10 * nn4
}
# 未加算分を足し込み、現在 scale の窓に整合させる
if (i0 <= scale) {
idx <- i0:scale
di <- as.numeric(idx - k - 1L)
yi <- y[idx]
di2 <- di * di
y10 <- y10 + sum(yi)
y1i <- y1i + sum(yi * di)
y1ii <- y1ii + sum(yi * di2)
y1iii <- y1iii + sum(yi * di2 * di)
y1iv <- y1iv + sum(yi * di2 * di2)
}
# 初期窓モーメントをキャッシュ
loc_sum0[iloc] <- y10
loc_sum1[iloc] <- y1i
loc_sum2[iloc] <- y1ii
loc_sum3[iloc] <- y1iii
loc_sum4[iloc] <- y1iv
# 初期中心でのトレンド推定と偏差
a0 <- c1 * y10 + c2 * y1ii + c3 * y1iv
d <- y[ic] - a0
f2 <- d * d
if (itemp > 0L) {
k1 <- as.numeric(k + 1L)
kk <- as.numeric(k)
i <- 1L
while (i <= itemp) {
next_ref <- ((i %/% i_refresh) + 1L) * i_refresh
if (next_ref > itemp) next_ref <- itemp
# リフレッシュ直前までは再帰更新で前進
if (i <= (next_ref - 1L)) {
for (j in i:(next_ref - 1L)) {
y00 <- y10
y0i <- y1i
y0ii <- y1ii
y0iii <- y1iii
y0iv <- y1iv
y10 <- y00 + y[j + scale] - y[j]
# 以下は窓を 1 つずらしたときのモーメント更新（order 4 用に拡張）
t1 <- y[j] * k1
t2 <- y[j + scale] * kk
y1i <- y0i - y00 + t1 + t2
t1 <- t1 * k1
t2 <- t2 * kk
y1ii <- y0ii - 2 * y0i + y00 - t1 + t2
t1 <- t1 * k1
t2 <- t2 * kk
y1iii <- y0iii - 3 * y0ii + 3 * y0i - y00 + t1 + t2
y1iv <- y0iv - 4 * y0iii + 6 * y0ii - 4 * y0i + y00 - t1 * k1 + t2 * kk
a0 <- c1 * y10 + c2 * y1ii + c3 * y1iv
d <- y[j + ic] - a0
f2 <- f2 + d * d
}
}
# next_ref の位置でモーメントを再構成（キャッシュ更新点）
if (next_ref <= itemp) {
iloc <- iloc + 1L
if (prev_scale == 0L) {
i0 <- 1L
y10 <- 0.0; y1i <- 0.0; y1ii <- 0.0; y1iii <- 0.0; y1iv <- 0.0
} else {
i0 <- prev_scale + 1L
y10 <- loc_sum0[iloc]
y1i <- loc_sum1[iloc] + y10 * nn1
y1ii <- loc_sum2[iloc] + 2 * loc_sum1[iloc] * nn1 + y10 * nn2
y1iii <- loc_sum3[iloc] + 3 * loc_sum2[iloc] * nn1 + 3 * loc_sum1[iloc] * nn2 + y10 * nn3
y1iv <- loc_sum4[iloc] + 4 * loc_sum3[iloc] * nn1 + 6 * loc_sum2[iloc] * nn2 +
4 * loc_sum1[iloc] * nn3 + y10 * nn4
}
if (i0 <= scale) {
jj <- i0:scale
idx_y <- next_ref + jj
di <- as.numeric(jj - k - 1L)
yi <- y[idx_y]
di2 <- di * di
y10 <- y10 + sum(yi)
y1i <- y1i + sum(yi * di)
y1ii <- y1ii + sum(yi * di2)
y1iii <- y1iii + sum(yi * di2 * di)
y1iv <- y1iv + sum(yi * di2 * di2)
}
loc_sum0[iloc] <- y10
loc_sum1[iloc] <- y1i
loc_sum2[iloc] <- y1ii
loc_sum3[iloc] <- y1iii
loc_sum4[iloc] <- y1iv
a0 <- c1 * y10 + c2 * y1ii + c3 * y1iv
d <- y[next_ref + ic] - a0
f2 <- f2 + d * d
}
i <- next_ref + 1L
}
}
# 実際に評価した中心点数 = itemp + 1 = (n - scale + 1)
F2[k2] <- f2 / as.numeric(itemp + 1L)
prev_scale <- scale
}
F2
}
# ------------------------------------------------------------
# fastDMA: ユーザ向け関数
#
# 入力:
# x : 時系列（numeric vector）
# order : 0 / 2 / 4（局所トレンド近似の次数に対応）
# minbox : 最小スケール（窓長）。NULL の場合は order に応じて自動設定
# maxbox : 最大スケール。0 以下なら n/4 を採用
# boxratio : スケール増加倍率（幾何級数の比）
# integrate : TRUE のとき、x の平均を引いてから累積和（DFA と同様の積分ステップ）
# i_refresh : ローリング更新の再計算間隔（数値誤差を抑えるためのリフレッシュ）
#
# 出力:
# list(order, scales, F2, log10_scale, log10_F)
# - log10_F は sqrt(F2) の log10（= 0.5*log10(F2)）
# ------------------------------------------------------------
fastDMA <- function(x,
order = 0L,
minbox = NULL,
maxbox = 0L,
boxratio = 2^(3/16),
integrate = TRUE,
i_refresh = NULL) {
stopifnot(is.numeric(x), length(x) >= 10)
n <- length(x)
order <- as.integer(order)
if (!order %in% c(0L, 2L, 4L)) stop("order must be one of 0, 2, 4.")
# order に応じた、実用的な最小窓長のデフォルト
if (is.null(minbox)) {
minbox <- if (order == 0L) 5L else if (order == 2L) 7L else 9L
}
# order に応じて i_refresh のデフォルトを調整（order が高いほど更新式が重いので短めに）
if (is.null(i_refresh)) {
i_refresh <- if (order == 0L) 1000L else if (order == 2L) 200L else 20L
}
minbox <- as.integer(minbox)
maxbox <- as.integer(maxbox)
i_refresh <- as.integer(i_refresh)
# ---- 積分ステップ（任意）----
# 長時間相関解析では、平均除去してから累積和を取った系列 y を解析対象にすることが多い。
if (integrate) {
x0 <- x - mean(x, na.rm = TRUE)
y <- cumsum(x0)
} else {
y <- as.numeric(x)
}
# ---- 窓長の整形 ----
# Centered DMA では窓長は奇数が自然なので奇数化し、
# order に応じて最低限必要な窓長を下回らないようにする。
min_required <- if (order == 0L) 5L else if (order == 2L) 7L else 9L
minbox <- (minbox %/% 2L) * 2L + 1L
if (minbox < min_required) minbox <- min_required
# ---- 最大窓長の制約 ----
# スケールが大きすぎると推定点が少なくなり不安定になるため、ここでは n/4 を上限とする。
if (maxbox <= 0L || maxbox > (n %/% 4L) || minbox > maxbox) {
maxbox <- as.integer(n %/% 4L)
}
# 使用するスケール列を作成
scales <- fastDMA_rscale(minbox, maxbox, boxratio)
# ---- F2 の計算 ----
if (order == 0L) {
# order 0 はスケールごとに独立な簡潔カーネルなので vapply で回す
F2 <- vapply(scales, function(s) fastDMA0_F2_one(y, s, i_refresh), numeric(1L))
} else if (order == 2L) {
# order 2 はキャッシュを用いて全スケールを一括計算
F2 <- fastDMA2_F2_all(y, scales, i_refresh)
} else {
# order 4 も同様に一括計算
F2 <- fastDMA4_F2_all(y, scales, i_refresh)
}
# ---- 返り値 ----
list(
order = order,
scales = scales,
F2 = F2,
log10_scale = log10(scales),
log10_F = 0.5 * log10(F2)
)
}
# =========================================================
# fastDMA scaling analysis with automatic order correction
# + optional scale correction
#
# 目的:
# fastDMA によって得られた F(s) に対し，
# 指定したスケール区間 [s1, s2] でスケーリング指数を推定する．
#
# 主な機能:
# - 偶数次数（0,2,4）のみをサポート（奇数指定は自動補正）
# - スケール補正オプション
# scale_correction == TRUE のとき
# order 0,1 -> s / 1
# order 2,3 -> s / 1.93
# order 4,5 -> s / 2.74
# - log–log プロットと回帰直線を自動描画
# - 実行時間を計測・表示
#
# 描画仕様:
# - 点: 青（col = 4）
# - 回帰直線: 破線
# - s1, s2 の位置に縦破線（視覚的なガイド）
# =========================================================
estimate_scaling_fastDMA <- function(x, s1, s2,
order = 0L,
minbox = NULL,
maxbox = 0L,
boxratio = 2^(3/16),
integrate = TRUE,
i_refresh = NULL,
scale_correction = TRUE,
main = NULL,
xlab = expression(log[10](s)),
ylab = expression(log[10](F(s))),
pch = 16, cex = 0.8) {
# ----------------------------------------------------------
# 実行時間計測（開始）
# ----------------------------------------------------------
t_start <- proc.time()
# 入力の基本チェック
x <- as.numeric(x)
n <- length(x)
stopifnot(n >= 10)
# ----------------------------------------------------------
# order の自動補正
# fastDMA は偶数次数（0,2,4）のみ対応しているため，
# 奇数が指定された場合は直近の偶数に丸める
# ----------------------------------------------------------
order_input <- as.integer(order)
order_used <- switch(as.character(order_input),
"1" = 0L,
"3" = 2L,
"5" = 4L,
order_input)
if (!identical(order_input, order_used)) {
message(sprintf(
"Note: order=%d was replaced by order=%d (nearest supported even order).",
order_input, order_used
))
}
if (!order_used %in% c(0L, 2L, 4L)) {
stop("order must be one of 0, 1, 2, 3, 4, 5 (odd values are auto-corrected).")
}
# ----------------------------------------------------------
# スケーリング区間 [s1, s2] のチェック
# ここでは「補正前スケール」として扱う
# ----------------------------------------------------------
if (!is.finite(s1) || !is.finite(s2) || s1 <= 0 || s2 <= 0) {
stop("s1 and s2 must be positive finite numbers.")
}
if (s1 > s2) { tmp <- s1; s1 <- s2; s2 <- tmp }
# ----------------------------------------------------------
# fastDMA の実行
# ----------------------------------------------------------
out <- fastDMA(
x = x,
order = order_used,
minbox = minbox,
maxbox = as.integer(maxbox),
boxratio = boxratio,
integrate = integrate,
i_refresh = i_refresh
)
# fastDMA の出力
s_raw <- out$scales # 生のスケール
F <- sqrt(out$F2) # F(s) = sqrt(F2)
# ----------------------------------------------------------
# スケール補正（x 軸表示・回帰用）
# s1, s2 自体は補正しない
# ----------------------------------------------------------
corr_factor <- 1.0
if (isTRUE(scale_correction)) {
corr_factor <- if (order_input %in% c(0L, 1L)) {
1.0
} else if (order_input %in% c(2L, 3L)) {
1.93
} else if (order_input %in% c(4L, 5L)) {
2.74
} else {
1.0
}
if (!is.finite(corr_factor) || corr_factor <= 0) corr_factor <- 1.0
message(sprintf("Scale correction enabled: s <- s / %.3g (s1,s2 unchanged)", corr_factor))
}
# 補正後スケール（プロット・回帰用）
s <- s_raw / corr_factor
log10_s <- log10(s)
log10_F <- log10(F)
# ----------------------------------------------------------
# 回帰に使用する点の選択
# s1, s2 は補正前スケールで判定する
# ----------------------------------------------------------
ok <- is.finite(log10_s) &
is.finite(log10_F) &
(s_raw >= s1) & (s_raw <= s2)
if (sum(ok) < 2) stop("Not enough points within [s1, s2] to fit a line.")
# ----------------------------------------------------------
# log–log 空間での線形回帰
# ----------------------------------------------------------
fit <- lm(log10_F[ok] ~ log10_s[ok])
slope <- unname(coef(fit)[2])
intercept <- unname(coef(fit)[1])
r2 <- summary(fit)$r.squared
# ----------------------------------------------------------
# 実行時間計測（終了）
# ----------------------------------------------------------
t_end <- proc.time()
elapsed_sec <- as.numeric((t_end - t_start)[["elapsed"]])
# ----------------------------------------------------------
# プロット
# ----------------------------------------------------------
if (is.null(main)) {
if (isTRUE(scale_correction)) {
main <- sprintf("fastDMA (order=%d): fit in [%g, %g] (s/%.3g)",
order_used, s1, s2, corr_factor)
} else {
main <- sprintf("fastDMA (order=%d): fit in [%g, %g]", order_used, s1, s2)
}
}
plot(log10_s, log10_F,
pch = pch, cex = cex, col = 4,
xlab = xlab, ylab = ylab, main = main)
# 回帰に使用した点を強調
points(log10_s[ok], log10_F[ok],
pch = pch, cex = cex + 0.2, col = 4)
# ----------------------------------------------------------
# s1, s2 に対応する縦破線（視覚的ガイド）
# ----------------------------------------------------------
x_v <- log10(c(s1, s2) / corr_factor)
y_center <- approx(log10_s, log10_F, xout = x_v)$y
y_rng <- range(log10_F, na.rm = TRUE)
half_h <- (y_rng[2] - y_rng[1]) / 10
for (j in 1:2) {
segments(x_v[j],
y_center[j] - half_h,
x_v[j],
y_center[j] + half_h,
lty = 2, col = "gray40", lwd = 1.5)
}
# ----------------------------------------------------------
# 回帰直線
# ----------------------------------------------------------
xx <- seq(min(log10_s[ok]), max(log10_s[ok]), length.out = 200)
yy <- intercept + slope * xx
lines(xx, yy, lwd = 2, lty = 2)
# ----------------------------------------------------------
# 注釈（推定結果 + 実行時間）
# ----------------------------------------------------------
txt <- if (isTRUE(scale_correction)) {
sprintf(
"order = %d (input=%d)\nscale: s/%.3g (s1,s2 raw)\nslope = %.5f\nR^2 = %.4f\nn = %d\nelapsed = %.3f s",
order_used, order_input, corr_factor, slope, r2, sum(ok), elapsed_sec
)
} else {
sprintf(
"order = %d (input=%d)\nslope = %.5f\nR^2 = %.4f\nn = %d\nelapsed = %.3f s",
order_used, order_input, slope, r2, sum(ok), elapsed_sec
)
}
usr <- par("usr")
text(usr[1] + 0.02 * diff(usr[1:2]),
usr[4] - 0.02 * diff(usr[3:4]),
labels = txt, adj = c(0, 1))
# ----------------------------------------------------------
# 返り値（解析結果をまとめて返す）
# ----------------------------------------------------------
invisible(list(
out = out,
fit = fit,
slope = slope,
intercept = intercept,
r2 = r2,
n_used = sum(ok),
idx_used = which(ok),
s_range_input = c(s1, s2),
order_input = order_input,
order_used = order_used,
scale_correction = isTRUE(scale_correction),
scale_divisor = corr_factor,
elapsed_sec = elapsed_sec
))
}
[おまけ] Pythonでもやってみる
Pythonの場合は、以下のコードを実行して関数を定義します。
"""
fastDMA: Unified Centered DMA (order = 0 / 2 / 4) + scaling exponent estimation
This is a direct Python port of the provided R implementation, preserving:
- odd scale generation
- integrate option (mean removal + cumulative sum)
- i_refresh "rebuild" strategy to limit floating-point drift
- order=0 single-scale rolling sum kernel
- order=2 and order=4 cache-based multi-scale kernels (as in the R code)
- estimate_scaling_fastDMA(): odd order auto-correction, optional scale correction, log-log fit, plot
Dependencies: numpy, matplotlib
(Optional but used here for regression) numpy.polyfit
"""
from __future__ import annotations
import math
import time
from dataclasses import dataclass
from typing import Dict, Iterable, List, Optional, Tuple, Union
import numpy as np
import matplotlib.pyplot as plt
# =========================================================
# Scale generator (odd scales)
# =========================================================
def fastDMA_rscale(minbox: int, maxbox: int, boxratio: float) -> np.ndarray:
"""
Generate an increasing sequence of odd scales from minbox to maxbox
using geometric progression with ratio boxratio, then rounding and forcing odd.
Mirrors the R function fastDMA_rscale().
"""
minbox = int(minbox)
maxbox = int(maxbox)
if minbox <= 0:
raise ValueError("minbox must be positive.")
if maxbox < minbox:
return np.array([minbox], dtype=int)
# planned length (log estimate)
rslen = int(math.floor(math.log10(maxbox / minbox) / math.log10(boxratio) + 1.0))
if rslen < 1:
return np.array([minbox], dtype=int)
rs = np.zeros(rslen, dtype=int)
rs[0] = minbox
idx = 1
ir = 1
while idx < rslen and rs[idx - 1] < maxbox:
rw = int(math.floor(minbox * (boxratio ** ir) + 0.5)) # round half up-ish
if rw > rs[idx - 1]:
rs[idx] = (rw // 2) * 2 + 1 # force odd
idx += 1
ir += 1
idx -= 1
if idx >= 0 and rs[idx] > maxbox:
idx -= 1
if idx < 0:
return np.array([minbox], dtype=int)
return rs[: idx + 1]
# =========================================================
# Order 0: one-scale kernel
# =========================================================
def fastDMA0_F2_one(y: np.ndarray, scale: int, i_refresh: int) -> float:
"""
Compute F2 for one scale using order-0 centered DMA (local mean).
Direct port of fastDMA0_F2_one().
y: 1D array
scale: odd window length (2k+1)
i_refresh: rebuild interval; if <=0 or > itemp -> pure rolling update
"""
y = np.asarray(y, dtype=float)
n = y.size
scale = int(scale)
if scale < 1 or scale > n:
raise ValueError("scale must be in [1, len(y)].")
if scale % 2 == 0:
raise ValueError("scale must be odd for centered DMA.")
k = (scale - 1) // 2
ic = k # 0-based index of center within window
itemp = n - scale # number of shifts
y10 = float(np.sum(y[0:scale]))
c1 = 1.0 / scale
a0 = c1 * y10
d = y[ic] - a0
f2 = d * d
if itemp > 0:
if i_refresh is None:
i_refresh = 0
i_refresh = int(i_refresh)
if i_refresh <= 0 or i_refresh > itemp:
for i in range(1, itemp + 1):
y10 = y10 + y[i + scale - 1] - y[i - 1]
a0 = c1 * y10
d = y[i + ic] - a0
f2 += d * d
else:
i = 1
while i <= itemp:
next_ref = ((i // i_refresh) + 1) * i_refresh
if next_ref > itemp:
next_ref = itemp
# rolling up to next_ref - 1
if i <= next_ref - 1:
for j in range(i, next_ref):
y10 = y10 + y[j + scale - 1] - y[j - 1]
a0 = c1 * y10
d = y[j + ic] - a0
f2 += d * d
# rebuild at next_ref
if next_ref <= itemp:
start = next_ref
end = next_ref + scale
y10 = float(np.sum(y[start:end]))
a0 = c1 * y10
d = y[next_ref + ic] - a0
f2 += d * d
i = next_ref + 1
return f2 / float(itemp + 1)
# =========================================================
# Order 2 coefficients
# =========================================================
def fastDMA2_cmat(k: int) -> Dict[str, float]:
"""
Return coefficients for order-2 local quadratic centered DMA:
a0 = c1 * sum(y) + c2 * sum(y*d^2)
Direct port of fastDMA2_cmat().
"""
k = float(k)
k2 = k * k
den = 8 * k2 * k + 12 * k2 - 2 * k - 3
c1 = (9 * k2 + 9 * k - 3) / den
c2 = -15.0 / den
return {"c1": c1, "c2": c2}
# =========================================================
# Order 4 coefficients
# =========================================================
def fastDMA4_cmat(k: int) -> Dict[str, float]:
"""
Return coefficients for order-4 local quartic centered DMA:
a0 = c1*y10 + c2*y1ii + c3*y1iv
Direct port of fastDMA4_cmat().
"""
k = float(k)
k2 = k * k
k3 = k2 * k
k4 = k3 * k
den = 180 + 72 * k - 800 * k2 - 320 * k3 + 320 * k4 + 128 * k * k4
c1 = (180 - 750 * k - 525 * k2 + 450 * k3 + 225 * k4) / den
c2 = (1575 - 1050 * k - 1050 * k2) / den
c3 = 945.0 / den
return {"c1": c1, "c2": c2, "c3": c3}
# =========================================================
# Order 2: all scales (cache-based)
# =========================================================
def fastDMA2_F2_all(y: np.ndarray, scales: Union[np.ndarray, List[int]], i_refresh: int) -> np.ndarray:
"""
Compute F2 for all scales for order=2, using cache-based refresh and
reuse of previous-scale cached moments. Direct port of fastDMA2_F2_all().
"""
y = np.asarray(y, dtype=float)
n = y.size
scales = np.asarray(scales, dtype=int)
nr = scales.size
i_refresh = int(i_refresh)
if i_refresh < 1:
i_refresh = 1
if i_refresh > n:
i_refresh = n
cache_len = n // i_refresh + 1
loc_sum0 = np.zeros(cache_len, dtype=float) # Σ y
loc_sum1 = np.zeros(cache_len, dtype=float) # Σ y*d
loc_sum2 = np.zeros(cache_len, dtype=float) # Σ y*d^2
F2 = np.zeros(nr, dtype=float)
prev_scale = 0
for k2_idx in range(nr):
scale = int(scales[k2_idx])
k = (scale - 1) // 2
ic = k # 0-based
itemp = n - scale
cc = fastDMA2_cmat(k)
c1 = cc["c1"]
c2 = cc["c2"]
iloc = 0 # 0-based cache slot index
if prev_scale == 0:
i0 = 1 # R's 1L
y10 = 0.0
y1i = 0.0
y1ii = 0.0
nn1 = 0.0
nn2 = 0.0
else:
kb = (prev_scale - 1) // 2
nn1 = float(kb - k)
nn2 = nn1 * nn1
i0 = prev_scale + 1 # R's index in window coordinates
y10 = loc_sum0[iloc]
y1i = loc_sum1[iloc] + y10 * nn1
y1ii = loc_sum2[iloc] + 2.0 * loc_sum1[iloc] * nn1 + y10 * nn2
# add missing contributions for indices i0..scale (R's 1..scale window coordinates)
if i0 <= scale:
idx = np.arange(i0, scale + 1, dtype=float) # 1..scale
di = idx - (k + 1) # (idx - k - 1L)
yi = y[(idx.astype(int) - 1)] # map 1-based positions to 0-based y indices
y10 += float(np.sum(yi))
y1i += float(np.sum(yi * di))
y1ii += float(np.sum(yi * di * di))
# cache initial window moments
loc_sum0[iloc] = y10
loc_sum1[iloc] = y1i
loc_sum2[iloc] = y1ii
a0 = c1 * y10 + c2 * y1ii
d = y[ic] - a0
f2 = d * d
if itemp > 0:
k1 = float(k + 1)
kk = float(k)
i = 1
while i <= itemp:
next_ref = ((i // i_refresh) + 1) * i_refresh
if next_ref > itemp:
next_ref = itemp
if i <= next_ref - 1:
for j in range(i, next_ref):
# j is shift index (R: j), window is y[j : j+scale-1] in 0-based
y00 = y10
y0i = y1i
y0ii = y1ii
# rolling Σy
y10 = y00 + y[j + scale - 1] - y[j - 1]
temp1 = y[j - 1] * k1
temp2 = y[j + scale - 1] * kk
y1i = y0i - y00 + temp1 + temp2
y1ii = y0ii - 2.0 * y0i + y00 - temp1 * k1 + temp2 * kk
a0 = c1 * y10 + c2 * y1ii
d = y[j + ic] - a0
f2 += d * d
if next_ref <= itemp:
iloc += 1
if prev_scale == 0:
i0 = 1
y10 = 0.0
y1i = 0.0
y1ii = 0.0
else:
i0 = prev_scale + 1
y10 = loc_sum0[iloc]
y1i = loc_sum1[iloc] + y10 * nn1
y1ii = loc_sum2[iloc] + 2.0 * loc_sum1[iloc] * nn1 + y10 * nn2
if i0 <= scale:
jj = np.arange(i0, scale + 1, dtype=float)
idx_y = next_ref + jj # R: next_ref + jj (1-based in y)
di = jj - (k + 1)
yi = y[(idx_y.astype(int) - 1)]
y10 += float(np.sum(yi))
y1i += float(np.sum(yi * di))
y1ii += float(np.sum(yi * di * di))
loc_sum0[iloc] = y10
loc_sum1[iloc] = y1i
loc_sum2[iloc] = y1ii
a0 = c1 * y10 + c2 * y1ii
d = y[next_ref + ic] - a0
f2 += d * d
i = next_ref + 1
F2[k2_idx] = f2 / float(itemp + 1)
prev_scale = scale
return F2
# =========================================================
# Order 4: all scales (cache-based)
# =========================================================
def fastDMA4_F2_all(y: np.ndarray, scales: Union[np.ndarray, List[int]], i_refresh: int) -> np.ndarray:
"""
Compute F2 for all scales for order=4, cache-based.
Direct port of fastDMA4_F2_all().
"""
y = np.asarray(y, dtype=float)
n = y.size
scales = np.asarray(scales, dtype=int)
nr = scales.size
i_refresh = int(i_refresh)
if i_refresh < 1:
i_refresh = 1
if i_refresh > n:
i_refresh = n
cache_len = n // i_refresh + 1
loc_sum0 = np.zeros(cache_len, dtype=float)
loc_sum1 = np.zeros(cache_len, dtype=float)
loc_sum2 = np.zeros(cache_len, dtype=float)
loc_sum3 = np.zeros(cache_len, dtype=float)
loc_sum4 = np.zeros(cache_len, dtype=float)
F2 = np.zeros(nr, dtype=float)
prev_scale = 0
for k2_idx in range(nr):
scale = int(scales[k2_idx])
k = (scale - 1) // 2
ic = k
itemp = n - scale
cc = fastDMA4_cmat(k)
c1, c2, c3 = cc["c1"], cc["c2"], cc["c3"]
iloc = 0
if prev_scale == 0:
i0 = 1
y10 = y1i = y1ii = y1iii = y1iv = 0.0
nn1 = nn2 = nn3 = nn4 = 0.0
else:
kb = (prev_scale - 1) // 2
nn1 = float(kb - k)
nn2 = nn1 * nn1
nn3 = nn2 * nn1
nn4 = nn2 * nn2
i0 = prev_scale + 1
y10 = loc_sum0[iloc]
y1i = loc_sum1[iloc] + y10 * nn1
y1ii = loc_sum2[iloc] + 2.0 * loc_sum1[iloc] * nn1 + y10 * nn2
y1iii = loc_sum3[iloc] + 3.0 * loc_sum2[iloc] * nn1 + 3.0 * loc_sum1[iloc] * nn2 + y10 * nn3
y1iv = (
loc_sum4[iloc]
+ 4.0 * loc_sum3[iloc] * nn1
+ 6.0 * loc_sum2[iloc] * nn2
+ 4.0 * loc_sum1[iloc] * nn3
+ y10 * nn4
)
if i0 <= scale:
idx = np.arange(i0, scale + 1, dtype=float)
di = idx - (k + 1)
yi = y[(idx.astype(int) - 1)]
di2 = di * di
y10 += float(np.sum(yi))
y1i += float(np.sum(yi * di))
y1ii += float(np.sum(yi * di2))
y1iii += float(np.sum(yi * di2 * di))
y1iv += float(np.sum(yi * di2 * di2))
loc_sum0[iloc] = y10
loc_sum1[iloc] = y1i
loc_sum2[iloc] = y1ii
loc_sum3[iloc] = y1iii
loc_sum4[iloc] = y1iv
a0 = c1 * y10 + c2 * y1ii + c3 * y1iv
d = y[ic] - a0
f2 = d * d
if itemp > 0:
k1 = float(k + 1)
kk = float(k)
i = 1
while i <= itemp:
next_ref = ((i // i_refresh) + 1) * i_refresh
if next_ref > itemp:
next_ref = itemp
if i <= next_ref - 1:
for j in range(i, next_ref):
y00 = y10
y0i = y1i
y0ii = y1ii
y0iii = y1iii
y0iv = y1iv
y10 = y00 + y[j + scale - 1] - y[j - 1]
# update moments (order-4 recurrence)
t1 = y[j - 1] * k1
t2 = y[j + scale - 1] * kk
y1i = y0i - y00 + t1 + t2
t1 = t1 * k1
t2 = t2 * kk
y1ii = y0ii - 2.0 * y0i + y00 - t1 + t2
t1 = t1 * k1
t2 = t2 * kk
y1iii = y0iii - 3.0 * y0ii + 3.0 * y0i - y00 + t1 + t2
y1iv = y0iv - 4.0 * y0iii + 6.0 * y0ii - 4.0 * y0i + y00 - t1 * k1 + t2 * kk
a0 = c1 * y10 + c2 * y1ii + c3 * y1iv
d = y[j + ic] - a0
f2 += d * d
if next_ref <= itemp:
iloc += 1
if prev_scale == 0:
i0 = 1
y10 = y1i = y1ii = y1iii = y1iv = 0.0
else:
i0 = prev_scale + 1
y10 = loc_sum0[iloc]
y1i = loc_sum1[iloc] + y10 * nn1
y1ii = loc_sum2[iloc] + 2.0 * loc_sum1[iloc] * nn1 + y10 * nn2
y1iii = (
loc_sum3[iloc]
+ 3.0 * loc_sum2[iloc] * nn1
+ 3.0 * loc_sum1[iloc] * nn2
+ y10 * nn3
)
y1iv = (
loc_sum4[iloc]
+ 4.0 * loc_sum3[iloc] * nn1
+ 6.0 * loc_sum2[iloc] * nn2
+ 4.0 * loc_sum1[iloc] * nn3
+ y10 * nn4
)
if i0 <= scale:
jj = np.arange(i0, scale + 1, dtype=float)
idx_y = next_ref + jj
di = jj - (k + 1)
yi = y[(idx_y.astype(int) - 1)]
di2 = di * di
y10 += float(np.sum(yi))
y1i += float(np.sum(yi * di))
y1ii += float(np.sum(yi * di2))
y1iii += float(np.sum(yi * di2 * di))
y1iv += float(np.sum(yi * di2 * di2))
loc_sum0[iloc] = y10
loc_sum1[iloc] = y1i
loc_sum2[iloc] = y1ii
loc_sum3[iloc] = y1iii
loc_sum4[iloc] = y1iv
a0 = c1 * y10 + c2 * y1ii + c3 * y1iv
d = y[next_ref + ic] - a0
f2 += d * d
i = next_ref + 1
F2[k2_idx] = f2 / float(itemp + 1)
prev_scale = scale
return F2
# =========================================================
# User-facing fastDMA()
# =========================================================
def fastDMA(
x: Union[np.ndarray, List[float]],
order: int = 0,
minbox: Optional[int] = None,
maxbox: int = 0,
boxratio: float = 2 ** (3 / 16),
integrate: bool = True,
i_refresh: Optional[int] = None,
) -> Dict[str, np.ndarray]:
"""
Python equivalent of R fastDMA().
Returns dict: {order, scales, F2, log10_scale, log10_F}
"""
x = np.asarray(x, dtype=float)
if x.size < 10:
raise ValueError("x must have length >= 10.")
order = int(order)
if order not in (0, 2, 4):
raise ValueError("order must be one of 0, 2, 4.")
if minbox is None:
minbox = 5 if order == 0 else (7 if order == 2 else 9)
if i_refresh is None:
i_refresh = 1000 if order == 0 else (200 if order == 2 else 20)
minbox = int(minbox)
maxbox = int(maxbox)
i_refresh = int(i_refresh)
# integrate step
if integrate:
x0 = x - np.nanmean(x)
y = np.cumsum(x0)
else:
y = x.copy()
# window formatting: odd, min required
min_required = 5 if order == 0 else (7 if order == 2 else 9)
minbox = (minbox // 2) * 2 + 1
if minbox < min_required:
minbox = min_required
n = y.size
if maxbox <= 0 or maxbox > (n // 4) or minbox > maxbox:
maxbox = n // 4
scales = fastDMA_rscale(minbox, maxbox, boxratio)
if order == 0:
F2 = np.array([fastDMA0_F2_one(y, int(s), i_refresh) for s in scales], dtype=float)
elif order == 2:
F2 = fastDMA2_F2_all(y, scales, i_refresh)
else:
F2 = fastDMA4_F2_all(y, scales, i_refresh)
log10_scale = np.log10(scales.astype(float))
log10_F = 0.5 * np.log10(F2)
return {
"order": np.array([order], dtype=int),
"scales": scales,
"F2": F2,
"log10_scale": log10_scale,
"log10_F": log10_F,
}
# =========================================================
# Scaling exponent estimation (with order auto-correction)
# =========================================================
def estimate_scaling_fastDMA(
x: Union[np.ndarray, List[float]],
s1: float,
s2: float,
order: int = 0,
minbox: Optional[int] = None,
maxbox: int = 0,
boxratio: float = 2 ** (3 / 16),
integrate: bool = True,
i_refresh: Optional[int] = None,
scale_correction: bool = True,
main: Optional[str] = None,
pch: str = "o",
markersize: float = 4.0,
) -> Dict[str, object]:
"""
Python equivalent of estimate_scaling_fastDMA().
- odd order auto-corrected to nearest supported even order (1->0, 3->2, 5->4)
- scale correction factors: 0/1 -> /1, 2/3 -> /1.93, 4/5 -> /2.74
- fits log10(F) ~ log10(s) on points where raw scale s_raw in [s1,s2]
- plots log-log points and fitted line, with guide segments at s1,s2
Returns dict including slope, intercept, r2, indices used, elapsed time, etc.
"""
t_start = time.perf_counter()
x = np.asarray(x, dtype=float)
n = x.size
if n < 10:
raise ValueError("x must have length >= 10.")
order_input = int(order)
if order_input == 1:
order_used = 0
elif order_input == 3:
order_used = 2
elif order_input == 5:
order_used = 4
else:
order_used = order_input
if order_used not in (0, 2, 4):
raise ValueError("order must be one of 0,1,2,3,4,5 (odd values are auto-corrected).")
if not (np.isfinite(s1) and np.isfinite(s2) and s1 > 0 and s2 > 0):
raise ValueError("s1 and s2 must be positive finite numbers.")
if s1 > s2:
s1, s2 = s2, s1
out = fastDMA(
x=x,
order=order_used,
minbox=minbox,
maxbox=int(maxbox),
boxratio=boxratio,
integrate=integrate,
i_refresh=i_refresh,
)
s_raw = out["scales"].astype(float)
F = np.sqrt(out["F2"].astype(float))
corr_factor = 1.0
if scale_correction:
if order_input in (0, 1):
corr_factor = 1.0
elif order_input in (2, 3):
corr_factor = 1.93
elif order_input in (4, 5):
corr_factor = 2.74
else:
corr_factor = 1.0
s = s_raw / corr_factor
log10_s = np.log10(s)
log10_F = np.log10(F)
ok = np.isfinite(log10_s) & np.isfinite(log10_F) & (s_raw >= s1) & (s_raw <= s2)
if np.sum(ok) < 2:
raise ValueError("Not enough points within [s1, s2] to fit a line.")
# linear regression in log10 space
x_fit = log10_s[ok]
y_fit = log10_F[ok]
slope, intercept = np.polyfit(x_fit, y_fit, deg=1)
y_hat = intercept + slope * x_fit
ss_res = float(np.sum((y_fit - y_hat) ** 2))
ss_tot = float(np.sum((y_fit - float(np.mean(y_fit))) ** 2))
r2 = 1.0 - ss_res / ss_tot if ss_tot > 0 else float("nan")
elapsed_sec = time.perf_counter() - t_start
# plot
if main is None:
if scale_correction:
main = f"fastDMA (order={order_used}): fit in [{s1:g}, {s2:g}] (s/{corr_factor:.3g})"
else:
main = f"fastDMA (order={order_used}): fit in [{s1:g}, {s2:g}]"
plt.figure(figsize=(7, 5))
plt.plot(log10_s, log10_F, pch, markersize=markersize, label="all scales")
plt.plot(log10_s[ok], log10_F[ok], pch, markersize=markersize + 1.0, label="used for fit")
# guide segments at s1,s2 (converted to corrected scale for x position)
x_v = np.log10(np.array([s1, s2], dtype=float) / corr_factor)
# y_center via interpolation on all points
y_center = np.interp(x_v, log10_s, log10_F)
y_rng = np.nanmax(log10_F) - np.nanmin(log10_F)
half_h = y_rng / 10.0 if np.isfinite(y_rng) and y_rng > 0 else 0.1
for xv, yc in zip(x_v, y_center):
plt.plot([xv, xv], [yc - half_h, yc + half_h], linestyle="--", linewidth=1.5)
# fitted line over fit range
xx = np.linspace(np.min(x_fit), np.max(x_fit), 200)
yy = intercept + slope * xx
plt.plot(xx, yy, linestyle="--", linewidth=2, label="fit")
plt.title(main)
plt.xlabel(r"$\log_{10}(s)$")
plt.ylabel(r"$\log_{10}(F(s))$")
plt.legend(loc="best")
# annotation
txt = (
f"order = {order_used} (input={order_input})\n"
+ (f"scale: s/{corr_factor:.3g} (s1,s2 raw)\n" if scale_correction else "")
+ f"slope = {slope:.5f}\n"
+ f"R^2 = {r2:.4f}\n"
+ f"n = {int(np.sum(ok))}\n"
+ f"elapsed = {elapsed_sec:.3f} s"
)
ax = plt.gca()
xmin, xmax = ax.get_xlim()
ymin, ymax = ax.get_ylim()
plt.text(xmin + 0.02 * (xmax - xmin), ymax - 0.02 * (ymax - ymin), txt,
ha="left", va="top")
plt.tight_layout()
plt.show()
return {
"out": out,
"slope": float(slope),
"intercept": float(intercept),
"r2": float(r2),
"n_used": int(np.sum(ok)),
"idx_used": np.where(ok)[0],
"s_range_input": (float(s1), float(s2)),
"order_input": int(order_input),
"order_used": int(order_used),
"scale_correction": bool(scale_correction),
"scale_divisor": float(corr_factor),
"elapsed_sec": float(elapsed_sec),
}
# =========================================================
# Example (equivalent to the R example)
# =========================================================
if __name__ == "__main__":
# Example: x is long-range-correlated data (here: just white noise, same as your R example)
N = 10000
rng = np.random.default_rng(0)
x = rng.normal(size=N)
# DMA2 example: estimate scaling exponent and plot
res_fit = estimate_scaling_fastDMA(
x, s1=15, s2=len(x) / 20,
order=2,
scale_correction=True
)
print("Estimated scaling exponent (slope):", res_fit["slope"])
実行例は以下です。
# Example: x is long-range-correlated data (here: just white noise, same as your R example)
N = 100000
rng = np.random.default_rng(0)
x = rng.normal(size=N)
# DMA2 example: estimate scaling exponent and plot
res_fit = estimate_scaling_fastDMA(
x, s1=15, s2=len(x) / 20,
order=2,
scale_correction=True
)
print("Estimated scaling exponent (slope):", res_fit["slope"])
※ もし記事の中で「ここ違うよ」という点や気になるところがあれば、気軽に指摘していただけると助かります。質問や「このテーマも取り上げてほしい」といったリクエストも大歓迎です。必ず対応するとは約束できませんが、できるだけ今後の記事で扱いたいと思います。それと、下のはてなブログランキングはあまり信用できる指標ではなさそうですが (私のブログを読んでいる人は、実際とても少ないです)、押してもらえるとシンプルに励みになります。気が向いたときにポチッとしていただけたら嬉しいです。
