ケィオスの時系列解析メモランダム

時系列解析、生体情報解析などをやわらかく語ります

【フラクタル時系列解析の基礎】理想と現実のあいだで考えるフラクタルゆらぎ：ブラウン運動と非整数ブラウン運動

非整数ブラウン運動 確率過程

ブラウン運動と非整数ブラウン運動（fractional Brownian motion）について、しっかり理解したい人向けに、基礎の基礎から説明してみたいと思います。

1. 生物学者ブラウンの発見から、数学的ブラウン運動へ

　ブラウン運動という言葉は、もともと生物学者ロバート・ブラウン（Robert Brown, 1773–1858）が、水中の微粒子が不規則に動く様子を観察したことに由来します。この「なんだかランダムに動いている」現象を、数学的に洗練させたものが、数学的に定義されるブラウン運動です。

　数学の世界で定義される1次元ブラウン運動には、次のような特徴があります。

  • 時間は連続
  • グラフはどこを拡大してもギザギザが続く
  • しかし、線は途切れず連続した曲線になっている

さらに重要なのが、拡大しても似た揺らぎが現れるという性質です。 具体的には、ブラウン運動のグラフの一部について

  • 横幅（時間幅）を a 倍
  • 縦幅（変動幅）を a^ {0.5} 倍

（たとえば、横を２倍、縦を 2^ {0.5} = \sqrt{2} 倍）すると、見た目としては元のブラウン運動とよく似た変動が現れます。ここでいう「よく似ている」とは、完全に同じ形になるという意味ではありません。変動の大きさや粗さといった統計的な特徴が同じ、という意味です。この特徴をグラフで描いたのが下のアニメです。

数学的ブラウン運動のグラフは、ある部分の横幅（時間幅）を a 倍し、縦幅を a^ {0.5} 倍すると、元のグラフと似た変動が見られる。数学的世界では、何度でも（無限回）拡大できる。

　上図では拡大を途中でやめていますが、この拡大操作は、数学で想定している世界では何度でも繰り返せます。なぜなら、数学的ブラウン運動は、無限小の極限まで考えられる理想的な存在だからです。

　このブラウン運動の性質をさらに一般化したものが、非整数ブラウン運動（fractional Brownian motion）です。非整数ブラウン運動のグラフは、ある部分の横幅を a 倍したときに、縦幅を

\displaystyle{ a^H }

倍すると似た変動が現れるようなものを考えます。 H = 0.5 のときが、ブラウン運動ですので、 H を自由に変えてみようということです。 ここで H は ハースト（Hurst）指数 と呼ばれ、

\displaystyle{ 0 < H < 1 }

の範囲にあります。 H の値によって以下のような違いが見られます。

  • \displaystyle{H = 0.5}：通常のブラウン運動です。
  • \displaystyle{H \gt 0.5}：ブラウン運動よりも上下に広がりやすい。ブラウン運動より「なめらか」に見えますが、ギザギザはあるので微分はできません。
  • \displaystyle{H \lt 0.5}：ブラウン運動よりも上下の広がりにくい。ブラウン運動よりも、ギザギザが目立つように見えます。

\displaystyle{H = 0.8} の非整数ブラウン運動

\displaystyle{H = 0.2} の非整数ブラウン運動

　数学的なブラウン運動や非整数ブラウン運動は、ルベーグ積分や確率解析、非整数微積分といった数学を用いて厳密に議論されるため、「なんだか難しそう」と感じる人も多いと思います。数学を専門としない人にとって、こうした数学の基礎を一から学ぶのはなかなか大変です（正直に言えば、私にも簡単ではありません）。ですので、ここでは難しい話はいったん聞かなかったことにして、 H の値が変わると、グラフの広がり具合やギザギザの度合いが変わる――その点だけを印象として持ってもらえれば十分です。

2. 実世界で観測した時系列は「離散的」

　フラクタル時系列解析では、「現実のデータにも、ブラウン運動や非整数ブラウン運動と似た特徴が現れる」と考えます。ただし、ここで最初に強調しておきたい点があります。

理想と現実は同じではありません。

