時系列の長時間相関（スケーリング）を調べるとき、トレンド除去を含む detrended moving-average analysis (DMA) は、DFA の次に良く使われている手法です。実は、DMAには、DFAよりも優れた部分が複数あり、上位互換と言える方法です。しかし、知名度（主には著者の有名さ）により「世界はDFA最高」という残念な風潮があります。
そこで、皆さんにもDMA、とくに、私が開発したDMAの高速アルゴリズムを使ってもらえるように、ここでは、DMA（centered DMA）を高速に計算するための C プログラム sgdma を公開し、使い方を説明します。前回のブログで、この高速アルゴリズムのRスクリプトを紹介しましたが、残念ながら、Rでは、超高速とは言えません。Cプログラムをコンパイルして使えば、Rに比べて、１万倍以上速く処理できます。
以下で紹介するプログラムは、時系列を含むデータファイルを入力とし、各スケール に対する変動関数
を計算し、log10 スケールと log10 変動の 2 列を CSV 形式で出力します。このプログラムを試してみたい人は、以下のコードをコピーし、“sgdma.c” のような名前のファイルとして保存してCコンパイラで実行形式ファイルを作成してください。
/*
sgdma.c
Unified Savitzky-Golay-style detrended moving-average scaling analysis
(centered DMA; detrending orders d = 0, 2, 4)
Original algorithms:
DMA0 / DMA4: Ken Kiyono, 26 Sep. 2015
DMA2: Ken Kiyono, 20 Nov. 2016
Modified by: Ken Kiyono, 27 Jul. 2024
Public release update: 17 Jan. 2026
Notes:
- Some memory-handling / implementation style was inspired by "Numerical Recipes in C"
(Press, Teukolsky, Vetterling, Flannery).
- This implementation uses fast rolling-moment updates (O(N) per scale).
- Optional periodic "refresh" recomputation mitigates floating-point drift.
Citation request:
If you use this code (or a modified version) in published research, please cite:
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).
License:
CC BY 4.0 (Creative Commons Attribution 4.0 International).
You may use/copy/redistribute/modify for any purpose, including commercial use,
provided proper credit is given and changes are indicated.
Output (CSV):
Three columns:
log10(s_eff), log10(F(s)), Emax
where:
s_eff = s / C(d) if scale correction enabled (-S 1, default)
s_eff = s if scale correction disabled (-S 0)
and Emax is a diagnostic: maximum absolute difference (in trend estimate units)
between rolling-update and refresh-recomputed trend (when refresh is enabled).
Default correction constants:
C(0)=1.00, C(2)=1.93, C(4)=2.74
Options:
-d ORDER detrending order: 0, 2, or 4 (default: 0)
-S FLAG scale correction: 1=use s/C(d) (default), 0=use s
-s SKIP skip first SKIP lines (default: 0)
-c COLUMN select 1-based column index (default: 1)
-l MINBOX minimum box width (odd enforced; default depends on d)
-u MAXBOX maximum box width (default: N/4)
-a asymmetric ends (supported for d=0,2; ignored for d=4)
-i disable integration (default: integrate demeaned series)
-r REFRESH refresh interval (0 disables refresh; default depends on d)
-h help
Important:
- The legacy "single-scale trend output mode" is intentionally removed.
*/
#include
#include
#include
#include
#include
#include
/* -------------------------
Small utilities
------------------------- */
static const char *g_progname = "sgdma";
static void die(const char *msg)
{
fprintf(stderr, "%s: %s\n", g_progname, msg);
exit(EXIT_FAILURE);
}
static void *xmalloc(size_t nbytes)
{
void *p = malloc(nbytes);
if (!p) die("Out of memory (malloc).");
return p;
}
static void *xrealloc(void *p, size_t nbytes)
{
void *q = realloc(p, nbytes);
if (!q) die("Out of memory (realloc).");
return q;
}
/* Force odd >= 1 */
static long force_odd(long x)
{
if (x < 1) x = 1;
return (x / 2) * 2 + 1;
}
/* Order-dependent defaults */
static long default_minbox_for_order(int d)
{
if (d == 0) return 5;
if (d == 2) return 7;
return 9; /* d==4 */
}
static long default_refresh_for_order(int d)
{
/* chosen to balance drift suppression and overhead */
if (d == 0) return 2000;
if (d == 2) return 200;
return 20; /* d==4 */
}
static double default_scale_correction_for_order(int d)
{
if (d == 0) return 1.0;
if (d == 2) return 1.93;
return 2.74; /* d==4 */
}
/* -------------------------
CLI / usage
------------------------- */
static void print_usage(void)
{
fprintf(stderr,
"usage: %s [OPTIONS ...] inputfile\n"
"OPTIONS:\n"
" -d ORDER detrending order: 0, 2, or 4 (default: 0)\n"
" -S FLAG scale correction: 1=use s/C(d) (default), 0=use s\n"
" -s SKIP skip first SKIP lines (default: 0)\n"
" -c COLUMN select column index (1-based) (default: 1)\n"
" -l MINBOX minimum box width (odd enforced; default depends on d)\n"
" -u MAXBOX maximum box width (default: N/4)\n"
" -a asymmetric ends (supported for d=0,2; ignored for d=4)\n"
" -i disable integration (default: integrate demeaned series)\n"
" -r REFRESH refresh interval (0 disables refresh; default depends on d)\n"
" -h print this help\n"
"\n"
"Default correction constants C(d): C(0)=1.00, C(2)=1.93, C(4)=2.74\n"
"Output (CSV): log10(s_eff), log10(F(s)), Emax\n",
g_progname
);
}
/* -------------------------
Scale generator (exact behavior preserved)
------------------------- */
static long generate_scales_exact(long minbox, long maxbox, double boxratio, long **out_scales)
{
/* minbox is assumed odd and >= required minimum */
minbox = force_odd(minbox);
if (maxbox < minbox) maxbox = minbox;
const double ratio = (double)maxbox / (double)minbox;
if (!(ratio > 0.0)) die("Invalid scale bounds.");
long rslen = (long)(log10(ratio) / log10(boxratio) + 1.0);
if (rslen < 1) rslen = 1;
long *rs = (long *)xmalloc((size_t)rslen * sizeof(long));
long n = 1;
rs[0] = minbox;
for (long ir = 1; (n <= rslen) && (rs[n - 1] < maxbox); ir++) {
const long rw = (long)(minbox * pow(boxratio, (double)ir) + 0.5);
if (rw > rs[n - 1]) {
const long odd = force_odd(rw);
if (n < rslen) rs[n++] = odd;
else break;
}
}
if (n >= 1 && rs[n - 1] > maxbox) n--;
if (n < 1) n = 1;
*out_scales = rs;
return n;
}
/* -------------------------
Fast tokenizer: parse selected column
- Accepts separators: space/tab/comma
- Ignores non-numeric rows
------------------------- */
static int parse_column_inplace(const char *line, int selcol, double *out_value)
{
const unsigned char *p = (const unsigned char *)line;
int col = 0;
while (*p) {
while (*p == (unsigned char)' ' || *p == (unsigned char)'\t' || *p == (unsigned char)',') p++;
if (*p == '\0' || *p == '\n' || *p == '\r') break;
col++;
if (col == selcol) {
char *endptr = NULL;
errno = 0;
const double v = strtod((const char *)p, &endptr);
if (endptr == (const char *)p) return 0; /* no conversion */
if (errno != 0) return 0; /* overflow/underflow etc. */
*out_value = v;
return 1;
}
/* skip token */
while (*p &&
*p != (unsigned char)' ' && *p != (unsigned char)'\t' &&
*p != (unsigned char)',' &&
*p != (unsigned char)'\n' && *p != (unsigned char)'\r') {
p++;
}
}
return 0;
}
/* -------------------------
Context
------------------------- */
typedef struct {
/* data */
double *x; /* raw series */
double *y; /* working series (integrated demeaned by default) */
long n;
/* options */
int order; /* 0,2,4 */
int without_integration; /* -i */
int asymmetric_ends; /* -a (only d=0,2) */
long refresh_interval; /* -r (0 disables refresh) */
int use_scale_correction; /* -S */
double scale_correction; /* C(d) */
} DMAContext;
/* -------------------------
Load series (streaming, large-buffered)
------------------------- */
static long load_series_fast(DMAContext *ctx, const char *filename, int skip_lines, int selcol)
{
FILE *fp = fopen(filename, "r");
if (!fp) {
fprintf(stderr, "%s: cannot open file '%s'\n", g_progname, filename);
return 0;
}
/* large stdio buffer: reduces syscalls on huge files */
{
static char io_buf[1 << 20]; /* 1 MiB */
setvbuf(fp, io_buf, _IOFBF, sizeof(io_buf));
}
static char line_buf[1 << 20]; /* 1 MiB per line (safety for wide CSVs) */
for (int i = 0; i < skip_lines; i++) {
if (!fgets(line_buf, (int)sizeof(line_buf), fp)) break;
}
long cap = 1 << 14; /* 16384 initial */
long n = 0;
double *tmp = (double *)xmalloc((size_t)cap * sizeof(double));
while (fgets(line_buf, (int)sizeof(line_buf), fp)) {
double v;
if (!parse_column_inplace(line_buf, selcol, &v)) continue;
if (n >= cap) {
/* geometric growth to reduce realloc frequency */
cap = (cap < (1L << 20)) ? (cap * 2) : (cap + (1L << 20));
tmp = (double *)xrealloc(tmp, (size_t)cap * sizeof(double));
}
tmp[n++] = v;
}
fclose(fp);
if (n <= 0) {
free(tmp);
fprintf(stderr, "%s: no numeric data read from '%s'\n", g_progname, filename);
return 0;
}
ctx->x = (double *)xrealloc(tmp, (size_t)n * sizeof(double)); /* shrink-to-fit */
ctx->y = (double *)xmalloc((size_t)n * sizeof(double));
ctx->n = n;
if (!ctx->without_integration) {
/* demean and integrate */
double mean = 0.0;
for (long i = 0; i < n; i++) mean += ctx->x[i];
mean /= (double)n;
ctx->y[0] = ctx->x[0] - mean;
for (long i = 1; i < n; i++) ctx->y[i] = ctx->y[i - 1] + (ctx->x[i] - mean);
} else {
memcpy(ctx->y, ctx->x, (size_t)n * sizeof(double));
}
return n;
}
/* -------------------------
Center-trend coefficients by order
(computed per scale)
------------------------- */
typedef struct {
double c1, c2, c3; /* used as: a0 = c1*S0 + c2*S2 + c3*S4 (depending on order) */
} Coeff;
static Coeff coeff_dma0(long k)
{
Coeff c;
c.c1 = 1.0 / (double)(2 * k + 1);
c.c2 = 0.0;
c.c3 = 0.0;
return c;
}
static Coeff coeff_dma2(long k)
{
const double kk = (double)k;
const double k2 = kk * kk;
const double den = 8.0 * k2 * kk + 12.0 * k2 - 2.0 * kk - 3.0;
Coeff c;
c.c1 = (9.0 * k2 + 9.0 * kk - 3.0) / den;
c.c2 = -15.0 / den;
c.c3 = 0.0;
return c;
}
static Coeff coeff_dma4(long k)
{
const double kk = (double)k;
const double k2 = kk * kk;
const double k3 = k2 * kk;
const double k4 = k3 * kk;
const double k6 = k2 * k4;
const double den = 180.0 + 72.0 * kk - 800.0 * k2 - 320.0 * k3 + 320.0 * k4 + 128.0 * k6;
Coeff c;
c.c1 = (180.0 - 750.0 * kk - 525.0 * k2 + 450.0 * k3 + 225.0 * k4) / den;
c.c2 = (1575.0 - 1050.0 * kk - 1050.0 * k2) / den;
c.c3 = 945.0 / den;
return c;
}
/* -------------------------
Asymmetric ends (d=0,2)
------------------------- */
/* DMA0 end trend uses running average length that grows toward the center */
static double dma0_end_coeff(long m, long scale)
{
/* c1 = 1 / (((scale-1)/2) + m + 1) (original definition) */
return 1.0 / (double)(((scale - 1) / 2) + m + 1);
}
/* SG-like coefficients for DMA2 ends (copied from refined DMA2 implementation)
coeff array length: (2*k_half+1), indexed by (idx + k_half) for idx in [-k_half, k_half]
*/
static void sg2_fill(double *coeff, long m, long n, long k_half)
{
for (long i = 0; i < 2 * k_half + 1; i++) coeff[i] = 0.0;
const double mm = (double)m;
const double nn = (double)n;
double mn = 2.0 * mm + nn;
const double n2 = nn * nn;
const double n3 = n2 * nn;
const double m2 = mm * mm;
const double m3 = m2 * mm;
const double d = (mn + 3.0) * (mn - 3.0) * (mn - 1.0) * (mn + 1.0) * (mn + 5.0);
mn = mm * nn;
for (long idx = -m; idx <= k_half; idx++) {
const double kx = (double)idx;
const double kk2 = kx * kx;
const double val =
(90.0 + 720.0 * kk2 - 240.0 * mm + 576.0 * kx * mm + 960.0 * kk2 * mm +
48.0 * m2 + 1728.0 * kx * m2 + 960.0 * kk2 * m2 + 576.0 * m3 +
1152.0 * kx * m3 + 288.0 * m2 * m2 - 120.0 * nn - 432.0 * kx * nn -
960.0 * kk2 * nn - 528.0 * mn - 2304.0 * kx * mn - 1920.0 * kk2 * mn -
864.0 * m2 * nn - 2304.0 * kx * m2 * nn - 576.0 * m3 * nn +
84.0 * n2 + 576.0 * kx * n2 + 240.0 * kk2 * n2 + 720.0 * mm * n2 +
1152.0 * kx * mm * n2 + 720.0 * m2 * n2 - 72.0 * n3 - 144.0 * kx * n3 -
144.0 * mm * n3 + 18.0 * n2 * n2) / d;
if (idx >= -k_half && idx <= k_half) coeff[idx + k_half] = val;
}
}
/* -------------------------
Core estimators
------------------------- */
static double estimate_f2_dma0(const DMAContext *ctx, long scale, double *out_emax)
{
const double *restrict y = ctx->y;
const long n = ctx->n;
const long k = (scale - 1) / 2;
const Coeff cf = coeff_dma0(k);
const long nsteps = n - scale;
double f2 = 0.0;
double maxdiff = 0.0;
/* asymmetric ends (optional) */
if (ctx->asymmetric_ends) {
double ysumL = 0.0, ysumR = 0.0;
for (long i = 0; i < k; i++) {
ysumL += y[i];
ysumR += y[n - 1 - i];
}
for (long i = 0; i < k; i++) {
ysumL += y[k + i];
ysumR += y[n - 1 - k - i];
const double aL = dma0_end_coeff(i, scale) * ysumL;
const double aR = dma0_end_coeff(i, scale) * ysumR;
const double dL = y[i] - aL;
const double dR = y[n - 1 - i] - aR;
f2 += dL * dL + dR * dR;
}
}
/* center part (rolling sum) */
double S0 = 0.0;
for (long j = 0; j < scale; j++) S0 += y[j];
{
const double a0 = cf.c1 * S0;
const double d0 = y[k] - a0;
f2 += d0 * d0;
}
for (long p = 0; p < nsteps; p++) {
S0 += y[p + scale] - y[p];
if (ctx->refresh_interval > 0 && (p % ctx->refresh_interval) == 0) {
const double S0_old = S0;
double S0_new = 0.0;
const long base = (p + 1);
for (long j = 0; j < scale; j++) S0_new += y[base + j];
S0 = S0_new;
const double diff = fabs(S0_new - S0_old);
if (diff > maxdiff) maxdiff = diff;
}
const long center = (p + 1) + k;
const double a0 = cf.c1 * S0;
const double d = y[center] - a0;
f2 += d * d;
}
/* Emax in trend units */
if (out_emax) *out_emax = cf.c1 * maxdiff;
/* normalization consistent with refined DMA0 code */
if (ctx->asymmetric_ends) return f2 / (double)(nsteps + 1);
return f2 / (double)n;
}
static double estimate_f2_dma2(const DMAContext *ctx, long scale, double *out_emax)
{
const double *restrict y = ctx->y;
const long n = ctx->n;
const long k = (scale - 1) / 2;
const Coeff cf = coeff_dma2(k);
const long nsteps = n - scale;
double f2 = 0.0;
double maxdiff = 0.0;
/* asymmetric ends (optional) */
if (ctx->asymmetric_ends) {
double *coeff = (double *)xmalloc((size_t)(2 * k + 1) * sizeof(double));
for (long i = 1; i <= k; i++) {
const long m = i - 1;
sg2_fill(coeff, m, scale, k);
double aL = 0.0, aR = 0.0;
for (long j0 = 0; j0 <= m + k; j0++) {
const long cidx = (j0 - m);
const double w = coeff[cidx + k];
aL += w * y[j0];
aR += w * y[n - 1 - j0];
}
const double dL = y[m] - aL;
const double dR = y[n - 1 - m] - aR;
f2 += dL * dL + dR * dR;
}
free(coeff);
}
/* initial moments on first window */
double S0 = 0.0, S1 = 0.0, S2 = 0.0;
for (long j = 0; j < scale; j++) {
const double t = (double)(j - k);
const double v = y[j];
S0 += v;
S1 += v * t;
S2 += v * (t * t);
}
{
const double a0 = cf.c1 * S0 + cf.c2 * S2;
const double d0 = y[k] - a0;
f2 += d0 * d0;
}
for (long p = 0; p < nsteps; p++) {
const double x_out = y[p];
const double x_in = y[p + scale];
const double S0_old = S0;
const double S1_old = S1;
const double S2_old = S2;
/* rolling updates (refined DMA2 recurrence) */
S0 = S0_old - x_out + x_in;
S1 = S1_old - S0_old + x_out * (double)(k + 1) + x_in * (double)k;
S2 = S2_old - 2.0 * S1_old + S0_old
- x_out * (double)(k + 1) * (double)(k + 1)
+ x_in * (double)k * (double)k;
if (ctx->refresh_interval > 0 && (p % ctx->refresh_interval) == 0) {
const double S0r_old = S0;
const double S2r_old = S2;
double S0n = 0.0, S1n = 0.0, S2n = 0.0;
const long base = (p + 1);
for (long j = 0; j < scale; j++) {
const double t = (double)(j - k);
const double v = y[base + j];
S0n += v;
S1n += v * t;
S2n += v * (t * t);
}
S0 = S0n; S1 = S1n; S2 = S2n;
/* compare trend (not raw moments) to get a meaningful drift measure */
const double diff = fabs(cf.c1 * (S0n - S0r_old) + cf.c2 * (S2n - S2r_old));
if (diff > maxdiff) maxdiff = diff;
}
const long center = (p + 1) + k;
const double a0 = cf.c1 * S0 + cf.c2 * S2;
const double d = y[center] - a0;
f2 += d * d;
}
if (out_emax) *out_emax = maxdiff;
if (ctx->asymmetric_ends) return f2 / (double)(nsteps + 1);
return f2 / (double)n;
}
static double estimate_f2_dma4(const DMAContext *ctx, long scale, double *out_emax)
{
const double *restrict y = ctx->y;
const long n = ctx->n;
const long k = (scale - 1) / 2;
const Coeff cf = coeff_dma4(k);
const long nsteps = n - scale;
double maxdiff = 0.0;
/* initial moments */
double S0 = 0.0, S1 = 0.0, S2 = 0.0, S3 = 0.0, S4 = 0.0;
for (long j = 0; j < scale; j++) {
const double t = (double)(j - k);
const double v = y[j];
const double t2 = t * t;
const double t3 = t2 * t;
const double t4 = t2 * t2;
S0 += v;
S1 += v * t;
S2 += v * t2;
S3 += v * t3;
S4 += v * t4;
}
double f2;
{
const double a0 = cf.c1 * S0 + cf.c2 * S2 + cf.c3 * S4;
const double d0 = y[k] - a0;
f2 = d0 * d0;
}
const double kp1 = (double)(k + 1);
const double k0 = (double)k;
for (long p = 0; p < nsteps; p++) {
const double x_out = y[p];
const double x_in = y[p + scale];
/* store old values for recurrences and drift diagnostics */
const double S0_old = S0;
const double S1_old = S1;
const double S2_old = S2;
const double S3_old = S3;
const double S4_old = S4;
S0 = S0_old - x_out + x_in;
S1 = S1_old - S0_old + x_out * kp1 + x_in * k0;
S2 = S2_old - 2.0 * S1_old + S0_old
- x_out * (kp1 * kp1) + x_in * (k0 * k0);
S3 = S3_old - 3.0 * S2_old + 3.0 * S1_old - S0_old
+ x_out * (kp1 * kp1 * kp1) + x_in * (k0 * k0 * k0);
S4 = S4_old - 4.0 * S3_old + 6.0 * S2_old - 4.0 * S1_old + S0_old
- x_out * (kp1 * kp1 * kp1 * kp1) + x_in * (k0 * k0 * k0 * k0);
if (ctx->refresh_interval > 0 && (p % ctx->refresh_interval) == 0) {
double S0n = 0.0, S1n = 0.0, S2n = 0.0, S3n = 0.0, S4n = 0.0;
const long base = (p + 1);
for (long j = 0; j < scale; j++) {
const double t = (double)(j - k);
const double v = y[base + j];
const double t2 = t * t;
const double t3 = t2 * t;
const double t4 = t2 * t2;
S0n += v;
S1n += v * t;
S2n += v * t2;
S3n += v * t3;
S4n += v * t4;
}
S0 = S0n; S1 = S1n; S2 = S2n; S3 = S3n; S4 = S4n;
const double diff = fabs(cf.c1 * (S0n - S0_old) +
cf.c2 * (S2n - S2_old) +
cf.c3 * (S4n - S4_old));
if (diff > maxdiff) maxdiff = diff;
}
const long center = (p + 1) + k;
const double a0 = cf.c1 * S0 + cf.c2 * S2 + cf.c3 * S4;
const double d = y[center] - a0;
f2 += d * d;
}
if (out_emax) *out_emax = maxdiff;
/* DMA4 normalization consistent with refined DMA4 code */
return f2 / (double)(nsteps + 1);
}
static double estimate_f2(const DMAContext *ctx, long scale, double *out_emax)
{
const long n = ctx->n;
if (scale > n) die("Scale exceeds series length.");
if ((scale % 2) == 0) die("Scale must be odd.");
if (ctx->order == 0) {
if (scale < 5) die("For d=0, scale must be >= 5.");
return estimate_f2_dma0(ctx, scale, out_emax);
} else if (ctx->order == 2) {
if (scale < 7) die("For d=2, scale must be >= 7.");
return estimate_f2_dma2(ctx, scale, out_emax);
} else if (ctx->order == 4) {
if (scale < 9) die("For d=4, scale must be >= 9.");
return estimate_f2_dma4(ctx, scale, out_emax);
}
die("Unsupported order. Use -d 0, -d 2, or -d 4.");
return 0.0;
}
/* -------------------------
Main
------------------------- */
int main(int argc, char **argv)
{
g_progname = (argc > 0 && argv[0]) ? argv[0] : "sgdma";
/* CLI defaults */
int skip_lines = 0;
int selcol = 1;
long minbox = -1; /* order-dependent if not specified */
long maxbox = 0;
const double boxratio = pow(2.0, 1.0 / 8.0);
DMAContext ctx;
memset(&ctx, 0, sizeof(ctx));
ctx.order = 0;
ctx.use_scale_correction = 1;
ctx.without_integration = 0;
ctx.asymmetric_ends = 0;
ctx.refresh_interval = -1; /* sentinel: not set */
if (argc < 2) {
print_usage();
return EXIT_FAILURE;
}
/* parse options */
int argi = 1;
while (argi < argc && argv[argi][0] == '-') {
const char *opt = argv[argi];
if (strcmp(opt, "-h") == 0) {
print_usage();
return EXIT_SUCCESS;
} else if (strcmp(opt, "-d") == 0) {
if (++argi >= argc) die("Missing value after -d.");
ctx.order = atoi(argv[argi]);
if (!(ctx.order == 0 || ctx.order == 2 || ctx.order == 4)) die("ORDER must be 0, 2, or 4.");
} else if (strcmp(opt, "-S") == 0) {
if (++argi >= argc) die("Missing value after -S.");
ctx.use_scale_correction = (atoi(argv[argi]) != 0);
} else if (strcmp(opt, "-s") == 0) {
if (++argi >= argc) die("Missing value after -s.");
skip_lines = atoi(argv[argi]);
if (skip_lines < 0) skip_lines = 0;
} else if (strcmp(opt, "-c") == 0) {
if (++argi >= argc) die("Missing value after -c.");
selcol = atoi(argv[argi]);
if (selcol < 1) die("COLUMN must be >= 1.");
} else if (strcmp(opt, "-l") == 0) {
if (++argi >= argc) die("Missing value after -l.");
minbox = atol(argv[argi]);
} else if (strcmp(opt, "-u") == 0) {
if (++argi >= argc) die("Missing value after -u.");
maxbox = atol(argv[argi]);
} else if (strcmp(opt, "-a") == 0) {
ctx.asymmetric_ends = 1;
} else if (strcmp(opt, "-i") == 0) {
ctx.without_integration = 1;
} else if (strcmp(opt, "-r") == 0) {
if (++argi >= argc) die("Missing value after -r.");
ctx.refresh_interval = atol(argv[argi]);
if (ctx.refresh_interval < 0) ctx.refresh_interval = 0;
} else {
print_usage();
die("Unknown option.");
}
argi++;
}
if (argi >= argc) {
print_usage();
die("Missing input filename.");
}
const char *filename = argv[argi];
/* apply order-dependent defaults */
if (minbox < 0) minbox = default_minbox_for_order(ctx.order);
minbox = force_odd(minbox);
if (ctx.refresh_interval < 0) {
ctx.refresh_interval = default_refresh_for_order(ctx.order);
}
/* note: -r 0 explicitly disables refresh */
ctx.scale_correction = default_scale_correction_for_order(ctx.order);
/* for d=4: ignore -a but accept for CLI compatibility */
if (ctx.order == 4 && ctx.asymmetric_ends) {
fprintf(stderr, "%s: note: -a is ignored for d=4.\n", g_progname);
ctx.asymmetric_ends = 0;
}
/* load data */
const long n = load_series_fast(&ctx, filename, skip_lines, selcol);
if (n <= 0) return EXIT_FAILURE;
fprintf(stderr, "Input file: %s\n", filename);
fprintf(stderr, "Total data length = %ld\n", n);
fprintf(stderr, "***** centered DMA (d=%d) *****\n", ctx.order);
/* maxbox default and validation */
const long min_required = default_minbox_for_order(ctx.order);
if (minbox < min_required) minbox = min_required;
if (maxbox <= 0 || maxbox > n / 4 || minbox > maxbox) {
maxbox = n / 4;
}
if (maxbox < minbox) maxbox = minbox;
/* scales */
long *scales = NULL;
const long scales_len = generate_scales_exact(minbox, maxbox, boxratio, &scales);
if (scales_len <= 0) die("Failed to generate scales.");
double *F2 = (double *)xmalloc((size_t)scales_len * sizeof(double));
double *Emax = (double *)xmalloc((size_t)scales_len * sizeof(double));
/* timing */
const clock_t t0 = clock();
for (long i = 0; i < scales_len; i++) {
double emax = 0.0;
F2[i] = estimate_f2(&ctx, scales[i], &emax);
Emax[i] = emax;
}
const clock_t t1 = clock();
const double ms = 1000.0 * (double)(t1 - t0) / (double)CLOCKS_PER_SEC;
fprintf(stderr, "Processing time: %.3f [ms]\n", ms);
/* output */
for (long i = 0; i < scales_len; i++) {
const double s = (double)scales[i];
const double s_eff = (ctx.use_scale_correction ? (s / ctx.scale_correction) : s);
/* use log10() directly (faster and clearer than log()/log(10)) */
const double log10s = log10(s_eff);
const double log10F = 0.5 * log10(F2[i]); /* log10(sqrt(F2)) */
fprintf(stdout, "%.15g,%.15g,%.15g\n", log10s, log10F, Emax[i]);
}
/* cleanup */
free(ctx.x);
free(ctx.y);
free(scales);
free(F2);
free(Emax);
return EXIT_SUCCESS;
}
1. プログラムの概要
上で与えたコードでは、
- centered DMA のゆらぎ関数
を、スケール
ごとに計算
- トレンド除去次数を d=0 / d=2 / d=4 から選択
- 大きな入力でも速度が落ちにくい
出力は 3列（コンマ区切り）：
- 第1列：
- 第2列：
- 第3列：高速アルゴリズムの累積誤差の評価結果
- 第1列：
という処理をします。DMAについては、別の記事を参照してください。
とにかにDMAを使いたい人は、まず、Cプログラムをコンパイルしてください。
本プログラムは R や Python のようなインタプリタ言語で書かれたスクリプトではなく、C 言語で書かれています。そのため、ソースコード（.c ファイル）をそのまま実行することはできません。あらかじめ コンパイラを用いて、CPU が直接実行できる 実行ファイル（sgdma や sgdma.exe）に変換する必要があります。
R や Python では、コードをそのまま読み込んで逐次解釈・実行しますが、C 言語では、
「書く → コンパイルする → 実行する」という手順を踏みます。このひと手間がある代わりに、C プログラムは高速で、長い時系列や大量データでも安定して処理できるという利点があります。
一度コンパイルしてしまえば、その後は 実行ファイルをコマンドとして繰り返し使うだけです。DMA 解析を本格的に回したい場合や、長時間相関解析をバッチ処理したい場合には、この点が大きなメリットになります。
2. コンパイル方法
本プログラムは 標準 C（C99 相当）で書かれており、Linux / macOS / Windows のいずれでも、一般的な C コンパイラでコンパイルできます。ここでは代表的な gcc と Visual C++（cl.exe） の例を示します。
2.1 gcc（Linux / macOS / WSL など）
最小コマンドで、とにかく動かしたい場合は、コマンドラインで、
gcc sgdma.c -lm -o sgdma
を実行してください。これだけでコンパイルできます。
-lm：数学ライブラリのリンク（pow,logなどを使用）
より高速で、安全にコンパイルするための推奨コマンドは、
gcc -O3 -std=c99 -Wall -Wextra -pedantic sgdma.c -lm -o sgdma
です。
-O3：最適化（高速化）-std=c99：C99 規格を指定-Wall -Wextra -pedantic：警告をできるだけ厳しく出す（推奨）-o sgdma：出力ファイル名を指定
補足
clangを使う場合も、ほぼ同じオプションでコンパイルできます。
2.2 Visual C++（Windows / cl.exe）
Windows 環境では、Visual Studio に付属する Visual C++ コンパイラ（cl.exe） を使用できます。Developer Command Prompt for Visual Studio（または Developer PowerShell）を起動して実行してください。
最小コマンドで、とにかく動かしたい場合は、コマンドラインで、
cl sgdma.c
を実行してください。これだけで、コンパイルできて、
sgdma.exe
が生成されます。
- 数学関数（
pow,logなど）のための特別なライブラリ指定は 不要です。 - 最適化や警告設定を指定しなくても、標準設定で問題なく動作します。
より高速で、安全にコンパイルするための推奨コマンドは、
cl /O2 /std:c11 /W4 sgdma.c /Fe:sgdma.exe
です。
/O2：最適化（高速化）/std:c11：C11 モード （Visual C++ では C99 より C11 指定が一般的で、C99 相当の機能も含まれます）/W4：警告レベルを高める/Fe:sgdma.exe：出力ファイル名を指定
cl はDeveloper Command Prompt から実行するのが確実です。
3. 入力ファイル形式
解析する時系列データの準備（データファイルの作成）について説明します。データは、テキストファイル（CSV / 空白区切り / タブ区切り）で、1 行につき 1 サンプルにしてください。複数の時系列を複数列にまとめて保存しておいても構いません。入力ファイルに、列が複数ある場合、解析する際に目的の列を -c で指定できます。ヘッダなど数値でない行が先頭にあっても構いません。そのようなヘッダ行は、-s で指定すればスキップできます。
入力データ（ファイルの中身）は、以下のような感じです。
time,signal 0.00,0.12 0.01,0.15 0.02,0.11 ...
この場合、解析したい列が、signal の 2 列目なので -c 2 を指定します。ヘッダ行を１行読み込まないので -s 1 とします（省略しても自動的にスキップしますが）。
4. 最小の実行例
解析では以下のコマンドを実行します。ここでは、ファイル "data.csv" に、１列の時系列データが含まれており、その時系列を入力データとして解析して、DMAの結果を "out.csv" に出力します。
./sgdma data.csv > out.csv
Windowsのコマンドラインでは、
sgdma.exe data.csv > out.csv
とします。
さらなる応用として、入力データのヘッダ 1 行を飛ばして 2 列目を解析する場合は、
./sgdma -s 1 -c 2 data.csv > out.csv
Windowsのコマンドラインでは、
sgdma.exe -s 1 -c 2 data.csv > out.csv
とします。
5. オプションの意味と選び方
sgdma のオプションは次の通りです。
5.1 -d ORDER：トレンド除去の次数
-d 0 # DMA0（移動平均のみ） -d 2 # DMA2（2次トレンド除去） -d 4 # DMA4（4次トレンド除去）
目安：
- d=0：最も単純。オリジナルのDMAに対応します。
- d=2：実務で使いやすい標準的な選択肢。3次多項式で記述されるトレンドまで除去できます。
- d=4：5次多項式で記述される高次のトレンドまで除去可能。短い時系列（1000点以下）や短時間のスケーリング領域（100点以下）の解析には向きません。
たとえば、DMA2 を使う場合、
./sgdma -d 2 data.csv > out.csv
あるいは、Windowsのコマンドラインでは、
sgdma -d 2 data.csv > out.csv
とします。
5.2 -S FLAG：スケール補正を使うか
このプログラムは、次数 d に応じたスケール補正をデフォルトで適用します。
出力第1列は ですが、ここで
-S 1（デフォルト）：-S 0：
です。
定数はコード中の既定値：
補正なしで出したいときは、
./sgdma -d 2 -S 0 data.csv > out.csv
あるいは、Windowsのコマンドラインでは、
sgdma -d 2 -S 0 data.csv > out.csv
とします。
5.3 -i：積分をしない（入力をそのまま使う）
DMAのデフォルトでは、「平均除去後に積分（累積和）」を行い、解析対象系列 を作ります。しかし、非整数ブラウン運動のような拡散的時系列を分析するときは、積分（累積和）する必要はありません（ただし、スケーリング指数の解釈には注意してください）。この積分の有無は、
- デフォルト（積分あり）：
-iを付けない - 積分しない：
-iを付ける（入力をそのまま解析）
で設定できます。
積分なしで解析する場合は、
./sgdma -i data.csv > out.csv
とします。
5.4 -l MINBOX と -u MAXBOX：スケール範囲
解析するスケールの下限、上限の指定のオプションは、
-l：最小スケール（窓幅）-u：最大スケール（窓幅）
です。いずれの値も 奇数に強制されます。デフォルトの最小値は d に依存します：
- d=0 → 5
- d=2 → 7
- d=4 → 9
最大スケールはデフォルト N/4（N はデータ長）です。
4次のDMAで、解析するスケールを 31 から 5000 に設定する場合は、
./sgdma -d 4 -l 31 -u 5000 data.csv > out.csv
5.5 -r REFRESH：rolling update の再計算間隔（数値ドリフト対策）
このプログラムは、スライド窓の更新を O(1) で行います（高速）。ただし、非常に長い系列では、加減算の繰り返しにより浮動小数のドリフトが蓄積することがあります。-r はそれを抑えるために、一定間隔で総和を再計算します。
デフォルト値（d に依存）：
- d=0 → 2000
- d=2 → 200
- d=4 → 20
d=2 でより安定側に振りたい場合は、デフォルトの -r 200 を変更し、
./sgdma -d 2 -r 50 data.csv > out.csv
のように -r 50 の値を小さくしてください。
DMA結果の両対数プロットが、単調増加ではない場合は、数値的な丸目誤差の蓄積による計算ミスですので、-r の値を小さくしてみてください。両対数プロットが、単調増加になれば問題は解決です。
5.6 -a：端点の非対称処理（d=0,2 のみ）
-a を付けると、端点（両端）についても非対称窓でトレンド推定を行い、端点寄与を含めた を計算します。
- d=0,2：有効
- d=4：現在の実装では無効（指定しても無視され、stderr に注意が出ます）
時系列の端点まで分析したいとき（短い時系列を解析するとき）は、
./sgdma -d 2 -a data.csv > out.csv
を試してください。
6. 出力の見方
出力は常にコンマ区切り（CSV）の 2 列です。
log10(s_eff),log10(F(s))
- 第1列：
（
-S 1の場合は (s/C(d))、-S 0の場合は (s)） - 第2列：
この 2 列を両対数プロットし、直線領域を回帰すればスケーリング指数を推定できます。
7. このブログを参考にした研究成果の発表について
ここに掲載したプログラム（または改変版）を用いて論文等を書く場合は、以下の論文を引用してください。
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).
※ もし記事の中で「ここ違うよ」という点や気になるところがあれば、気軽に指摘していただけると助かります。質問や「このテーマも取り上げてほしい」といったリクエストも大歓迎です。必ず対応するとは約束できませんが、できるだけ今後の記事で扱いたいと思います。それと、下のはてなブログランキングはあまり信用できる指標ではなさそうですが (私のブログを読んでいる人は、実際とても少ないです)、押してもらえるとシンプルに励みになります。気が向いたときにポチッとしていただけたら嬉しいです。
