今回から、間欠性（Intermittency）、あるいは 分散不均一性（Heteroscedasticity） と呼ばれる性質をもつ非ガウス時系列を対象に、多重スケール確率分布解析を用いて評価する実務的な方法を解説していきます。
ここで扱う時系列は、対数正規カスケード過程により、ガウスノイズの局所的な変動幅（振幅）を変調したものです。変動幅の変調とは、下図に示すように、時系列の上下方向の揺らぎの大きさを、場所によって広げたり狭めたりすることを意味します。下図では、部分区間をステップごとに 2 分割し、区間の数を
と増やしています。ただし、この図が示しているのは、あくまで時系列の局所標準偏差の大きさのみです。実際の解析では、このように局所標準偏差が時間とともに変化するガウス過程を、数値的に生成した時系列データを用います。
ここでの問題意識は、下図に示す左右の非ガウス過程の違いを、どのように定量的に特徴づけるかという点にあります。非ガウス過程は「弱定常」の枠組みでは適切に記述できないため、古典的なパワースペクトル解析は有効な手段とはなりません。
視覚的には、下図の左右の時系列が異なる性質をもっていることは、直感的に理解できると思います。しかしながら、両者の パワースペクトル も、各時刻の値 の 確率分布 も、実際には完全に一致しています。このような状況で役に立つのが「多重スケール確率分布解析」（multiscale PDF analysis）です。多重スケール確率分布解析では、時系列を粗視化（局所平均化）し、上のカスケードステップを逆に（細かい構造から粗い構造へ）観察することで、分布の非ガウス性の変化を評価します。
の分布はまったく同じ非ガウス分布になっている。
局所的な粗視化（局所平均化）と、値のばらつきの非一様性との関係を直感的に理解するために、身近な例として高校のテスト成績を考えてみましょう。
同じ学年の生徒が同一のテストを受け、その結果をクラスごとに集計するとします。高校によっては、前年までの成績に基づいてクラス分けを行う場合もあれば、成績とは無関係に、ほぼランダムに生徒をクラスに分ける場合もあります。いずれの場合でも、テストの点数は高校全体としても、各クラス内でもばらつきます。多くの生徒は平均的な点数を取りつつ、数人は飛び抜けて高得点、あるいは低得点を取るでしょう。ここでは、高校全体としての点数分布は、どちらの方法でクラス分けした場合でもほぼ同じであると仮定します（教育効果の違いなどは無視できるとします）。
このとき、40 人程度のクラスごとに 平均点（局所平均） を計算すると、状況は大きく変わります。過去の成績に基づいてクラス分けをしている高校では、成績の良いクラスと悪いクラスの差が、そのままクラス平均の違いとして現れます。すなわち、クラスごとに見た平均値のばらつきは大きくなります。
一方で、生徒をランダムにクラス分けしている高校では、どのクラスにも成績の良い生徒と悪い生徒が混在するため、クラス平均はどこも似たような値に収束します。クラス内の点数分布が大きくばらついていたとしても、クラス平均のばらつきは小さくなり、中心極限定理により、クラス平均の分布は正規分布に従うことが期待されます。
この違いは、局所的に似た性質のデータが集まっているか（クラスターが存在するか）、それとも 局所的に異なるデータが混在しているか の違いに対応しています。非ガウス性や分散の非一様性をもつ時系列では、まさにこのように「場所（時間）ごとにばらつきの大きさが異なる」構造が現れます。
粗視化とは、このような局所平均の振る舞いを通して、データの中に隠れたクラスター構造や変動の偏りが存在するかどうかを見極める操作と理解することができます。
多重スケール確率分布解析の流れ
(1) 粗視化
多重スケール確率分布解析では、解析したい時系列を粗視化（局所平均化）します。このときに、時系列に含まれるトレンド成分を取り除くことが不可欠です。以下で与えるRスクリプトには、トレンド除去の処理が既に含まれていますが、今回はトレンド除去の詳しい説明は省略します。
粗視化（局所平均化）は、窓幅（局所のデータ点数）を表すスケール s で、局所あるいは局所平均を計算することです。通常は、窓を１点ずつスライドさせますが、下図では窓（あるいはブロック）のサイズを見やすくするために、オバーラップなしで区間を区切って粗視化しています。
(2) 粗視化時系列を標準化して非ガウスパラメタを推定
次に、粗視化時系列を平均０、分散１に標準化して、非ガウスパラメタ を計算します。
は、以前に紹介した相乗対数正規モデルの非ガウス性を表すパラメタで、
のとき正規（ガウス）分布、
で大きくなるほど、裾が厚く、中央が尖った非ガウス分布を意味します。
下図右上は、実際に間欠性ゆらぎのサンプル時系列に対して、スケール s を変化させ、非ガウスパラメタ を推定した結果です。そして下段は、粗視化時系列の分布（薄青実線）に、相乗対数正規モデルの分布（青破線）を重ねて描いたものです。このグラフでは縦軸が対数スケールになっていますので、正規（ガウス）分布は２次関数（黒破線）になり、相乗対数正規モデルの分布はテント型の形状に見えます。
ここで使った時系列は対数正規カスケード過程ですので、理論的に、
となり、この依存性が上図右上の破線で描かれています。実際の粗視化時系列も理論に近い依存性になっています。さらに、下段に示した分布形状も、モデルのものと良く一致しています。
のスケール依存性（破線は理論値）。下段：
での粗視化時系列の確率密度関数の推定結果（薄青実線）に対して相乗対数正規モデルをフィットしたもの（青破線）と正規（ガウス）分布（黒破線）。縦軸が対数スケールになっていることに注意。縦軸が対数スケールの場合、正規（ガウス）分布は２次関数になる。
とにかくRでやってみる
数式をたくさん書いて理論を説明しても、学生がやる気をなくすので、まずは数値的に生成したサンプル時系列を使って多重スケール確率分布解析を試せるようにしておきます。 とはいえ、最低限のポイントは言っておきます。
- 平均、分散ではなく、分布の形を見たいので粗視化したデータは平均０、分散１に標準化
- 非ガウス分布のモデルとして、今回は「相乗対数正規分布」を適用し、モデルの非ガウスパラメタ
で非ガウス性を定量化
- トレンド成分がない時系列でも、トレンド除去をした方が確率密度関数の推定結果は改善
これらのポイントを押さえた解析を行うRスクリプトの例が以下です。ここでは、対数正規カスケード過程のサンプル時系列を数値的に生成して、それを分析しています。実行すると上の図と同じものが出てくるはずです。
############################################################
# Simulation of an intermittent (non-Gaussian) time series
# and multiscale analysis (SG-detrended data only)
#
# Objectives (analyze SG-detrended series only):
# 1) Generate x[n] = (lognormal cascade weights) * (white noise)
# 2) Center x[n] by the median and build the integrated series y[n]
# 3) Apply Savitzky–Golay (SG) detrending (window length = scale s),
# then crop endpoints
# 4) Compute increments Δ_s y[n] and standardize them
# 5) Estimate λ^2 from q-th order moments (SG only)
# 6) Plot PDFs: empirical KDE + Castaing-type model PDF (dashed),
# and a standard Gaussian PDF (gray dashed)
#
# Plot layout:
# - Top row: 2 panels (left: x[n] spanning 2 columns, right: λ^2(s))
# - Bottom row: 3 panels (PDFs at 3 scales)
# - KDE: light blue with alpha
# - Model: blue dashed (lty=2)
# - Gaussian: gray dashed (lty=2)
############################################################
# -----------------------------
# 0) Package
# -----------------------------
if (!requireNamespace("signal", quietly = TRUE)) install.packages("signal")
library(signal)
# -----------------------------
# 1) Colors
# -----------------------------
col_kde <- rgb(0, 0.5, 1, alpha = 0.3) # light blue (KDE)
col_mod <- 4 # blue (model)
# -----------------------------
# 2) Utility: standardize to mean 0 and variance 1
# -----------------------------
standardize <- function(z) {
z <- z[is.finite(z)]
if (length(z) < 5) return(rep(NA_real_, length(z)))
mu <- mean(z)
sdv <- sd(z)
if (!is.finite(sdv) || sdv <= 0) return(rep(NA_real_, length(z)))
(z - mu) / sdv
}
# -----------------------------
# 3) Increment Δ_s y[n] = y[n+s] - y[n] (safe version)
# -----------------------------
diff_by_scale <- function(y, s) {
n <- length(y)
s <- as.integer(s)
if (!is.finite(s) || s < 1L) return(numeric(0))
if (n <= s) return(numeric(0))
y[(s + 1L):n] - y[1L:(n - s)]
}
# -----------------------------
# 4) SG detrending and endpoint cropping
# -----------------------------
sg_detrend_and_crop <- function(y, sg_order, sg_window) {
sg_window <- as.integer(sg_window)
if (sg_window %% 2 == 0) stop("sg_window must be odd.")
if (sg_window <= sg_order) stop("sg_window must be > sg_order.")
half <- (sg_window - 1L) / 2L
y_tr <- sgolayfilt(y, p = sg_order, n = sg_window, m = 0, ts = 1)
y_det <- y - y_tr
if (length(y_det) <= 2L * half + 1L) return(numeric(0))
y_det[(half + 1L):(length(y_det) - half)]
}
# -----------------------------
# 5) Generator: lognormal cascade weights × white noise
# -----------------------------
generate_cascade_series <- function(lambda2_total, n_level, n_base) {
lambda2_per_level <- lambda2_total / n_level
sigma_logw <- sqrt(lambda2_per_level)
weight <- rep(1.0, n_base)
for (lev in seq_len(n_level)) {
n_parent <- length(weight)
# logW ~ N(-lambda2_per_level, lambda2_per_level)
logw <- rnorm(2L * n_parent, mean = -lambda2_per_level, sd = sigma_logw)
w_left <- weight * exp(logw[1L:n_parent])
w_right <- weight * exp(logw[(n_parent + 1L):(2L * n_parent)])
weight <- as.numeric(rbind(w_left, w_right))
}
rnorm(length(weight)) * weight
}
# -----------------------------
# 6) Estimate lambda^2 from q-th order moments
# -----------------------------
estimate_lambda2_from_q <- function(z, q) {
z <- z[is.finite(z)]
if (length(z) < 10) return(NA_real_)
M_const <- (-digamma(1) + log(2)) / 2
az <- abs(z)
eps <- .Machine$double.eps
az <- pmax(az, eps)
if (abs(q) < 1e-3) {
return(-mean(log(az), na.rm = TRUE) - M_const)
} else if (abs(q - 2) < 1e-3) {
return(mean(z^2 * log(az), na.rm = TRUE) + M_const - 1)
} else {
mq <- mean(az^q, na.rm = TRUE)
if (!is.finite(mq) || mq <= 0) return(NA_real_)
2 / (q * (q - 2)) * (
log(pi) / 2 + log(mq) - log(2) * q / 2 - lgamma((q + 1) / 2)
)
}
}
# -----------------------------
# 7) KDE (density) on a common grid
# -----------------------------
kernel_pdf_on_grid <- function(z, x_grid,
bw = "nrd0", adjust = 1,
kernel = "gaussian", n = 2048) {
z <- z[is.finite(z)]
if (length(z) < 10) return(rep(NA_real_, length(x_grid)))
d <- density(z, bw = bw, adjust = adjust, kernel = kernel,
from = min(x_grid), to = max(x_grid), n = n, na.rm = TRUE)
approx(d$x, d$y, xout = x_grid, rule = 2, yleft = 0, yright = 0)$y
}
# -----------------------------
# 8) Castaing-type model PDF for standardized z
# -----------------------------
castaing_pdf_std <- function(z_grid, lambda2, n_u = 2001, n_sd = 8) {
if (!is.finite(lambda2) || lambda2 <= 0) return(rep(NA_real_, length(z_grid)))
sd_u <- sqrt(lambda2)
mu_u <- -lambda2
u_min <- mu_u - n_sd * sd_u
u_max <- mu_u + n_sd * sd_u
u <- seq(u_min, u_max, length.out = n_u)
du <- u[2] - u[1]
f_u <- dnorm(u, mean = mu_u, sd = sd_u)
coef_u <- (1 / exp(u)) * f_u * (1 / sqrt(2 * pi))
exp_factor_u <- exp(-2 * u)
w <- rep(1, n_u)
w[c(1, n_u)] <- 0.5
pdf <- numeric(length(z_grid))
z2 <- z_grid^2
for (i in seq_along(z_grid)) {
pdf[i] <- sum(w * coef_u * exp(-0.5 * z2[i] * exp_factor_u)) * du
}
pdf
}
# -----------------------------
# 9) Plot PDF: KDE + model + Gaussian (SG only), y-axis as 10^m
# -----------------------------
plot_pdf_kernel_with_model_sg <- function(dy_sg_std, lam2_sg, s,
xlim = c(-7, 7),
ylim = c(1e-4, 1),
bw = "nrd0", adjust = 1,
kernel = "gaussian",
n_grid = 1200,
col_kde, col_mod) {
x_grid <- seq(xlim[1], xlim[2], length.out = n_grid)
y_kde <- kernel_pdf_on_grid(dy_sg_std, x_grid,
bw = bw, adjust = adjust, kernel = kernel)
y_mod <- castaing_pdf_std(x_grid, lam2_sg)
# Standard Gaussian N(0,1)
y_gauss <- dnorm(x_grid, mean = 0, sd = 1)
# Avoid zeros for log scale
eps <- 1e-300
y_kde <- pmax(y_kde, eps)
y_mod <- pmax(y_mod, eps)
y_gauss <- pmax(y_gauss, eps)
# Plot (custom y-axis)
plot(x_grid, y_kde, type = "l", log = "y",
col = col_kde, lwd = 2,
xlim = xlim, ylim = ylim,
main = paste0("scale = ", s),
xlab = "Standardized increment Δ_s y / SD(Δ_s y)",
ylab = "Probability density",
yaxt = "n")
# Gaussian (gray dashed)
lines(x_grid, y_gauss, col = "gray50", lwd = 2, lty = 2)
# Model (blue dashed)
lines(x_grid, y_mod, col = col_mod, lwd = 2, lty = 2)
# y-axis ticks: 10^m
m_min <- floor(log10(ylim[1]))
m_max <- ceiling(log10(ylim[2]))
m_vec <- m_min:m_max
at_y <- 10^m_vec
axis(2,
at = at_y,
labels = as.expression(sapply(m_vec, function(m) bquote(10^.(m)))),
las = 1)
legend("topright",
legend = c("Empirical PDF",
"Gaussian",
"Model",sprintf("(lambda^2 = %.3f)", lam2_sg)),
col = c(col_kde, "gray50", col_mod,col_mod),
lwd = c(2,2,2,NA),
lty = c(1, 2, 2,NA),
bty = "n")
}
############################################################
# 10) Main: parameter settings
############################################################
# --- Generation parameters
lambda2_target <- 0.9
n_level <- 12
n_base <- 16
# --- Analysis parameters
q_moment <- 0.25
sg_order <- 2
# --- Scale list (odd; also used as SG window length)
scale_list <- floor(exp(seq(log(5), log(2^n_level), length.out = 20)) / 2) * 2 + 1
scale_list <- unique(as.integer(scale_list))
scale_list <- scale_list[scale_list >= 5]
# --- Scales for PDF panels (edit these if needed)
scale_show <- c(9, 99, 999)
############################################################
# 11) Run: generate -> integrate -> estimate lambda^2 (SG only) -> plot
############################################################
set.seed(1)
x_raw <- generate_cascade_series(lambda2_target, n_level, n_base)
N <- length(x_raw)
# Center by median (recommended)
x_centered <- x_raw - median(x_raw, na.rm = TRUE)
# Integrate
y_int <- cumsum(x_centered)
# Estimate lambda^2 across scales (SG only)
lambda2_with_sg <- rep(NA_real_, length(scale_list))
for (k in seq_along(scale_list)) {
s <- scale_list[k]
# SG detrend (window length = s)
y_det <- sg_detrend_and_crop(y_int, sg_order = sg_order, sg_window = s)
dy <- diff_by_scale(y_det, s)
dy_std <- standardize(dy)
lambda2_with_sg[k] <- estimate_lambda2_from_q(dy_std, q_moment)
}
# Reference theoretical curve
theory_lambda2 <- function(s) {
lambda2_target * (1 - log(s) / log(2) / n_level)
}
ylim_top <- lambda2_target * (1 - log(scale_list[1]) / log(2) / n_level)
# Get lambda^2 at a specific scale (or NA if missing)
get_lambda2_at_scale <- function(s, scale_list, lambda2_vec) {
idx <- which(scale_list == s)
if (length(idx) == 1) return(lambda2_vec[idx])
NA_real_
}
############################################################
# 12) Plot (top row: 2 panels; bottom row: 3 panels; no blank panel)
############################################################
layout(matrix(c(1, 1, 2,
3, 4, 5),
nrow = 2, byrow = TRUE))
par(las = 1, xaxs = "i", mar = c(5, 5, 2, 2),cex.axis=1.2,cex.lab=1.5)
# (1) Generated series x[n] (spans 2 columns)
x_max <- max(abs(x_raw), na.rm = TRUE)
plot(0:(N - 1), x_raw, type = "l", col = 4,
xlab = "n", ylab = "x[n]",
ylim = c(-x_max, x_max),
main = "Generated series x[n]")
# (2) lambda^2 vs scale (SG only)
plot(scale_list, lambda2_with_sg,
col = 4, pch = 16, log = "x",
xlab = "Coarse-graining scale s (window length)", ylab = expression(lambda^2),
main = sprintf("Scale dependence of non-Gaussianity (q = %.2f)", q_moment),
ylim = c(0, ylim_top))
lines(scale_list, theory_lambda2(scale_list), lty = 2)
abline(h=0,lty=2)
# (3)-(5) PDFs (KDE + Gaussian + model)
for (s in scale_show) {
y_det <- sg_detrend_and_crop(y_int, sg_order = sg_order, sg_window = s)
dy <- diff_by_scale(y_det, s)
dy_sg_std <- standardize(dy)
lam2_sg <- get_lambda2_at_scale(s, scale_list, lambda2_with_sg)
if (!is.finite(lam2_sg)) lam2_sg <- estimate_lambda2_from_q(dy_sg_std, q_moment)
plot_pdf_kernel_with_model_sg(
dy_sg_std = dy_sg_std,
lam2_sg = lam2_sg,
s = s,
xlim = c(-7, 7), ylim = c(1e-4, 1),
bw = "nrd0", adjust = 1, kernel = "gaussian",
col_kde=col_kde, col_mod=col_mod
)
}
間欠性のないIID過程との比較
間欠性の特徴を定量化できるのが多重スケール確率分布解析ですので、比較のために間欠性がない独立同分布（IID）過程の分析結果も示しておきます。上のRスクリプトを実行した後に、下のRスクリプトを実行すれば、下図が描かれると思います。ここでは、間欠性ゆらぎのサンプル時系列の順番をシャッフルして、IID過程のサンプルを作成しています。このことは、 の分布が、元の間欠性ゆらぎでも、IID過程でも全く同じであることを意味します。
間接性を示す対数正規カスケード過程のサンプル時系列と比べれば、IID過程の粗視化時系列が急速に正規（ガウス）分布に収束することがわかります（下図右上）。したがって、両者の違いは、非ガウスパラメタのスケール依存性を見れば、はっきりとわかります。
のスケール依存性（青は間欠性ゆらぎ、赤はIID過程）。下段：
での粗視化時系列の確率密度関数の推定結果（薄青実線）に対して相乗対数正規モデルをフィットしたもの（青破線）と正規（ガウス）分布（黒破線）。縦軸が対数スケールになっていることに注意。縦軸が対数スケールの場合、正規（ガウス）分布は２次関数になる。
############################################################
# Shuffled control analysis (starting from an existing x_raw)
#
# Purpose:
# After the original analysis has been completed for x_raw,
# shuffle x_raw (destroy temporal correlations while preserving
# the marginal distribution) and repeat the same SG-based
# multiscale analysis. Then compare the results with the original.
#
# Assumptions (already defined in the workspace):
# - x_raw : original generated series
# - scale_list, scale_show : scale settings
# - sg_order, q_moment : analysis parameters
# - lambda2_with_sg : lambda^2(s) from the original series (SG only)
# - Functions:
# sg_detrend_and_crop(), diff_by_scale(), standardize(),
# estimate_lambda2_from_q(), plot_pdf_kernel_with_model_sg()
#
# Output:
# A 2x2+ layout:
# (1) Shuffled x[n] (spans 2 columns)
# (2) lambda^2(s): original vs shuffled
# (3)-(5) PDFs (shuffled) at 3 selected scales
############################################################
# -----------------------------
# 0) Shuffle x_raw
# -----------------------------
N <- length(x_raw)
x_shuf <- x_raw[sample.int(N)]
# Same preprocessing as the original: center by median and integrate
x_shuf_centered <- x_shuf - median(x_shuf, na.rm = TRUE)
y_int_shuf <- cumsum(x_shuf_centered)
# -----------------------------
# 1) Estimate lambda^2 across scales (SG only) for the shuffled series
# -----------------------------
lambda2_with_sg_shuf <- rep(NA_real_, length(scale_list))
for (k in seq_along(scale_list)) {
s <- scale_list[k]
# SG detrending (window length = s) + endpoint cropping
y_det <- sg_detrend_and_crop(y_int_shuf, sg_order = sg_order, sg_window = s)
# Increments and standardization
dy <- diff_by_scale(y_det, s)
dy_std <- standardize(dy)
# lambda^2 estimate (same moment order q as the original)
lambda2_with_sg_shuf[k] <- estimate_lambda2_from_q(dy_std, q_moment)
}
# Helper: fetch lambda^2 at an exact scale, if present in scale_list
get_lambda2_at_scale <- function(s, scale_list, lambda2_vec) {
idx <- which(scale_list == s)
if (length(idx) == 1) return(lambda2_vec[idx])
NA_real_
}
# -----------------------------
# 2) Plot: shuffled series + lambda^2 comparison + PDFs (shuffled)
# -----------------------------
layout(matrix(c(1, 1, 2,
3, 4, 5),
nrow = 2, byrow = TRUE))
par(las = 1, xaxs = "i", mar = c(5, 5, 2, 2), cex.axis = 1.2, cex.lab = 1.5)
# (1) Shuffled x[n] (use the same y-range as the original for fair comparison)
x_max <- max(abs(x_raw), na.rm = TRUE)
plot(0:(N - 1), x_shuf, type = "l", col = 2,
xlab = "n", ylab = "x_shuf[n]",
ylim = c(-x_max, x_max),
main = "Shuffled series x_shuf[n]")
# (2) lambda^2(s): original vs shuffled
ylim_top <- max(c(lambda2_with_sg, lambda2_with_sg_shuf), na.rm = TRUE)
if (!is.finite(ylim_top)) ylim_top <- 1
plot(scale_list, lambda2_with_sg,
col = 4, pch = 16, log = "x",
xlab = "Coarse-graining scale s (window length)",
ylab = expression(lambda^2),
main = sprintf("Scale dependence of non-Gaussianity (q = %.2f)", q_moment),
ylim = c(0, ylim_top))
points(scale_list, lambda2_with_sg_shuf, col = 2, pch = 16)
lines(scale_list, theory_lambda2(scale_list), lty = 2)
abline(h = 0, lty = 2)
legend("topright", bty = "n",
legend = c("Original", "Shuffled"),
col = c(4, 2), pch = 16)
# (3)-(5) PDFs at selected scales (shuffled series only)
col_kde <- rgb(1, 0.5, 0.5, alpha = 0.3) # light blue (KDE)
col_mod <- 2 # blue (model)
for (s in scale_show) {
y_det <- sg_detrend_and_crop(y_int_shuf, sg_order = sg_order, sg_window = s)
dy <- diff_by_scale(y_det, s)
dy_sg_std <- standardize(dy)
lam2_sg <- get_lambda2_at_scale(s, scale_list, lambda2_with_sg_shuf)
if (!is.finite(lam2_sg)) lam2_sg <- estimate_lambda2_from_q(dy_sg_std, q_moment)
plot_pdf_kernel_with_model_sg(
dy_sg_std = dy_sg_std,
lam2_sg = lam2_sg,
s = s,
xlim = c(-7, 7), ylim = c(1e-4, 1),
bw = "nrd0", adjust = 1, kernel = "gaussian",
col_kde=col_kde, col_mod=col_mod
)
}
まとめ
多重スケール確率分布解析については、まだ、よくわからないと感じている方も多いのではないでしょうか。今後は、この手法について繰り返し丁寧に説明していきたいと思います。本手法は、私自身が提案してきた解析方法ですので、その考え方や有効性を多くの人に知ってもらい、使いこなしてほしいと思います。
※ もし記事の中で「ここ違うよ」という点や気になるところがあれば、気軽に指摘していただけると助かります。質問や「このテーマも取り上げてほしい」といったリクエストも大歓迎です。必ず対応するとは約束できませんが、できるだけ今後の記事で扱いたいと思います。それと、下のはてなブログランキングはあまり信用できる指標ではなさそうですが (私のブログを読んでいる人は、実際とても少ないです)、押してもらえるとシンプルに励みになります。気が向いたときにポチッとしていただけたら嬉しいです。