もし実世界の時系列が、数学的に理想化されたブラウン運動そのものであれば、どれだけ拡大しても、永遠に似たゆらぎが見え続けるはずです。しかし、実際のデータではそうなりません。

  • 拡大を何度か繰り返すと
  • あるところで点がバラバラに見え始め
  • それ以上拡大すると、元と似た構造は消えてしまいます

要するに、拡大しすぎるとスカスカ（点）になるのです。このこと示したのが下図で、上に示した理想のブラウン運動・非整数ブラウン運動との違いを理解してください。

ブラウン運動のような離散時系列。何度か拡大すると、少数の点が現れる。

　どこまでも拡大し続けられないのは当然で、実世界のデータは、離散時間でサンプリングされた離散時系列、つまり、有限個の離散的な点の集まりだからです。

　私たちが時系列解析で実際にやっているのは、無限小の極限を見ることではなく、有限のスケールの範囲で眺めて「この範囲ではブラウン運動っぽい」、「この範囲では非整数ブラウン運動っぽい」と判断しているにすぎません。

　ここを理解しておかないと、数学的な議論と実データ解析を混同してしまいます。

3. 理想は連続軌道、現実は離散データから始まる

　時系列解析や信号処理では、ランダムな信号の代表としてホワイトノイズが登場します。ホワイトノイズとは、時間的に相関のないランダムな数値列です。

　ここで少し注意が必要なのは、連続ブラウン運動の「差分」として定義されるホワイトノイズは、通常の意味でグラフを描くことができないという点です。

　数学的に理想化された立場では、まず連続時間のブラウン運動が定義されます。ブラウン運動は、どこを拡大してもギザギザが続くものの、曲線としては連続です。 ところが、そのブラウン運動を「時間で微分したもの」あるいは「無限小の時間差で引き算したもの」としてホワイトノイズを考えると、状況は一変します。

　無限小の時間間隔でブラウン運動の値の変化を見ようとすると、その変化量はどんどん不安定になり、有限な値として落ち着きません。その結果、ホワイトノイズは

  • 各時刻で値が定まっている「関数」ではなく
  • 連続した曲線として描けるものでもなく
  • 時間軸に沿って線で結ぶことができない存在

になります。言い換えると、ホワイトノイズは「この時刻ではこの値」と点を並べて描くようなものではなく、無限小の時間解像度で眺めると、そもそも“形”を持たないのです。

このため、数学的な意味でのホワイトノイズは、

  • グラフを描く対象ではなく
  • 積分された結果（＝ブラウン運動）を通してのみ意味を持つ

という、かなり抽象的な存在になります。

　だからこそ、数学の教科書では、

  1. まず連続で扱いやすいブラウン運動を定義し
  2. その「差分」あるいは「形式的な微分」としてホワイトノイズを導入する

という順序が取られるのです。

一方で、現実のデータ解析では、そもそも時間は離散的であり、観測点も有限です。そのため、私たちはホワイトノイズを「とびとびの時刻で、互いに相関のないランダムな数値列」として自然に扱うことができます。この点が、数学的理想と現実の解析との大きな違いです。

　そして、離散的ホワイトノイズを累積（積分）したものを「ブラウン運動のようなもの」と理解している人も多いと思います。ただし、厳密に言えば、それは、数学的に理想化されたブラウン運動そのものではありません。

　そのため、少し慎重な立場の人は、これを「ランダムウォーク」と呼びます。物理学では、離散時間のランダムウォークを、時間間隔を無限小にして連続化した極限が数学的ブラウン運動という立場をとることがあります。本来は、時間間隔を無限小にすると、分子の運動の世界になりますが、そこは無視して理想化する立場です。

4. おわりに：用語の「いい加減さ」と上手につきあう

　時系列解析の知識がある程度ついてくると、 離散時間の時系列なのに「ブラウン運動」とか「非整数ブラウン運動」とか、いい加減に用語を使ってしまうことがあります。しかし、最近になって、この曖昧さこそが入門者を混乱させる原因なのではないかと感じるようになりました。

　ですので、まず、次の点を意識して、違いを意識して考えてください。

  • 数学的理想の世界に住む、 ブラウン運動・非整数ブラウン運動
  • 実世界のデータで、それっぽく見えるブラウン運動・非整数ブラウン運動

これらは同一ではありません。

　実世界の話は実は単純で、「厳密には違うけれど、途中のスケール（ある範囲の大きさ）ではよく似ている」――ただそれだけのことです。

　この視点を持っておくと、フラクタル時系列解析やスケーリング解析が、特定のスケールの領域に注目する理由が分かるはずです。

[付録] Rで非整数ブラウン運動のサンプルパスを生成

　離散時系列としての非整数ブラウン運動（fractional Brownian motion; fBm）のサンプルを、実際に解析してみたい人のために、ここでは サンプルパス（時系列の実現値）を生成する R スクリプトを掲載します。 以下のコードは、非整数ブラウン運動の数学的定義に整合するよう設計されており、数値シミュレーションの標準的手法である Davies–Harte 法を用いています。

############################################################
# Fractional Brownian motion (fBm) sample path generator
#
# Definition matched:
#   - Gaussian process
#   - stationary increments
#   - self-similar with exponent H
#   - Var(B_H(t+dt)-B_H(t)) = dt^(2H)
#
# Output:
#   x: length N, with x[1] = 0 and time step dt
############################################################

# ------------------------------------------------------------
# Autocovariance of fractional Gaussian noise (fGn), k = 0..(N-1)
# gamma(k) = 0.5*(|k+1|^(2H) - 2|k|^(2H) + |k-1|^(2H))
# ------------------------------------------------------------
fgn_acvf <- function(N, H) {
  k <- 0:(N - 1)
  # Use abs() to be safe, though k is nonnegative here
  g <- 0.5 * (abs(k + 1)^(2 * H) - 2 * abs(k)^(2 * H) + abs(k - 1)^(2 * H))
  g
}

# ------------------------------------------------------------
# Generate fGn of length N using Davies–Harte
# Returns increments with Var(increment) = 1 for dt = 1
# ------------------------------------------------------------
fgn_daviesharte <- function(N, H, seed = NULL) {
  if (!is.null(seed)) set.seed(seed)
  if (!(H > 0 && H < 1)) stop("H must satisfy 0 < H < 1.")
  if (N < 2) stop("N must be >= 2.")

  n2 <- 2 * N

  # Circulant embedding vector (length 2N)
  g <- fgn_acvf(N, H)                 # gamma(0..N-1)
  cvec <- c(g, 0, rev(g[2:N]))        # length = 2N

  # Eigenvalues via FFT of the circulant first row
  lam <- Re(fft(cvec))

  # Numerical guard: tiny negatives can appear from floating error
  tol <- 1e-12 * max(lam)
  lam[lam < 0 & lam > -tol] <- 0
  if (any(lam < 0)) {
    stop("Davies–Harte failed: negative eigenvalues (try different N or check numeric stability).")
  }

  # Build complex Gaussian vector Z with Hermitian symmetry in distribution
  # Z[1] and Z[N+1] are real; others are complex with N(0,1/2)+iN(0,1/2)
  Z <- complex(length.out = n2)
  Z[1]       <- rnorm(1)
  Z[N + 1]   <- rnorm(1)
  if (N > 1) {
    Z[2:N]         <- (rnorm(N - 1) + 1i * rnorm(N - 1)) / sqrt(2)
    Z[(N + 2):n2]  <- Conj(rev(Z[2:N]))
  }

  # Scale and inverse FFT
  W <- sqrt(lam / n2) * Z
  y <- fft(W, inverse = TRUE) / n2

  # fGn increments: take first N real parts
  inc <- Re(y[1:N])

  inc
}

# ------------------------------------------------------------
# fBm generator: cumulative sum of fGn increments
# x[1] = 0, x has length N, time step dt
# ------------------------------------------------------------
fbm_daviesharte <- function(N, H, dt = 1, seed = NULL) {
  if (dt <= 0) stop("dt must be positive.")
  inc <- fgn_daviesharte(N - 1, H, seed = seed)  # N-1 increments to make N points
  inc <- inc * (dt^H)                            # scale increments to time step dt
  x <- c(0, cumsum(inc))
  x
}

　非整数ブラウン運動そのものは連続時間・連続値の確率過程として定義されますが、実データ解析や数値実験では、有限長・離散時間の時系列として扱う必要があります。本スクリプトは、そのような解析用の離散サンプルを生成することを目的としています。

このスクリプトで何をしているか

　コードは大きく分けて、次の3つのステップから構成されています。

  1. 非整数ガウス雑音（fractional Gaussian noise; fGn）の自己共分散関数を定義する fBm の増分は fGn と呼ばれ、時間的に相関をもつガウス過程です。 この相関構造は、ハースト指数 H によって決まり、
\displaystyle{ \gamma(k) = \frac{1}{2}\Bigl(|k+1|^{2H} - 2|k|^{2H} + |k-1|^{2H}\Bigr) }

という明示的な形で与えられます。

  1. Davies–Harte 法により、指定した共分散構造をもつ fGn を高速に生成する Davies–Harte 法は、循環行列埋め込みと FFT を用いて、 理論通りの共分散をもつガウス時系列を O(N \log N) で生成できる方法です。 数値誤差により固有値がわずかに負になる場合があるため、コード中では安全のためのチェックも行っています。

  2. 生成した fGn を累積和（離散積分）して fBm を得る 離散時系列としての fBm は、fGn を時間方向に累積したものとして構成できます。 このとき、時間刻み (dt) に対して増分の分散が dt^ {2H} となるよう、 増分を dt^ H 倍しています。

注意点

  • ここで生成されるのは、数学的に理想化された連続時間 fBm そのものではなく、 それを有限長・離散時間で近似したサンプルパスです。
  • 拡大を無限に繰り返しても同じ構造が現れるわけではなく、 あくまで 有限のスケール範囲で fBm 的な性質を示す時系列であることに注意してください。
  • この点を理解した上で用いると、DFA、DMA、スケーリング解析などの 数値実験用データとして活用できます。

実行例

　実際に動作を確認するには、まず上の関数定義部分を実行し、その後、以下の簡単なスクリプトを実行してください。 指定したハースト指数 H に対応する、非整数ブラウン運動らしいギザギザした軌道が描画されるはずです。

N  <- 4096
H  <- 0.8
dt <- 1.0
x  <- fbm_daviesharte(N, H, dt = dt)

plot((0:(N-1))*dt, x, type = "l", xlab = "t", ylab = "B_H(t)",
      main = sprintf("fBm sample path (H=%.2f, N=%d)", H, N))

[付録] Python で非整数ブラウン運動のサンプルパスを生成

　Python 版では、R 版と同様に Davies–Harte 法を用いて、非整数ブラウン運動（fractional Brownian motion; fBm）のサンプルパスを生成します。 数値的な考え方や処理の流れは R 版と完全に同一であり、言語だけを Python に置き換えた実装になっています。

　まず、fgn_acvf 関数では、非整数ブラウン運動の増分に対応する 非整数ガウス雑音（fractional Gaussian noise; fGn）の自己共分散関数を定義しています。 この共分散構造は、ハースト指数 H によって決まり、増分同士が長距離相関をもつことを明示的に反映しています。次に、fgn_daviesharte 関数では、この共分散構造をもつ fGn を 高速に生成します。

　fbm_daviesharte 関数では、生成した fGn を時間方向に 累積和（離散積分）することで、非整数ブラウン運動のサンプルパスを構成しています。 時間刻み dt に対して、増分の分散が

\displaystyle{ \mathrm{Var}[B_H(t+dt)-B_H(t)] = dt^{2H} }

となるよう、増分を dt^ H 倍している点が重要です。

############################################################
# Fractional Brownian motion (fBm) sample path generator
# (Davies–Harte method)
#
# Definition matched:
#   - Gaussian process
#   - stationary increments
#   - self-similar with exponent H
#   - Var(B_H(t+dt)-B_H(t)) = dt^(2H)
#
# Requirements:
#   numpy only
############################################################

import numpy as np


# ------------------------------------------------------------
# Autocovariance of fractional Gaussian noise (fGn)
# gamma(k) = 0.5*(|k+1|^(2H) - 2|k|^(2H) + |k-1|^(2H))
# k = 0..N-1
# ------------------------------------------------------------
def fgn_acvf(N, H):
    k = np.arange(N)
    g = 0.5 * (
        np.abs(k + 1)**(2 * H)
        - 2 * np.abs(k)**(2 * H)
        + np.abs(k - 1)**(2 * H)
    )
    return g


# ------------------------------------------------------------
# Generate fGn of length N using Davies–Harte
# Returns increments with Var(increment) = 1 for dt = 1
# ------------------------------------------------------------
def fgn_daviesharte(N, H, seed=None):
    if seed is not None:
        np.random.seed(seed)

    if not (0 < H < 1):
        raise ValueError("H must satisfy 0 < H < 1.")
    if N < 2:
        raise ValueError("N must be >= 2.")

    n2 = 2 * N

    # Circulant embedding vector (length 2N)
    g = fgn_acvf(N, H)                  # gamma(0..N-1)
    cvec = np.concatenate([g, [0.0], g[1:][::-1]])

    # Eigenvalues via FFT
    lam = np.real(np.fft.fft(cvec))

    # Numerical guard against tiny negative values
    tol = 1e-12 * np.max(lam)
    lam[(lam < 0) & (lam > -tol)] = 0.0
    if np.any(lam < 0):
        raise RuntimeError(
            "Davies–Harte failed: negative eigenvalues "
            "(try different N or check numeric stability)."
        )

    # Build complex Gaussian vector Z
    Z = np.zeros(n2, dtype=np.complex128)

    Z[0]     = np.random.randn()
    Z[N]     = np.random.randn()

    if N > 1:
        re = np.random.randn(N - 1)
        im = np.random.randn(N - 1)
        Z[1:N] = (re + 1j * im) / np.sqrt(2)
        Z[N+1:] = np.conj(Z[1:N][::-1])

    # Scale and inverse FFT
    W = np.sqrt(lam / n2) * Z
    y = np.fft.ifft(W)

    # fGn increments
    inc = np.real(y[:N])

    return inc


# ------------------------------------------------------------
# fBm generator: cumulative sum of fGn increments
# x[0] = 0, length N, time step dt
# ------------------------------------------------------------
def fbm_daviesharte(N, H, dt=1.0, seed=None):
    if dt <= 0:
        raise ValueError("dt must be positive.")

    inc = fgn_daviesharte(N - 1, H, seed=seed)
    inc = inc * (dt ** H)

    x = np.empty(N)
    x[0] = 0.0
    x[1:] = np.cumsum(inc)

    return x

実行例

　続く実行例では、サンプル長 N、ハースト指数 H、時間刻み dt を指定し、 生成した非整数ブラウン運動のサンプルパスを Matplotlib を用いて描画しています。

import numpy as np
import matplotlib.pyplot as plt

# Parameters
N  = 4096
H  = 0.8
dt = 1.0

# Generate fBm sample path
x = fbm_daviesharte(N, H, dt=dt)

# Time axis
t = np.arange(N) * dt

# Plot
plt.plot(t, x)
plt.xlabel("t")
plt.ylabel(r"$B_H(t)$")
plt.title(f"fBm sample path (H={H:.2f}, N={N})")
plt.show()

※ もし記事の中で「ここ違うよ」という点や気になるところがあれば、気軽に指摘していただけると助かります。質問や「このテーマも取り上げてほしい」といったリクエストも大歓迎です。必ず対応するとは約束できませんが、できるだけ今後の記事で扱いたいと思います。それと、下のはてなブログランキングはあまり信用できる指標ではなさそうですが (私のブログを読んでいる人は、実際とても少ないです)、押してもらえるとシンプルに励みになります。気が向いたときにポチッとしていただけたら嬉しいです。

ランキング参加中
知識
ランキング参加中
テクノロジー
ランキング参加中
数学・科学・工学