2023 年度科学研究費シンポジウム

統計科学・機械学習・情報数学の最前線⚫

科学研究費補助金 基盤研究(A) 20H00576「大規模複雑データの理論と方法論の革新的展開」研究代表者：青嶋 誠（筑波大学）によるシンポジウム⚫

日時 : 2024 年 1 月 26 日（金）〜1 月 27 日（土）⚫

場所 : 東北大学 青葉山キャンパス 情報科学研究科棟 2 階大講義室⚫

開催責任者 : 荒木 由布子 (東北大学）⚫

内容・目的 : 複雑で多様な形態，サイズ，性質およびドメインを持つデータが観測・測定されるようになり，データ科学に基づく研究は多様な分野で急速に進化しています。その成果は学術界だけではなく，例えば根拠に基づく意思決定の手法の１つとしても，広く社会に浸透し始めています。本シンポジウムでは，データ駆動型研究でも特に，統計科学・機械学習・情報数学における新たな方法論の提案や最近の研究動向の紹介，教育に関連する報告，注目されている手法の解説とともに，応用統計学，バイオ統計学，材料科学，生命科学，経済学，ファイナ
ンス，認知科学，社会科学，自然科学など，多岐にわたる科学分野における応用例の紹介や問題提起などを，広く対象としています。分野横断型のシンポジウ
ムとすることで，参加者が各分野の情報を共有し，幅広い視点からの意見交換により，問題解決や新たな研究の発展につながることを期待します。

講演プログラム
1 月 26 日（金）10:00 ~ 10:05 オープニング
10:05 ~ 10:35 大石 峰暉（東北大学）
Kick-one-out 法に基づいた回帰におけるクラスタリング
10:35 ~ 11:05 永井 勇（中京大学）
GMANOVA モデルにおける様々な仮定が不要な罰則無推定法とその解釈
11:05 ~ 11:35 松井 秀俊（滋賀大学）, 寺田 吉壱（大阪大学/理研 AIP）
可変ドメイン関数データに対する軌道の動的予測と気象データ分析への応用
11:35 ~ 13:00 昼休み
13:00 ~ 14:00 ポスターセッション− 若山 智哉（東京大学）, 菅澤 翔之助（慶應義塾大学), 小林 弦矢（明治大学）
Similarity-based Random Partition Distribution for Clustering Functional Data
− 竹下 和慶（大阪大学), 寺田 吉壱（大阪大学/理研 AIP）
深層学習を用いた function-on-scalar 回帰とその理論的性質− 岡崎 彰良（九州大学）, 川野 秀一（九州大学）
頑健クラスタリングに基づくマルチタスク学習回帰− Huang, Hongming（大阪大学）, 鈴木 讓（大阪大学）
An Efficient Procedure for Computing the Bayesian Network Structure Learning− Liu Lirui（大阪大学）, 鈴木 讓（大阪大学）
An Information Criterion Improving WBIC and SBIC
− 池上 孝則（放送大学）
データ量子統計の提案− 松野 大輔（千葉大学), 内藤 貫太（千葉大学)
Semiparametric Density Estimation with Bregman Divergence− 蔣 思澄（千葉大学）, 内藤 貫太（千葉大学）
一般化 Kendall’s tau 行列による主成分分析− 森本 直紀（東京理科大学）, 八木 文香（東京理科大学），瀬尾 隆（東京理科大学）
3-step 単調欠測データにおける T 2 型検定統計量を用いた平均ベクトルの検定− Yao Sai（東北大学）, 荒木 由布子（東北大学）
Modal Interval Regression for Visualizing Data-frequent Areas− Bat-Erdene, A.（Iwate University）, Kawasaki, S（Iwate University）
Characterization of Asymptotic Distribution in Martingale Convergence of theSupercritical Branching Processes with Poissonian Offsprings

14:00 ~ 14:10 休憩
14:10 ~ 15:00 特別講演 1 山口 類（愛知県がんセンター, 名古屋大学）
タンパク質結合能予測のためのマルチモーダル機械学習モデルの開発
15:00 ~ 15:30 竹内 努 (名古屋大学, 統計数理研究所), 本武陽一 (一橋大学), 西澤 淳(岐阜聖徳学園大学, 名古屋大学), 村上 広椰 (名古屋大学), 加納 龍生(名古屋大学), 馬 海霞 (名古屋大学), 内田 舜也 (名古屋大学), Suchetha COORAY (国立天文台科学研究部), 河野 海 (名古屋大学)Application of the Topological Data Analysis to Cosmology（位相的データ解析の宇宙論への応用）
15:30 ~ 15:45 休憩
15:45 ~ 16:15 吉羽 要直（東京都立大学）
代表的な接合関数の裾従属係数と裾次数16:15 ~ 16:45 江村 剛志（統計数理研究所）
コピュラ・グラフィック推定量―理論と応用―
16:45 ~ 17:35 特別講演 2 服部 聡（大阪大学）
医学研究における不完全データの感度解析
18:00 ~ 20:00 懇親会（東北大学青葉山キャンパス内 Buddy's Table）
1 月 27 日（土）9:30 ~ 10:00 石井 一夫（公立諏訪東京理科大学）, 金打 聖菜 (公立諏訪東京理科大学)
野菜の価格と気象との関連性について
10:00 ~ 10:30 本田 敏雄（一橋大学）
Sparse Quantile Regression via ℓ_0 penalty
10:30 ~ 10:40 休憩
10:40 ~ 11:30 特別講演 3 今泉 允聡（東京大学/理化学研究所）
インコンテキスト学習の統計学的解析
11:30 ~ 12:50 昼休み
12:50 ~ 13:20 徳田悟（九州大学）
対数周辺尤度とその陰的正則化のスケーリング挙動
13:20 ~ 13:50 川久保 友超（千葉大学）, 各務 和彦（名古屋市立大学）
所得不平等度の分解の推定
13:50 ~ 14:00 休憩14:00 ~ 14:30 中北 誠（理化学研究所 AIP）, 鳥谷部 智規（慶應義塾大学）,中妻 照雄（慶應義塾大学）
Bayesian Analysis of Stochastic Conditional Duration Models with Intraday andIntra-deferred Future Seasonalities in High-frequency Commodity Market
14:30 ~ 15:00 渡名喜 庸蔵（大阪大学）
微小ノイズをもつ空間 2 次元線形放物型確率偏微分方程式モデルの係数パラメータの推定
15:00 ~ 15:10 休憩
15:10 ~ 16:00 特別講演 4 松田 安昌（東北大学）
VAR Models for Surface Time Series
16:00 ~ 16:05 クロージング

Kick-one-out 法に基づいた回帰におけるクラスタリング
大石 峰暉
東北大学 データ駆動科学・AI 教育研究センター
n 次元目的変数ベクトル y, n × k 説明変数行列 Z に対する以下のモデルを考える.
y = Zγ + ε.

(1)

ここで, γ = (γ1 , . . . , γk )′ は k 次元回帰係数ベクトル, ε は平均ベクトル 0n , 共分散行列 σ 2 In を持つ n 次元
誤差ベクトルであり, rank(Z) = k ≤ n を仮定する. Kick-One-Out 法 (Zhao et al., 1986; Bai et al., 2018)に基づいた k 個の説明変数のクラスタリング法を提案する.
まずは通常の KOO 法について説明する. KOO 法は変数選択を行う手法で, フルモデル (すべての説明
変数を用いたモデル) を基準として, 説明変数ごとにその変数が必要がどうかを判定する. c0 をフルモデルの良さ, cj (j ∈ {1, . . . , k}) をフルモデルから第 j 説明変数を除いたモデルの良さとする. ここで, モデルの良さはモデル選択規準などで測ることができる. このとき, 第 j 説明変数は cj > c0 であれば必要, cj ≤ c0であれば不要と判断する. つまり, KOO 法で選ばれる説明変数は以下のように定義される.

{j ∈ {1, . . . , k} | cj > c0 } .
KOO 法では (k + 1) 回の計算で変数選択ができるため, 説明変数が多いときにも適用可能である. また,様々な設定の下で変数選択の一致性が理論的に保証されている (e.g., Bai et al., 2018; Oda & Yanagihara,2020).
次に, 提案手法について説明する. 先に述べた通り, KOO 法はフルモデルと変数を 1 つ除いたモデルを比較し, その変数が必要かどうかを判定した. 提案手法では, フルモデルと 1 組の変数を結合したモデルを比較し, その結合を受理するかどうかを判定することでクラスタリングを行う. ここでのフルモデルとは, どの変数も結合されていないモデルを意味する. フルモデルにおける γ の最小二乗推定量, および当てはめ値はそれぞれ以下で与えられる.

γ̂ = (Z ′ Z)−1 Z ′ y,ŷ = Z γ̂.

例えばフルモデルの良さを一般化 Cp 規準 (Atkinson, 1980) で測るとき, フルモデルの良さ c0 は以下で与
えられる.

c0 =

1
∥y − ŷ∥2 + αk,
s2

s2 =

∥y − ŷ∥2.
n−kここで, α はモデルの当てはまりに対する罰則の強さを表す正の値である. 次に結合の候補を考える. E を結合の候補となる添え字の組の集合とする. つまり, E は {1, . . . , k}2 の部分集合で, (j, ℓ) ∈ E ⇒ j < ℓ とする. 任意の (j, ℓ) ∈ E に対し, 第 j, ℓ 変数を結合したモデルとは, モデル (1) に γj = γℓ の制約を課したモデルを意味する. この制約の下では, モデル (1) は以下のように書き換えられる.

y = Zjℓ γjℓ + ε.
ここで, Zjℓ は Z から第 ℓ 列ベクトルを除き第 j 列に加えた n × (k − 1) 行列, γjℓ は γ から第 ℓ 成分を除いた (k − 1) 次元ベクトルである. このモデルに対して, γjℓ の最小二乗推定量, および当てはめ値はそれぞれ以下で与えられる.
′′γ̂jℓ = (ZjℓZjℓ )−1 Zjℓy,ŷjℓ = Zjℓ γ̂jℓ .

よって 1 組の変数 (j, ℓ) を結合したモデルの良さ cjℓ は, 一般化 Cp 規準を用いると以下で与えられる.

cjℓ =

1
∥y − ŷjℓ ∥2 + α(k − 1).
s2

以上のように定義された c0 , cjℓ に対し, 最適な結合を以下で定義する.

Ê = {(j, ℓ) ∈ E | cjℓ ≤ c0 } .
モデル選択規準はその値が小さいほど良いモデルであると判断できる. したがって, cjℓ < c0 であれば 2 つの変数の結合がモデルを改善することからその結合を受理し, cjℓ > c0 であればモデルが改悪されることか
らその結合を受理しない. このようなルールに従って変数の組ごとにそれらを結合するかどうかを判定し,その結果を統合することで最終的なクラスタリング結果を得ることができる.
本発表では, 特に真のクラスターの選択に関する一致性に着目し, 提案手法と Ohishi et al. (2021) による一般化 fused Lasso を用いた手法を数値的に比較した. 一般化 Cp 規準の α として以下を用いたとき, 提案
手法は真の結合の選択に関する一致性を有することが数値的に確認された.
√n−k(n − k) n − k − 3√α=· k 1/4 log n.
+ β, β =n−k−2(n − k − 2) n − k − 4
上記の α は Oda & Yanagihara (2020) で示された変数選択の一致性を有する値である.

参考文献
Atkinson, A. C. (1980). A note on the generalized information criterion for choice of a model. Biometrika,
67, 413–418.
Bai, Z. D., Fujikoshi, Y. & Hu, J. (2018). Strong consistency of the AIC, BIC, Cp and KOO methods inhigh-dimensional multivariate linear regression. Technical Report TR-No. 18-09. Hiroshima Statistical
Research Group. Hiroshima.
Oda, R. & Yanagihara, H. (2020). A fast and consistent variable selection method for high-dimensional
multivariate linear regression with a large number of explanatory variables. Electron. J. Stat., 14, 1386–
1412.
Ohishi, M., Fukui, K., Okamura, K., Itoh, Y. & Yanagihara, H. (2021). Coordinate optimization forgeneralized fused Lasso. Comm. Statist. Theory Methods, 50, 5955–5973.
Zhao, L. C., Krishnaiah, P. R. & Bai, Z. D. (1986). On detection of the number of signals in presence ofwhite noise. J. Multivariate Anal., 20, 1–25.

GMANOVA モデルにおける様々な仮定が不要な罰則無推定法とその解釈
中京大学 教養教育研究院 永井 勇
n 個の各個体に対して，時間と共に測定されたデータは経時測定データと呼ばれ，非常に多くの分野で分析されている．このような経時測定データにおいて，全ての個体で p 回測定しその測定時点が揃っている場合，Potthoﬀ and Roy (1964) により提案された次の一
般化多変量分散分析（GMANOVA）モデルでの分析がよくされている;

Y = 1n µ′ X ′ + AΞX ′ + E ，(1)

ここで，Y は各行が各個体の経時測定データからなる n × p 既知行列，1n は全ての要素が

1 の n 次元のベクトル，µ は q 次元未知ベクトル，X は個体内計画行列と呼ばれる測定時点から作られる p × q 既知行列 (詳細は後述)，A は個体間計画行列と呼ばれる測定時点に無関係な各個体の特徴を表す説明変数からなる n × k 既知行列で各列で中心化されているとし，Ξ は k × q 未知行列であり，E は E[E] = 0n 0′p で Var[vec(E)] = Σ ⊗ In の n × p 誤
差行列であり，0r は全ての要素が 0 の r 次元ベクトル，Σ は未知の p × p 正則行列である．

Y の経時的な変動 (経時変動) を推定する際には，個体内計画行列 X の i 行目として，i番目の測定時点 ti (i = 1, . . . , p; t1 < · · · < tp ) の関数からなる行列が用いられる．例えばX の i 行目として (tiq−1 , . . . .t0i ) を用いて未知の µ や Ξ を推定することで，全ての個体で共通の経時変動 や 各説明変数に対応する経時変動 の切片と各次数の係数がそれぞれ推定
でき，それらの結果として目的の経時変動を測定時点 ti の (q − 1) 次多項式で推定できる．
また，X の i 行目として ti の別の関数からなるものを用いて，µ や Ξ を推定することで，用いた関数の重み付き和の形で全体の経時変動を推定できる．
このモデル (1) において，µ や Ξ の推定は次のリスクを最小にすることで行われる;

tr{(Y − 1n µ′ X ′ − AΞX ′ )Σ−1 (Y − 1n µ′ X ′ − AΞX ′ )′ }．

(2)

このリスクを最小にする µ や Ξ を求める際に，従来は rank(A) = k かつ rank(X) = q を仮定している．これらの仮定は暗に k ≤ n や経時変動の推定に用いる ti の関数への制限をしている．一方で，これらの仮定が満たされない場合，以下のような手法がある:

 リスク (2) に罰則を加えたものを最小化する罰則付き推定法 (例えば Katayama &Imori (2014) 参照) を用いる
 Yanagihara ら (2023) などの情報量規準を用いて，A や X からいくつかの列を選択し，選択したものだけを推定に用いる

 Moore-Penrose 型一般逆行列 (例えば Srivastava (2002) A.6) を用いる
本講演では，永井 (2021) のアイデアを用い，rank(A) < k かつ rank(X) < q の下で，リスク (2) を最小にする推定量を提案する．さらに，本講演で提案する推定量の解釈などについても触れる．数値実験を通じた比較なども講演で触れる予定である．

引用文献:[1] Katayama, S. & Imori, S. (2014) Lasso penalized model selection criteria for highdimensional multivariate linear regression analysis. J. Multivariate Anal., 132, 138–
150.

[2] 永井 勇 (2021) 高次元小標本における多変量線形回帰モデルでの推定法，2021 年度統計関連学会連合大会 講演.
[3] Potthoﬀ, R. F. & Roy, S. N. (1964) A generalized multivariate analysis of variancemodel usefule specially for growth curve problems. Biometrika, 51, 313–326.
[4] Srivastava, M. S. (2002) Methods of Multivariate Statistics, Wiley-Interscience.
[5] Yanagihara, H., Nagai, I., Fukui, K. & Hijikawa, Y. (2023) Modiﬁed Cp criterion
in widely applicable models. Smart Innov. Syst. Tec., 352, Intelligent Decision
Technologies 2023: Proceedings of the 15th KES International Conference on Intelligent Decision Technologies KES-IDT-2023 (eds. I. Czarnowski, R. J. Howlett & L.
C. Jain), 173-182.

可変ドメイン関数データに対する軌道の動的予測と気象データ分析への応用
滋賀大学 データサイエンス学部
大阪大学 大学院基礎工学研究科 / 理化学研究所 AIP

松井 秀俊
寺田 吉壱

1. 概要
観測個体がそれぞれ経時的に測定された1000

データを関数化処理し，得られた関数データ集合を対象として分析を行う方法は関数960

Pressure

ざまな分野で応用されている (Kokoszka &980

データ解析とよばれ，医学や工学などさま

900

920

940

Reimherr, 2017)．関数データ解析では，観測される関数データの定義域は観測個体間で同じという状況が多い．これに対して本
研究では，定義域の終端となる時点が個体ごとに異なるデータを扱う．このようなデー
0
20
40
60
タを，ここでは可変ドメイン関数データとTimeよぶ．たとえば，日本付近を通過した一部の台風についての発生から消滅までの中心
図 1: 台風の中心気圧（出典：デジタル台風）．
気圧の推移は台風が消滅するまでの時間が台風の発生時点を 0 で揃えている．
台風毎に異なるため, 従来の関数データとして捉えることができない（図 1）．また,
一般的な傾向として，台風は消滅までの時間が長いものほど中心気圧の最小値は低く，短いものほど高い．可変ドメイン関数データは，その特性が定義域の長さに応じて異なることから，定義域の長さを考慮に入れた分析が行われている (Geller et al., 2012; Johns et al., 2019)．
本研究では，1 つの特徴についての可変ドメイン関数データが得られているとき, 新たな観測個体について経時測定データが途中まで観測されている状況で，その先の軌道および終端の時点を予測する方法について検討する．図 1 の例では，これまでに赤や灰の曲線に対応する台
風が観測され，さらに青の曲線に対応する台風が実線部分まで観測されたとき，その先の破線部分の軌道と，その終端時点を予測するという問題を考える．そのために，可変ドメイン関数
データに対して従来の関数データに対する動的予測の考え方 (Chiou, 2012) を適用する．提案
手法の有効性を，人工データおよび実データの分析を通して検証する．

2. 可変ドメイン関数データの動的予測
いま，始点を 0，終点を確率変数 M ∈ (0, 1] としたランダム区間 I = [0, M ] 上のランダム関
数 Y (t) (t ∈ I) を考える．終端 M = m が与えられた下で, Y (t) は区間 [0, m] 上の平均二乗連
続な確率過程となることを仮定する. さらに，X を共変量を表す確率変数ベクトル，その実現

値を x とする．このとき, Y (t) の平均関数および共分散関数を次で表す．

µ(t | m, x) = E[Y (t) | M = m, X = x],γ(t, s | m, x) = E [{Y (t) − µ(t | m, x)} {Y (s) − µ(s | m, x)} | M = m, X = x] .
あるランダム関数 Y の軌道が，ある時点 τ ∈ I について S = [0, τ ] の区間でのみ観測された状況を考える．S の区間での軌道を表すランダム関数を YS , 実際に観測された軌道を yS と記
す．このとき，τ まで観測されているランダム関数 Y の終端 M を，次の条件付き期待値で予測する．

c = E [M | M > τ, YS = yS , X = x] .
M
実際には，終端が m > τ なる観測データに対して，区間 S での関数データ y(t) および共変量x を説明変数，対応するデータの終端 m を目的変数として，ノンパラメトリック回帰などを用いる．そして，S 上でのみ観測されている Y の終端 M を予測する．その予測値を mb と記す.
次に，τ より先の区間 T = (τ, m]b における Y (t) の軌道 YT を，次で予測する．
hi
YbT = E YT M = m,b Y S = yS , X = x .
この条件付き期待値を計算するにあたり，M = m と X = x が与えられた下でランダム関数

Y (t) が平均関数 µ(t | m, x)，共分散関数 γ(t, s | m, x) のガウス過程に従うと仮定する．このと⊤ ⊤き，Y (t) を離散時点で評価し，区間 S, T で分割したベクトル y = (y ⊤S , y T ) は次の平均ベク
トル，分散共分散行列をもつ多変量正規分布に従うとする．
!
!
!!
yS
b x)
µS (m,
ΓSS (m,b x) ΓST (m,
b x)
∼N
,
.
yT
µT (m,b x)
ΓT S (m,b x) ΓT T (m,b x)
そして，YT (t) を離散時点で評価したベクトル y T を，次の事後平均によって予測する．

b T = µT (m,b x) .
b x) y 0,S − µS (m,yb x) + ΓT S (m,b x)Γ−1
SS (m,
平均ベクトルや分散共分散行列は，Johns et al. (2019) の方法を利用して推定する．
提案した手法を台風のデータ分析に適用し，途中まで観測された台風の中心気圧の推移から，その先の軌道と，台風が消滅する時点を予測する．結果については，当日報告する．

参考文献
Chiou, J. M. (2012). Dynamical functional prediction and classiﬁcation, with application to traﬃc ﬂow prediction. Ann. Appl. Statist. 6, 1588–1614.
Gellar, J., Colantuoni, E., Needham, D. M., and Crainiceanu, C. (2014). Variabledomain functional regression for modeling ICU data. J. Amer. Stat. Assoc. 109,
1425–1439.
Johns, J. T., Crainiceanu, C., Zipunnikov, V, and Gellar, J. (2019). Variable-domainfunctional principal component analysis. J. Comput. Graph. Stat. 28, 993–1006.
Kokoszka, P. and Reimherr, M. (2017). Introduction to functional data analysis.
CRC Press, Boca Raton.

Similarity-based Random Partition Distribution for
Clustering Functional Data

1

若山智哉

菅澤翔之助

小林弦矢

東京大学経済学研究科

慶應義塾大学経済学部

明治大学商学部

はじめに
過去数十年にわたるモバイルデバイスの急速な普及により、人口データは現代社会の理解や政策立案において重要な役割を担っている。具体的には、都市計画、交通計画、医療サービスの提供、広範囲な政策立案など、多岐にわたる分野に影響を及ぼしている [1, 3, 5]。この発表では、人口データを深く理解するための確率分割を提案し、その応用について説明した。

2

データと背景クラスタリングは人口データから洞察を引き出すために有効である。クラスタ内の共通点を特定し、地域特性を明らかにすることで、ライドシェアサービスや街頭広告などの産業と、都市工学
や人文科学などの研究領域の両方に貢献することができる。特に、モデルベースのクラスタリン
グは、ノンパラメトリックベイズ手法と相性が良い。ディリクレ過程 (DP) に基づくクラスタリ
ングは、自動的にクラスタ数を決定し、空間構造を組み込むことも容易である。[2] の研究では、

Pitman–Yor 過程の中で既存のクラスターを新しい項目に割り当てる際に地理的な近接性を統合しており、広く応用されている。特に、人口データ解析でも隣接情報が重要であり、図 1 に示すように、隣接する地区間には高い相関が見られる。
にもかかわらず、我々のデータセットへの既存のクラスタリング手法の適用には、その多変量性と時間的情報の 2 つの特性が障壁となっている。多変量クラスタリング問題に取り組む際、クラスタ数を固定しない手法を単純に適用することはリスクを伴う。高次元データはその複雑性から、容

60000
40000
20000

図 1 2019 年 1 月 29 日 14 時の東京都の中心地区における人口マップ

易に分類可能であり、これは教師あり分類タスクでは役立つ特徴である [4]。一方で、この高い分類可能性は、教師なしクラスタリングではクラスタ数が過剰になる可能性がある。特に、この性質は DP が小さいクラスターサイズを好む傾向を推進し、余計なクラスターを作ることになり、デー
タの理解を妨げることになる。また、人口データの解析では、時間的な構造も重要である。地域に
よっては平日と休日では人口の増減が大きく異なるため、グループ構造は日毎に変化ておりし、典型的な平日の傾向だけに注目すると重要な洞察を見逃してしまう可能性がある。

3

方針
本研究では一般化ディリクレ過程 (GDP) を基にした新たな確率分割手法、SGDP を提案する。

この手法は、観測された時間毎の人口フローを関数（確率過程の実現値）としてクラスタリングすることを目的としている。SGDP では、一般化パラメトリゼーションを用いることで過剰なクラスタの形成リスクを低減し、同時に地理的隣接性を考慮しており、この方法の空間的な性質についても理論的な裏付けがある。さらに、観測値とクラスタ平均との間の差異をガウス過程でモデル化する。これは、ある時刻での観測された人口がクラスタ平均から逸脱する場合、その影響が連続する時刻にも及ぶという考え方に基づいており、単なる独立したノイズとして処理するよりも、過剰なクラスタ形成を抑制する効果がある。また、この方法は時間的な変動にも柔軟に対応できるため、各地域の時間によって変わるパターンをより詳細に捉え、理解することが可能になる。

4

結果

データ解析の結果や手法の有効性はポスターセッションで議論した.

参考文献[1] Amir Ahmadi-Javid, Oded Berman, and Pooya Hoseinpour. Location and capacity planning of facilities with general service-time distributions using conic optimization. arXiv
preprint arXiv:1809.00080, 2018.
[2] David B Dahl, Ryan Day, and Jerry W Tsai. Random partition distribution indexed bypairwise information. Journal of the American Statistical Association, 112(518):721–732,
2017.
[3] Antonio Páez and Darren M Scott. Spatial statistics for urban analysis: A review oftechniques with examples. GeoJournal, 61:53–67, 2004.
[4] Tomoya Wakayama and Masaaki Imaizumi. Fast convergence on perfect classiﬁcation forfunctional data. Statistica Sinica, 34(4), 2024.
[5] Mingshu Wang and Lan Mu. Spatial disparities of uber accessibility: An exploratory
analysis in atlanta, usa. Computers, Environment and Urban Systems, 67:169–175, 2018.

深層学習を用いた function-on-scalar 回帰とその理論的性質
竹下 和慶 1 , 寺田 吉壱 1,2
1 大阪大学大学院基礎工学研究科, 2 理化学研究所革新知能統合研究センター

1. はじめにある領域や区間上で連続的に観測されたデータをランダムな関数の実現値として捉えたデータ
解析は関数データ解析 (Functional Data Analysis; FDA) と呼ばれている. 本研究では, 入力デー
タ（予測変数）が従来の多変量データで, 出力（応答変数）が関数データであるような FunctionOn-Scalar (FOS) 回帰と呼ばれる関数回帰分析を考える.
多くの既存研究は, 単純で扱いやすい線形な FOS 回帰に焦点を当てており, 非線形な FOS 回帰についての研究結果は少ない. 例えば, Scheipl et al. (2015) は予測変数に関して非線形な関数型
加法モデル (Functional Addive Mixed Model; FAMM) を提案している. FAMM では, 応答変数
が予測変数の一部を引数とした非線形関数の和として表現できると仮定している. この手法は, 従来の線形 FOS 回帰よりも高い柔軟性を有するモデルである. しかし, 仮定した加法モデルがデータに適合しない場合には推定精度が低下してしまうという問題がある.
Luo et al. (2023) は Scheipl et al. (2015) の問題点を解消するため, より柔軟なモデリングが可能となるニューラルネットワークに基づく手法を提案した. また, 従来の万能近似定理を拡張し,予測変数に関して連続な関数に対してこのモデルの万能近似性を示している. この手法の利点は、従来の FAMM が基底関数を用いて真の関数を直接的に推定するのに対し、この手法では予測変数の次元数に基づいた複数の二変数関数を推定するのみで良いという点である. このことによって、FAMM よりも計算コストが低く、予測変数の次元が比較的大きい場合でも適用可能である.
しかし, この方法には, 2 つの問題点がある. １つ目は, 真の関数が予測変数に関する連続性を仮定している点である. 2 つ目は, 提案された推定関数自身が積分形で記述されているため実用上厳密にモデル化することが不可能な点である. そのため, 実際に手法を適用する際には, 時間と積分に関して離散化し基底関数を用いて推定している. Luo et al. (2023) の Theorem 3 では, 推定量の収束レートが導出されているが, 証明に重大なミスがあり, 正しい収束レートを導出できていな
い. 本発表では, 深層ニューラルネットワークを用いたより自然な FOS 回帰分析を提案し, その推定量の収束レートを導出する.
2. 提案手法と主結果
本研究では, Luo et al. (2023) の問題点を解決した深層ニューラルネットワークを用いたより直接的な推定量を提案する. また, 真の回帰関数 f ◦ が非等方的 Besov 空間 (Anisotropic Besov space)に含まれるという仮定の下で, 推定量の理論的性質を導出する. 本研究で扱う非等方的 Besov 空間
は, 方向によって滑らかさが異なるような複雑な関数を含む関数空間であり, Hölder 空間などを特別な場合として含んでいる. 従って, Luo et al. (2023) よりも柔軟な設定で議論することが可能である. また, 非等方的 Besov 空間において深層ニューラルネットワークは適応的な近似能力を有し, 通常のノンパラメトリック回帰問題の設定において, 深層ニューラルネットワークを用いた最小二乗型推定量は, minimax 最適性をもつことが知られている (Suzuki and Nitanda, 2021).
予測変数の空間は d 次元とする. このとき, 入力層のニューロンが d + 1 個, 層の個数が L 個, 各層のニューロンが W 個, 出力層のニューロンが 1 個の順伝搬深層 ReLU ニューラルネットワーク

のクラス Φ(L, W, S, B) を用いて f ◦ を推定することを考える.
nΦ(L, W, S, B) := (W (L) η(·) + b(L) ) ◦ · · · ◦ (W (1) x + b(1) )|W (L) ∈ R1×W , b(L) ∈ R,

W (1) ∈ RW ×(d+1) , b(1) ∈ RW , W (l) ∈ RW ×W , b(l) ∈ RW (1 < l < L),
)
L
X
(l)
(∥W (l) ∥0 + ∥b(l) ∥0 ) ≤ S, max ∥W∞∥ ∨ ∥b(l) ∥∞ ≤ Bll=1

ここで, ∥ · ∥0 は非ゼロ成分の個数を表し, ∥ · ∥∞ は各成分の絶対値の最大値とする. また, S はス
パース定数, B はノルム定数である.
予測変数と応答変数の訓練データ Dn = {(Xl , Yl )}nl=1 が与えられたもとで、f ◦ に対する推定量fˆn ∈ Φ(L, W, S, B) を次の目的関数の最小化によって定める.

fˆn =

n

Z 1

l=1

0

1X
argminf ∈Φ(L,W,S,B) n{Yl (t) − f¯(Xl , t)}2 dt

ここで、f は f のクリッピングであり, ある定数 F > 0 に対して, f = min{max{f, −F }, F } と定められる.
ニューラルネットワークの大きさを適切に定めることで, 従来の非線形回帰問題 (Suzuki and
Nitanda, 2021) と同様に, 次の推定誤差レートが得られる.

Z 1
ED n
0

fˆn (X, t) − f ◦ (X, t)



2
L2 (P

X)

dt ≲ n

−

2β̃
2β̃+1

(log n)3

ここで, ∥ · ∥L2 (PX ) は予測変数の訓練データが従う確率密度関数で重みづけた L2 ノルム, EDn [·] は訓練データ Dn の分布に対する期待値を表す. また, β̃ は非等方的 Besov 空間の各軸方向における−

2β̃

滑らかさパラメータの調和平均である. この結果により、収束レートは log 項を無視すると n 2β̃+1であり, 滑らかさのパラメータ β̃ が大きければ収束レートが早くなることがわかる. さらに, この収束レートは直接的に次元 d に依存していないため, 次元の呪いを緩和できる可能性が示唆される. 具体的には, 真の関数が予測変数の空間で複雑であったとしても, 時間方向に対して十分な滑らかさを有する場合には, 良い収束レートが期待できる. 提案手法とその理論的性質の詳細や数値
実験は当日報告する.
参考文献

Scheipl, F., Staicu, A.-M. and Greven, S. (2015). Functional additive mixed models. Journal ofComputational and Graphical Statistics, 24(2), 477–501.
Luo, Ruiyan and Qi, Xin (2023). Nonlinear function-on-scalar regression via functional universal
approximation. Biometrics, 79(4), 3319–3331
Taiji Suzuki and Atsushi Nitanda (2021). Deep learning is adaptive to intrinsic dimensionalityof model smoothness in anisotropic Besov space. Advances in Neural Information Processing
Systems 34, 3609–3621

頑健クラスタリングに基づくマルチタスク学習回帰
九州大学 大学院数理学府 岡崎 彰良
九州大学 大学院数理学研究院 川野 秀一

1.

はじめに

マルチタスク学習とは，タスクと呼ばれる複数の関連していると考えられるデータ集合が得られている際に，それらのタスクそれぞれに異なるモデルを設定し，モデル間に共通した構造の推定を行うことでモデル全体の推定精度を改善する統計手法である (Caruana, 1997)．マルチタスク学習を実行する際，データによってはタスク間に共通する構造が単一ではないと考えられる場合が存在する (Xue et al., 2007)．このような状況におけるマルチタスク学習は，共通構造を共有するタスク集合の分類と共通成分の推定を同時に行うため，モデルのパラメータに対するクラスタリングにより実現される．また，マルチタスク学習では，用いられるタスクす
べてが何らかの関連性を有していることを仮定したもとで実行されるが，実際にはどのタスクとも関連性を有さない外れ値タスクが含まれている状況が考えられる (Gong et al., 2012)．
Okazaki and Kawano (2023) では，各タスクに設定された回帰モデルに関する損失関数に対し，凸クラス
タリングに基づく正則化項を加えることで，凸最適化の枠組みにおいてクラスタリングによるマルチタスク学
習を提案した．しかし，この手法ではすべてのタスクがいずれかのクラスタに属するという仮定のもと推定が行われるため，外れ値タスクがタスク集合に含まれている場合にクラスタリングの精度が悪化する問題が存在する．
そこで本報告では，凸クラスタリングによるマルチタスク学習 (Okazaki and Kawano, 2023) において，外れ値タスクの影響を軽減すると同時に外れ値タスクの検出を行う手法を提案する．各タスクに対して外れ値を表すパラメータを導入し，グループスパース正則化 (例えば Yuan and Lin (2006)) によりそのパラメータの選択を行う．提案手法は，タスクのクラスタリングを担う凸クラスタリングにおける損失関数を，外れ値に対して頑健な多変量損失関数 (例えば Hampel et al. (1986)) で置き換えた手法として解釈でき，タスクが構成するクラスタ推定において外れ値タスクの影響に対して頑健であることが期待できる．提案手法への有効性を，モンテカルロ・シミュレーションおよび実データへの適用を通して検証する．

2.

凸クラスタリングに基づくマルチタスク学習回帰モデル

いま，T 個の関連すると考えられるタスクが与えられているとし，m (m = 1, . . . , T ) 番目のタスクに対し，nm 個の説明変数に関する p 次元データ {xmi ; i = 1, . . . , nm } および nm 個の目的変数に関するデータ{ymi ; i = 1, . . . , nm } の組 {(ymi , xmi ); i = 1, . . . , nm } が独立に得られているとする．また，m 番目のタスクに関する目的変数ベクトルを ym = (ym1 , . . . , ymnm )⊤ ∈ Rnm ，計画行列を Xm = (xm1 , . . . , xmnm )⊤ ∈ Rnm ×p
とおく．これに対し，T 個のタスクの組を {(ym , Xm ); m = 1, . . . , T } とする．
次に，T 個のタスクに対して，目的変数と説明変数の関係を次の T 個の線形回帰モデル

ym = Xm w m + ϵm ,m = 1, . . . , T(1)により捉える．ここで，wm = (wm1 , . . . , wmp )⊤ は m 番目のタスクに対応する回帰係数ベクトルであり，ϵmは各要素がそれぞれ独立に N (0, σ 2 ) に従う観測誤差ベクトルである．(1) 式に対して，Okazaki and Kawano(2023) では，次の最小化問題


T
∑


T

∑
∑
1
λ
1
∥wm − um ∥22 + λ2
min p∥ym − Xm wm ∥22 +rm,l ∥um − ul ∥2
wm ,um ∈R 
2nm
2 m=1
m=1
m=1,...,T(2)(m,l)∈Eを提案した．ここで，um ∈ Rp は m 番目のタスクが属するクラスタの重心に関するパラメータ，E はタスクの組の集合，λ1 , λ2 は正の値を取る正則化パラメータ，rm,l は m 番目と l 番目のタスク間の関係を表す係数であり 0 以上の値を取る．目的関数の第二項目および第三項目は凸クラスタリング (Hocking et al., 2011) に基づく正則化項であり，wm に関してクラスタリングを実行する．同じクラスタに属するタスクは um が等し

い値として推定される．すなわち，um は各タスクが属するクラスタにおける回帰係数ベクトルの共通成分を表す．

3.

外れ値タスクに頑健な凸クラスタリング回帰モデル(1) 式に対して，次の最小化問題
min

wm ,um ,om ∈Rp ,m=1,...,T

{ T
∑

T
1
λ1 ∑
∥ym − Xm wm ∥22 +∥wm − um − om ∥22
2n
2
m
m=1
m=1

+λ2

∑

rm,l ∥um − ul ∥2 + λ3

T
∑
m=1

(m,l)∈E
P (om , η)(3)

を考える．ここで，om ∈ Rp は m 番目の外れ値タスクに関するパラメータ，λ3 は正の値を取る正則化パラメータである．第四項目は om に関して変数選択を行う正則化項であり，η は縮小される値を調整するパラメータである．第四項目により，om ̸= 0 として推定された場合，m 番目のタスクは他タスクと共通性を有さない外れ値タスクとして解釈される．第四項目にグループ lasso 正則化 (Yuan and Lin, 2006) を用いた場合，最小化問題 (3) は次の最小化問題と等価となる．


T
∑


T

∑
∑
1
min p
∥ym − Xm wm ∥22 + λ1
hλ3 (wm − um ) + λ2
rm,l ∥um − ul ∥2 .

wm ,um ∈R , 
2nm
m=1
m=1

(4)

(m,l)∈E

m=1,...,T

ここで，hλ3 (·) は{
hλ3 (a) =12
2 ∥a∥2λ2

λ3 ∥a∥2 − 23∥a∥2 ≤ λ3 ,∥a∥2 > λ3 ,(5)で定義される多変量フーバー損失関数 (Hampel et al., 1986) である．すなわち，提案手法は回帰係数に対してクラスタリングを実行する第二項目および第三項目を頑健化した手法として考えることができる．
最小化問題に含まれるパラメータの推定は，最小化問題 (3) に対する交互方向乗数法により実行される．数値実験の結果については当日紹介する．

参考文献
Caruana, R. (1997). Multitask learning. Machine Learning, 28, 41–75.
Gong, P., Ye, J., and Zhang, C. (2012). Robust multi-task feature learning. KDD ’12: Proceedings of the
18th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, 895–903.
Hampel, F. R., Ronchetti, E. M., Rousseeuw, P. J., and Stahel, W. A. (1986). Robust Statistics: The
Approach Based on Inﬂuence Functions. John Wiley & Sons.
Hocking, T. D., Joulin, A., Bach, F., and Vert, J. P. (2011). Clusterpath an algorithm for clustering
using convex fusion penalties. In Proceedings of the 28th International Conference on Machine Learning,
745–752.
Okazaki, A. and Kawano, S.
(2023).
Multi-task learning regression via convex clustering.
arXiv:2304.13342.
Xue, Y., Liao, X., Carin, L., and Krishnapuram, B. (2007). Multi-task learning for classiﬁcation withdirichlet process priors. Journal of Machine Learning Research, 8(2), 35–63.
Yuan, M. and Lin, Y. (2006). Model selection and estimation in regression with grouped variables.
Journal of the Royal Statistical Society. Series B (Statistical Methodology), 68(1), 49–67.

An Efficient Procedure for Computing Bayesian Network
Structure Learning
Hongming Huang1 and Joe Suzuki1
1 Osaka University, Graduate School of Engineering Science, Department of Systems Innovation

1

Introduction

In recent years, there has been active progress in the Bayesian Network Structure Learning(BNSL). Finding the optimal Bayesian network structure is an NP-hard problem[2], so the proposed
techniques are only applicable to small networks.
Tomi Silander and Petri Myllymäki proposed a straightforward score-based method for finding
the optimal Bayesian network[1]. This state-of-the-art method continuously compares the scoresof various structures to select the most favorable outcome. However, as shown on the left side of
Figure 1, it needs to calculate and record both the local scores Q and the optimal local solutions
R, which requires two passes over all possible local structures. Then, based on these information,a calculation for the globally optimal structure is performed, making it impossible to eliminate
redundant data during the computation process.
In this paper, we introduce an efficient method for layered computation to reduce the spacecomplexity. Differing from existing methods, it only requires a single pass through the possible
local structures as shown on the right side of Figure 1 and has the capability to eliminate redundant information of the completed levels during the computational process. For example, when
computing the topmost red layer, it only requires Q and R from the immediately lower yellow
layer, and the information from the green layer and any layers below can be completely discarded
as shown in Figure 2. This method retains only the crucial data necessary for calculations at
higher structural levels, eliminating the need to hold onto redundant data.

Figure 1: The left side of the dividing line represents the existing method where Q and R needto traverse possible local structures separately. The right side of the dividing line represents the
proposed method where Q and R can be updated together, requiring only a single traversal.

1

Y, Z, W
Y, Z W, Y Z, W
YZW

Figure 2: An example for proposed method. The red nodes represent the cases to be computed,
the yellow ones indicate the data information required to compute the yellow nodes, and the green
ones denote the nodes and data information that can be deleted.

2

Background2.1

Marginal Likelihood

In the realm of statistical analysis, the precision of parameter estimation is paramount. Traditional methods often fall short when parameters are unknown, leading to uncertainty in the results.
To overcome this, we leverage Bayesian inference, specifically the marginal likelihood Q(X), which
considers parameter uncertainty and integrates prior knowledge effectively.
We consider a Bayesian statistical model characterized by a set of realized values X for arandom variable X, and a parameter set Θ ⊆ Rd . Given a sample size n and i.i.d. observations
x1 , . . . , xn ∈ X , the model is described by a family of probability distributions {p(⋅∣θ)}θ∈Θ and aprior distribution φ over Θ. The marginal likelihood is defined as:n

Q(X) ∶= ∫ φ(θ) ∏ p(xi ∣θ)dθ,Θ

i=1

which synthesizes the likelihood across all parameter values, weighted by the prior.

2.2

Bayesian Network Structure Learning

Suppose that given n tuples of examples
X (1) = xi,1 , X (2) = xi,2 , . . . , X (p) = xi,p ,where i = 1, 2, . . . , n, with respect to p variables X (1) , X (2) , . . . , X (p) .
A fundamental step in structure learning is finding the optimal parent sets for each node
within the network. The joint probability of the variables X (1) , . . . , X (p) in a Bayesian network is
factorized according to the network structure as follows:p

P (X (1) , . . . , X (p) ) = P (X1 )P (X2 ∣π2 )⋯P (Xi ∣πi )⋯P (Xp ∣πp ) = ∏ P (Xi ∣ πi ),i=1where Xi denotes the vertex at i-th position in the structure and πi denotes the parent set of thevertex Xi .
For a Directed Acyclic Graph, a key fact is that due to the absence of cyclic structures, there
must exist a universally bottom-most node. This node is defined as the network’s sink node.
2

Assuming S is a non-empty set of nodes, it contains a bottom-most node X. The optimal
parent set for X lies within S/{X}, and identifying a new sink node and its optimal parent set
continues recursively within this subset until completion. Set V = {X (1) , . . . , X (p) }, with themarginal likelihood Q(X), the purpose is to find a BN structure with the maximum posterior
probability:p
Q(Xi , πi ), πi ⊆ V /{Xp , ⋯, Xi }.
R(X (1) , . . . , X (p) ) = max ∏πi
Q(πi )
i=1

3

Existing method

For the state-of-the-art method proposed by Silander and Myllymäki, they firstly compute
Q(S) for all subsets S ⊆ {X (1) , . . . , X (p) }. Then, for each subset S, they will compute the optimal
parent set for each node in it:
Q(X, T )Q(X, S)
Q(X, π(X, S/Y ))Q(X, π(X, S))= max= max {
, max
} ,
T ⊆S
Q(π(X, S))Q(T )Q(S) Y ∈S Q(π(X, S/Y ))where π(X, S) is the parent set of X ∈/ S in S ⊆ {1, . . . , p}. Then, they constructed an orderedgraph to store the parent sets {π(X, S)}S for each X ∈ {1, . . . , p}.
In the second stage, they combine the results obtained from above and the process for finding
the sink node:⎧
Q(X, π(X, S/X))⎪⎪⎪ max R(S/X), ∣S∣ ≥ 1
R(S) = ⎨ X∈S,
Q(π(X, S/X))⎪⎪ 1,⎪∣S∣ = 0⎩which recursively obtain the optimal parent set and the global optimal structure. To this end,
they constructed an additional ordered graph to store and recursively obtain the {R(S)} values.

4

Proposed Method

The proposed procedure seeks the same score more efficiently, constructing an ordered graphjust once. We claim our method can dramatically save the memory as shown in Table 1.
For every subset S where ∣S∣ ≤ k, we assume that Q(S), R(S), and π(X, S/X) for each X in
S have already been computed, the process is then repeated for each subset S where ∣S∣ = k + 1.
The following is the description of the algorithm steps:
1. We firstly compute Q(S) for the calculation in Step 2.
2. For each X ∈ S, we define π(X, S/X) by S/X and π(X, S/{X, Y }), Y ∈ S/X, when
Q(S)
Q(S/X)and

Q(X, π(X, S/{X, Y }))Q(π(X, S/{X, Y }))are maximized, respectively. In this step, we will update the optimal parent set for eachnode in S. Then we store the information and value for Step 3 and upper level iterations.
Q(X, π(X, S/X))
to obtain R(S).
Q(π(X, S/X))
Through this way, the optimal global structure is always obtained for each subset S.

3. In Step 3, we will find the X ∈ S that maximizes R(S/X)

3

Computation
memorySilander-MyllymakiO(p2 2p )O(p2p )

Proposed method
O(p2 2p )√O(2p / p)

Table 1: Comparison of computation and memory requirements

5

Summary

We plan to conduct experiments in the future to verify that the spatial complexity of the proposed method is lower than that of the compared methods in practical applications. Furthermore,
although we believe that the time complexities of both methods are similar, it has been suggested
that the Branch & Bound technique was employed to reduce computational time in the comparedmethod[3]. Therefore, our future work will also include the use of optimization techniques todecrease the time complexity of our proposed method.

References[1] Silander, T., & Myllymaki, P. (2012). A simple approach for finding the globally optimal
Bayesian network structure. arXiv preprint arXiv:1206.6875.
[2] Chickering, M., Heckerman, D., & Meek, C. (2004). Large-sample learning of Bayesian networks is NP-hard. Journal of Machine Learning Research, 5, 1287-1330.
[3] Suzuki, J., & Kawahara, J. (2017). Branch and Bound for Regular Bayesian Network Structure
Learing. In UAI.
[4] Kitson, N. K., Constantinou, A. C., Guo, Z., Liu, Y., & Chobtham, K. (2023). A survey ofBayesian Network structure learning. Artificial Intelligence Review, 1-94.
[5] Rantanen, K., Hyttinen, A., & Järvisalo, M. (2017). Learning chordal Markov networks via
branch and bound. Advances in neural information processing systems, 30.
[6] Colombo, D., & Maathuis, M. H. (2014). Order-independent constraint-based causal structurelearning. J. Mach. Learn. Res., 15(1), 3741-3782.
[7] Spirtes, P., Glymour, C. N., & Scheines, R. (2000). Causation, prediction, and search. MITpress.

4

An Information Criterion improving WBIC and SBIC
LIU LIRUIandJoe Suzuki
Mathematical Science Department, Graduate School of Engineering Sciences, Osaka University
We introduce a novel Information Criterion (IC) aimed at enhancing the functionality of WBIC and sBIC,
designed to be eﬀective without the regularity constraints and yet remain stable. ICs are essential in model
selection, balancing model ﬁt and complexity. Many existing criteria, like the singular Bayesian information criterion(sBIC), presuppose regularity in the true distribution relative to the statistical model. sBIC, despite its theoreticaladvantages in capturing large sample behavior (Drton et al., 2017), faces limitations in broader applications due
to its dependency on maximum likelihood estimates and regularity assumptions (Suzuki, 2023). This led to thedevelopment of criteria such as the Watanabe Bayesian information criterion (WBIC) (Watanabe, 2013), whichfunctions in non-regular scenarios but struggles with large sample sizes and the redundant estimation of known
learning coeﬃcients. Addressing these gaps, our proposed Information Criterion enhances WBIC and sBIC’s utility,oﬀering a ﬂexible and robust approach to model selection without the regularity constraints.
∏n
We denote by θ̂n the maximum likelihood estimate, i.e., the θ that maximizes the likelihood i=1 p (xi | θ)
obtained from x1 , . . . , xn ∈ X , and assume that the learning coeﬃcient λ has been obtained somehow. The Singular
Bayesian Information Criterion(Drton et al., 2017) is deﬁned by
SBIC :=n∑()− log p xi | θ̂n + λ log n.

i=1Assume that p, q are regular. When Θ is compact, the estimator θ̂n converges in probability to θ∗ as n → ∞(Suzuki, 2023). Note that the performance of sBIC is guaranteed only in the regular case, and sBIC cannot beapplied to general cases.
We introduce a parameter β > 0, called the inverse temperature, and assume that it decreases with n asβ = β0 / log n for some constant β0 > 0. Given x1 , . . . , xn ∈ X , we deﬁne the posterior mean of a function
f : Θ → R with respect to β by
∫
n
∏
f (θ)
p(xi |θ)β φ(θ)dθEβ [f (θ)] :=

Θ

∫ ∏
n

i=1

(1)
′ β

′

p(xi |θ ) φ(θ )dθ

′

Θ i=1

The Watanabe Bayesian Information Criterion(Watanabe, 2013) is deﬁned as[ n]∑
W BICn := Eβ− log p (xi | θ) .
i=1

From (Watanabe, 2013), WBIC behaves as√n∑

λ log nW BICn =+ Un− log p (xi | θ∗ ) +β0
i=1

λ log n+ Op (1),2β0(2)

where U1 , U2 , . . . are random variables that converge in law to a normal distribution with mean zero.
We deﬁne the predictive distribution and the empirical loss(Watanabe, 2010) given x1 , . . . , xn ∈ X as follows:r (x | x1 , . . . , xn ) = E1 [p(x|θ)],and1∑{− log r (xi | x1 , . . . , xn )} ,
n i=1
n

Tn :=

x ∈ X,

1

respectively. We propose a novel information criterion given by
newIC := nTn + λ log n.

(3)

In comparison with WBIC and sBIC, we argue the merits of newIC. First, unlike sBIC that relies on maximum
likelihood estimates, newIC employs empirical loss. Assuming β0 = 1, equation (2) becomes
W BICn =n∑− log p (xi | θ∗ ) + λ log n + Op (√log n),

i=1

where, in (2), the third term divided by√log n is OP (1). As per Suzuki (2023), Tn can be expressed as
( )
(
]) 1
1
1
1
1 [√
{− log p (xi | θ∗ )} +tξn (u) | x1 , . . . , xn − EX [V(X)] + oP
Tn =
λ− E
,
n i=1
n
2
2
n
n
∑

where E[·] and V[·] are the posterior expectationand variance] in equation (1) with β = 1 respectively. Noting that λ[√tξn (u) | x1 , . . . , xn are Op (1) due to convergence to their expected values
is constant and that nEX [V(X)] and E
by the law of large numbers, and substituting them into the deﬁnition (3) of newIC, we obtain
newIC =n∑{− log p (xi | θ∗ )} + λ log n + OP (1).

(4)

i=1

If λ is not known, it can be estimated by λ̂, obtained as√
W BIC1 − W BIC2λ̂ == λ + OP (1/ log n),1/β1 − 1/β2

(5)where this formula is derived from equation (2). Substituting λ with λ̂, newIC can be approximated as
newIC ≈n∑{− log p (xi | θ∗ )} + λ̂ log n + Op (√log n).

(6)i=1

Comparing equations (2), (4), and (6), we provide the following statements:(I) If λ is not known, newIC and WBIC perform almost identically as n → ∞.
(II) If λ is known, newIC outperforms WBIC.
We applied newIC and WBIC to reduced-rank regression and Gaussian mixture models. Speciﬁcally, in former
experiments using λ̂ for newIC calculations, WBIC and newIC showed similar performance, with newIC consistently
larger by a constant. For Gaussian mixture models, employing λ for newIC. Note that newIC is independent of β,
it does not change with diﬀerent β0 and remained accurate, while WBIC was incorrect for β0 = 10.
Our paper introduces the innovative newIC for model selection. Distinguishing itself from sBIC, newIC operates
eﬀectively beyond the conﬁnes of regularity, making it applicable in a wider range of scenarios. When contrasted
with WBIC, newIC demonstrates enhanced performance, especially in cases where the learning coeﬃcient is known.
This advancement solidiﬁes newIC’s position as a versatile and reliable tool for model selection in diverse statisticalcontexts.

References[1] Suzuki, J. (2023). WAIC and WBIC with R Stan: 100 Exercises for Building Logic. Springer Nature.
[2] Drton, M., Plummer, M. (2017). A Bayesian information criterion for singular models. Journal of the Royal
Statistical Society Series B: Statistical Methodology, 79(2), 323-380.
[3] Watanabe, S. (2013). A widely applicable Bayesian information criterion. The Journal of Machine Learning
Research, 14(1), 867-897.
[4] Watanabe, S. (2009). Algebraic geometry and statistical learning theory (Vol. 25). Cambridge university press.
[5] Watanabe, S., Opper, M. (2010). Asymptotic equivalence of Bayes cross validation and widely applicable
information criterion in singular learning theory. Journal of machine learning research, 11(12).
[6] Aoyagi, M., Watanabe, S. (2005). Stochastic complexities of reduced rank regression in Bayesian estimation.
Neural Networks, 18(7), 924-933.

2データ量子統計の提案
放送大学大学院

池上 孝則

１. 最尤法の信頼性
外れ値や異常値を多く含むデータから統計量を安定的に算出するには，その方法論を適格に選択する必要がある。二次元データにおいて回帰式を決定する場合を例に挙げると，残差に着目した演算である最小
二乗法や最小絶対値法は外れ値等の影響を受けやすく処理に細心の配慮が必要であるが，母数の尤もらしさを示す尤度に着目した処理である最尤法は，求める情報に対する直接的な処理であって残差に着目した演算ではないため，不確かさを内包するデータ群の処理においても有効に機能する可能性がある。
しかし最尤法には，こうした期待感がある反面，数理の常識に照らして腑に落ちない点も少なくない。
に従うと仮定し，その確率密度関数 p(x| の逆関数を尤度と解釈した場

その一例として，データが合を考える。ここに，データ群における の尤度の同時確率によって を推定する尤度

関数Ｌを定義した場合，簡単のために係数を無視すると，
θ

＝

}

＝
＋

exp

＝

(1)

となる。式(1)は，その第１項は定数であることから，／ を中心とした分散が／n の正規分布であることを表している。すなわち，当該尤度関数の最尤推定値は標本平均に一致するが，その結果は平均値の妥当性を再確認するという消極的な情報を提供しているに過ぎない。また，尤度関数の分散が／n であるということはデータが増す毎に分散が減少することを意味しているが，その一方でデータには普遍的に程度のばらつきがあるとする統計量が付随しており，数理的に整合性がとれない。

２．重ね合わせによる尤度関数の検討
最尤法におけるこうした疑義の原因が積統合によって構築する尤度関数にあると推定し，ここでは和結合，すなわち尤度を表す関数の重ね合わせによる情報統合の可能性を検討する。ここに，母集団の分布がN(0,σ2)と仮定し，完全な重ね合わせを行うという意味で正規分布関数
畳み込み関数同士の畳み込み積分を行なうと，

は

ｈ(x)＝

＝

＝

＝

＝

σ

となる。
式(2)は

が分散 2σ2 の正規分布であることを示しているが，それと同時に同時に h(x)を二乗すると母集団の確率密度関数 f(x)に回帰するということをも意味している。この事実は，量子論において「波動関数の絶対値の二乗は電子の存在確率を表す」とする「ボルンの確率解釈」に通じるものがある。さらに，
尤度を虚数と仮定すれば尤度関数という虚の世界から確率密度関数という実の世界に回帰することになり，量子論の解釈に重ねることができる。こうした連想が以下に述べるデータ量子統計の骨格を作っている。
３．データ量子統計の世界観
データ量子統計とは，前述の正規分布同士の畳み込み積分の結果と量子論との相似性を踏まえ，尤度の処理を目的として構築した方法論である。ここに，量子論は「エネルギーには最小単位が存在し，その整数倍で分配される」とする革命的なアイデアに基づいて構成されているが，「データ量子」は「データは本1

来的にその信頼性に関して上限を有する情報」と捉えて導入した量であり，「個々のデータが有する不確かさの程度はそれらが属するデータ群の分散σ2 程度」と解釈する。
また，前述の畳み込み積分の結果を踏まえ，データ量子を表す関数としてデータ群を近似する正規分布関数の逆関数を対応させる。この際，尤度には規格化条件が課されないことから，正規分布関数における規格化の為の定数は省略する。更に，量子論における「波動関数の確率解釈」との整合性を図るという観点に加え，確率と尤度を明確に区別するために，尤度は純虚数として位置付ける。
すなわち，を仮定したデータ群

において，データ の尤度を表すデータ量子を とするとき，データ量子関数 を＝と定義する。ちなみに，データ量子関数は最尤法における尤度関数に相当する関数である。また，データ量子関数の絶対値の二乗で構成される関数Ψ，すなわち
Ψ＝を定義し，それを「データ密度関数」と呼ぶこととする。データ密度関数Ψは母世界における尤度の状態を表すデータ量子関数を実世界におけるデータの存在確率を表す関数に変換したものであり，離散的なデータを微積分が可能な連続関数で処理できることから，記述統計や検定の分野での活用が期待される。
４．データ量子統計に基づく統計量
データ量子統計における代表値は最尤法やベイズ推定と同様でデータ量子関数が最大となる最尤推定値
を代表値とするが，分散の定義は統計学とは異なる。すなわち，データ量子関数 の最大値がとき， と ＝

･exp(＝

（データ密度関数の場合Ψと ＝

･exp(

）
）との交点を

，

である
とすると，

＋

で定義し，混乱を避けるために「データ量子分散」と呼ぶことにする。なお，最尤推定値とデータ量子分
散は収束演算によって同時に決定され，共に外れ値や異常値に対して頑健な統計量である。
相関係数は，まず二次元データにおける 方向， 方向の最尤推定値およびデータ量子分散，を求め,次に共分散を求めるデータ量子関数 ，すなわち＝iに基づいて「データ量子共分散」
回帰式は，を算出し，相関係数ｒである場合には回帰係数はｂ

／

／

を求める。

，

で決定する。一般的な回帰式の場合は，まずデータ群と回帰式の関係を評価するデータ量子関数 を設定する。すなわち，点ら回帰式に下したｘ軸方向・ｙ軸方向の残差を＝，かとするとき，両者に比例関係がない場合は＝でデータ量子関数 を定義し， が最大になるように回帰係数を決定する。
５．まとめ
データ量子統計は，尤度は確率とは異質の情報との基本認識に基づき，量子論を手掛かりとして確立した方法論である。処理対象とするデータに統計モデルを仮定する必要がなく，不確かさを含むデータであっても画一的な手続きによって実数で信頼性の高い結果を安定的に導くことができる等の特徴を有するが，論理性に関する検証はまだ不十分な段階であり，またデータ処理の大半は収束演算を多用する数値演算であることから，先ずは方法論の実効性を重視する理工系の研究者での活用をお勧めしたい。
2

Semiparametric density estimation with Bregman divergence
千葉大・融合理工学府

松野 大輔

千葉大・理学研究院

内藤 貫太

はじめに：密度推定におけるセミパラメトリックな推定量を考える．提案する推定量は，パラメトリック密度推
定量とノンパラメトリックに構成した調整項との積として定義される．調整項は，関数間の距離尺度 Bregman
ダイバージェンスを局所化した基準の最小化を通して導出される．このようにして得た推定量の漸近的性質と実装について報告した．

Bregman ダイバージェンス：非負値関数 f , g 間における Bregman ダイバージェンスは，狭義凸関数
U : [0, ∞) → R，その凸共役 U ∗ (s) := supz∈R {zs − U (z)}，u := U ′ を用いて
Z
DU (f, g) :=[U (f (t)) − U (g (t)) − u (g (t)) {f (t) − g (t)}] dt
RdZ=[U ∗ (u (g (t))) − U ∗ (u (f (t))) − f (t) {u (g (t)) − u (f (t))}] dt
Rdと定義される．例えば，U (x) を x2 /2 や x log x − x とすれば，DU (f, g) は f , g 間の 2 乗距離や Kullback-Leibler (KL) ダイバージェンスとなる．また，U を β-function
Uβ (x) :=xβ+1x− , β ∈ (0, 1]
β (β + 1) β(1)とすると，パラメータが β = 1 と β ↘ 0 となるとき，DUβ (f, g) はそれぞれ 2 乗距離と KL ダイバージェンスに一致する．この Bregman ダイバージェンスを，密度の推定点 x ∈ Rd でZ
DU (f, g | x) :=Rd

KH (t − x) [U ∗ (u (g (t))) − U ∗ (u (f (t))) − f (t) {u (g (t)) − u (f (t))}] dtのように局所化する．ここで，KH (z) = |H|−1/2

K H −1/2 z ，K : Rd → [0, ∞) は原点対称な核関数，バンド幅行列 H は d × d 正定値対称行列である．
設定：n 個の Rd 値観測値 X1 , . . . , Xn が独立に密度 f (x) に従うとする．この密度 f (x) を観測値に基づいて推定する．提案する密度推定量は次の 3 段階で構成される：⟨Step 1⟩ パラメトリックモデル g (x, θ) とパラメータの推定量 θb を与える．ただし，θ ∈ Θ ⊂ Rp , p ≥ 1 と
する．






 

⟨Step 2⟩ 調整項 ξb(x) := arg minξ Ln x, ξ, θb を導入する．ここで，Ln x, ξ, θb は，DU f, g ·, θb ξ xの推定量の ξ に依存する項


Ln
である．

Z

b
x, ξ, θ :=n
    
 
 
1X
b
KH (t − x) U u g t, θ ξ dt −KH (Xi − x) u g Xi , θb ξ
n i=1
Rd
∗



⟨Step 3⟩ セミパラメトリック推定量 fb(x) := g x, θb · ξb(x) を得る．

漸近的結果：g (x, θ) による f の最良近似を与えるパラメータを θ0 とする．n → ∞ において θb − θ0 =√OP (1/ n) とし，g0 (x) := g (x, θ0 ) と略記する．θb 以外に，推定対象の f (x) と，提案推定量 fb(x) を定める
g, KH , U に対しても，収束や滑らかさ，可積分性についての仮定を置いた下で，fb(x) は次のような漸近正規
性を持つ．

(

q
n |H|

1/2

fb(x) − f (x) −µ2
2g0 (x))2 tr [HB (x)]d2
→ N 0, ∥K∥2 f (x) , as n → ∞, H → O.

ここで∂u′′ (f (x))
2f (x) g0 (x) ·B (x) := g0 (x) Hg0 ·f (x) − f (x) Hg02 (x) + 2 · ′u (f (x))
∂x



f (x)
g0 (x)·∂g0 (x)∂x⊤であり，Hh は関数 h の Hesse 行列とする．
実装：提案推定量 fb(x) の実装にあたって，g を正規密度，θb を平均と共分散行列の最尤推定量，K をガウシアンとする．U を (1) の Uβ として導かれた密度推定量

 n−1Zfbβ (x) := g x, θb ·



β
b
K
(X
−
x)
g
X
,
θ
H
i
i
i=1
 1+βKH (t − x) g t, θb
dt

Pn

Rd

(2)

を考える．この fbβ (x) が有限個の正規密度の線型結合として表されることが示される．
シミュレーション：推定対象の密度は，混合正規密度 (Wand and Jones (1993)) や歪正規密度 (Azzalini andCapitanio (2014)) とする．f がこのような混合正規密度，歪正規密度の場合，fbβ (x) の積分二乗誤差が陽に表現されることが示される．この事実に基づき，fbβ (x) を，Hjort and Glad (1995) で提案された関連推定量n
 1X
KH (Xi − x)bbfHG (x) := g x, θ ·n i=1 g X , θbiや，従来のノンパラメトリックカーネル密度推定量n1X
bKH (x − Xi )fKDE (x) :=n i=1と MISE で比較するシミュレーションを行った．また，fbβ (x) のバンド幅選択アルゴリズムについても提案した．このバンド幅選択アルゴリズムにより得られたバンド幅を用いた fbβ の精度についても，シミュレーションにより検証した．
シンポジウムにおいては，(2) の fbβ (x) の平均積分二乗誤差 (MISE) による推定精度を検証するシミュレー
ション結果について報告した．

参考文献
Azzalini, A. and Capitanio, A. (2014). The Skew-Normal and Related Families, volume 3 of Institute ofMathematical Statistics (IMS) Monographs. Cambridge University Press, Cambridge.
Hjort, N. L. and Glad, I. K. (1995). Nonparametric density estimation with a parametric start. Ann.
Statist., 23(3):882–904.
Wand, M. P. and Jones, M. C. (1993). Comparison of smoothing parameterizations in bivariate kernel
density estimation. J. Amer. Statist. Assoc., 88(422):520–528.

一般化 Kendall’s tau 行列による主成分分析
千葉大・融合理工学府 蔣思澄

千葉大・理学研究院

貫太内藤

はじめに：ロバストな主成分分析の手法として，Kendall’s tau 行列に基づく方法が知られている．この方法の理論的性質について調べた先行研究は数多く存在する．Kendall’s tau 行列は，多変量 Spatial Sign によって構成されるが，近年この多変量 Spatial Sign の一般化の研究に進捗が見られた．本報告では，この多変量 Spatial

Sign の一般化を反映した一般化 Kendall’s tau 行列を新たに定義し，その性質について議論した．
設定：X, Y , X1 , . . . , Xn は i.i.d. で期待値ベクトル µ ∈ Rd ，有限な正定値共分散行列 Σ を持つ分布 F に従うとする．d 次元ベクトル y の Spatial Sign は yS(y) = ||y|| 0,y ̸= 0,y=0

で定義される．ここで，||y||2 = y ⊤ y である．母集団での Kendall’s tau 行列 K ∈ Rd×d は

hi⊤
K = E S (X − Y ) S (X − Y )で定義される (Han and Liu (2018))．ξ を非負動径関数とし，これを重み関数とする．このとき，一般化

Kendall’s tau 行列 Kξ は
hi
2Kξ = E ξ (||X − Y ||) (X − Y )(X − Y )⊤ = E [h(X, Y )]で定義される．ここで，h(X, Y ) = ξ (||X − Y ||) (X − Y )(X − Y )⊤ であり，ξ(r) = 1/r のとき，Kξ = K
2である．Kξ のスペクトル分解を Kξ = ΓΛΓ⊤ とする．ここで，Γ = [γ1 . . . γd ] , ΓΓ⊤ = Γ⊤ Γ = Id , Λ = diag (λ1 , . . . , λd ) = 

λ1

O
..

.

O





λd

cξ は
であり，固有値については，単根性 λ1 > λ2 > · · · > λd > 0 を仮定する．Kξ の推定量 K
cξ =
K

n h
i
X
1
2
ξ (||Xi − Xj ||) (Xi − Xj )(Xi − Xj )⊤n(n − 1) i,jで定義される．

b は外れ値の影響を受けやすいので，Σ
b による主成分分析はロバストではないことはよ
問題：標本共分散行列 Σ
く知られている．一方，F が楕円分布のとき，Σ と Kξ が同じ直交行列で対角化できることを示すことができcξ の固有ベクトルが，Σ
b に基づく主成分分析に代わるロバストな主成分分析を提供するこる．このことから，Kとが期待される．本報告では Raymaekers and Rousseeuw (2019) で紹介された重み関数 ξ と K に対応する重み関数を含む 6 つの ξ の集合 Ξ に対して，ξ の選択基準を構築し，どの ξ がどのような意味で良いのか調べた．
特に，本報告では固有空間の類似度と Kξ の固有ベクトルの影響関数に基づく ξ の選択基準を考えた．
固有空間の類似度に基づく選択基準：d × d 正定値行列 A と B に対して，A と B のスペクトル分解を

A = Y ΛA Y ⊤ と B = W ΛB W ⊤ とし，直交行列を Y = [y1 · · · yd ], W = [w1 · · · wd ] とする．部分行列を
Yk = [y1 · · · yk ], Wk = [w1 · · · wk ], k = 1, . . . , d で定義する．k 次元までの A と B の固有空間の類似度を

1 X X ⊤ 2
1
Sk (A, B) =yi wj = Tr Yk Yk⊤ Wk Wk⊤
k i=1 j=1kkkで定義する (Marden (1999))．この定義より，0 ≤ Sk ≤ 1 である．または Sk (A, B) = 1 のとき，A と B が同cξ の類似度について，n → ∞ のとき，
じ固有空間を持つことがわかる．Kξ と K

Sk





k
d
X
X

cξ = 1 − 1 Kξ , Knk i=1√が得られる．ここで，W1 = (1/ n)γj⊤ W1 γi

j=k+1

!2 

λj − λi + oP

 
1
n

Pn

i=1 {2EY [h (Xi , Y )] − 2Kξ } である．データに基づき ξ を選択する

場合，



cξ
Gk Kξ , K



1X X
=
k i=1
k

d

j=k+1



γj⊤ W1 γiλj − λi!2cξ の変動を考慮し，k を固定したとき，最適な重み ξbopt をを考えることにする．特に，Gk Kξ , K
h 
i
cξ
ξbopt = arg min E Gk Kξ , Kξ∈Ξ

と定め，ξ を選択する．
影響関数に基づく選択基準：X を X1 , . . . , Xn から得られた標本平均となる．M (y) = y − X⊤b −1 y − XΣとする．このとき，

Xlarge = arg

maxX∈{X1 ,...,Xn }M (X)と定める．このとき，ロバストネスの観点から，Kξ の第一固有ベクトル γ1 の影響関数 (Hampel et al. (2005))
である IFξ (x, γ1 , F ) を x = Xlarge で評価した ||IFξ (Xlarge , γ1 , F )|| を用いる．特に，

ξ opt = arg min ||IFξ (Xlarge , γ1 , F )||ξ∈Ξで，ξ を選択する．
適用例：ξ の選択に関連するシミュレーション実験の結果および実データへの適用結果については，シンポジウ
ムにて報告した．

参考文献
Hampel, F., Ronchetti, E., Rousseeuw, P., and Stahel, W. (2005). Robust Statistics: The Approach Basedon Inﬂuence Functions. John Wiley and Sons, Inc.
Han, F. and Liu, H. (2018). Eca: High-dimensional elliptical component analysis in non-gaussian distributions. Journal of the American Statistical Association, 113(521):252–268.
Marden, J. I. (1999). Some robust estimates of principal components. Statistics and Probability Letters,
43(4):349–359.
Raymaekers, J. and Rousseeuw, P. (2019). A generalized spatial sign covariance matrix. Journal ofMultivariate Analysis, 171:94–111.

3-step 単調欠測データにおける T 2 型検定統計量を用いた
平均ベクトルの検定
東京理科大学大学院 理学研究科
東京理科大学 理学部
東京理科大学 理学部

1

森本 直紀
八木 文香
瀬尾 隆

はじめに

3-step 単調欠測データにおける平均ベクトルの検定問題 H0 : µ = µ0 vs. H1 : µ ̸= µ0 について考えた．ただし，µ0 は既知ベクトルとした．平均ベクトルに対する T 2 型検定統計量 Q はd µ]−1 (bQ = (bµ − µ0 )′ Cov[eµ − µ0 ) =

3
∑

Qi

(1)

i=1

と表され，Q を構成する統計量 Q1 , Q2 , Q3 の間に独立性を仮定せず相関を考慮することにより，Q の正確な特性関数の漸近展開及び近似上側パーセント点を導出した．なお，µ0 = 0 としても一般性を失d µ] はb は µ の最尤推定量，µe は Σ を既知としたときの µ の最尤推定量，Cov[eわない．(1) において，µCov[eµ] の推定量である (µ と Σ の最尤推定量については，Jinadasa, and Tracy (1992) を参照)．また，
Fujikoshi (2000) に基づいて，Q のカイ二乗近似を改良する変換統計量を提案する．単調欠測データにおける平均ベクトルの検定問題については，統計量の間の相関を考慮すると計算が複雑になるため，独立性を仮定したもとでの研究が進められてきた．相関を考慮したときの検定統計量 Q の正確な帰無分布に関しては，Yagi, Seo, and Hanusz (2019) で 2-step の場合まで議論されており，本研究はその研究を
3-step の場合に拡張したもので，より正確な検定統計量の帰無分布を与えた．
3-step 単調欠測データを以下のように設定した．


x′12
x′13
x′11


 ..
..
..  n .
.
.  1
 x′′′
xN1 3 xN1 2
 N1 1 ′
xN1 +1,1 x′N1 +1,2
∗ 


..
..  n ..
 .
.
.  2


′
 x′
xN2 2
∗ 
 N2 1

 ′

∗
∗ 
xN2 +1,1

 .
..
.. 
 .
 n3
.
. 
 .

′
xN3 1
∗
∗
| {z } | {z } | {z }


p1

p2

p3

i.i.d

xi1 ∼ Np1 (µ1 , Σ11 ), i = 1, 2, · · · , N3 ,i.i.dxi2 ∼ Np2 (µ2 , Σ22 ), i = 1, 2, · · · , N2 ,i.i.dxi3 ∼ Np3 (µ3 , Σ33 ), i = 1, 2, · · · , N1 ,i.i.d

xiP3 = (x′i1 , x′i2 , x′i3 )′ ∼ NP3 (µ, Σ), i = 1, 2, · · · , N1 ,i.i.d

xiP2 = (x′i1 , x′i2 )′ ∼ NP2 (µ12 , Σ2 ), i = 1, 2, · · · , N2 ,
Pj =

j
∑
i=1

pi , N j =j∑ni , j = 1, 2, 3,

i=1

ただし，pi , ni (i = 1, 2, 3) をそれぞれ単調欠測データの各 step の次元数及びデータ数とした．また，
3-step 単調欠測データの平均ベクトルと分散共分散行列をそれぞれ

 
(
)
µ1
µ12
 
µ = µ2  =
,
µ3
µ3
とした．



)
(
Σ11 Σ12 Σ13
Σ2
Σ12(3) ′Σ = Σ12 Σ22 Σ23  =Σ′12(3) Σ33
Σ′13 Σ′23 Σ33

2

検定統計量について
統計量 Q1 , Q2 , Q3 の間の相関を考慮すると Q の特性関数の漸近展開としてP− 23

ϕ(t) = (1 − 2it)

P3
1 ∑
−3
+
βj (1 − 2it)− 2 −j + O(n1 2 )
n1

2

(2)

j=0

が得られ，この結果を用いて以下の Q の近似上側 100α% 点 qAEC (α) と Q のカイ二乗近似を改良する変換統計量 QY C , QZC を提案した．

)
β2
2 (
β0 −
χ2P3 (α) χ2P3 (α),
n1 P3P3 + 2)(
1
QY C = (n1 γ1 + γ2 ) log 1 +
Q ,
n1 γ1
{
(
)}
1
QZC = (n1 γ1 + γ2 ) 1 − exp −Q
n1 γ1

qAEC (α) = χ2P3 (α) −

(3)
(4)
(5)

ここで，(2) ∼ (5) の各式において

1{p2 (4p1 + p2 − 2p3 + 2) p1 (p1 − 2p2 − 2p3 + 2) }p3 (4p1 + 4p2 + p3 + 2) ++,4
1 + r2
1 + r2 + r3
p2 (p1 − p3 ) p1 (p2 + p3 )β1 = p3 (p1 + p2 ) +−,1 + r2
1 + r2 + r3
1{
p2 (p2 + 2p3 + 2) p1 (p1 + 2p2 + 2p3 + 2) }β2 =p3 (p3 + 2) ++,4
1 + r2
1 + r2 + r3
P3 (P3 + 2)
(P3 + 2)β02β0γ1 =, γ2 =
=
γ1 ,
4β2
2β2
P3
n3
n2
, r3 =r2 =n1n1
β0 = −であり，χ2P3 (α) は自由度 P3 のカイ二乗分布の上側 100α% 点である．(4) 及び (5) の変換統計量は，(2)で与えられる Q の正確な特性関数から得られるものであり，極限 (n1 −→ ∞) を考えたときに，Q の分−2布関数は O(n1−1 ) となるのに対して，変換統計量の分布関数には O(n−11 ) の項は現れず O(n1 ) となるため，Q と比較して，より小さな n1 でカイ二乗分布へ収束する．
最後に，いくつかのパラメータの下でモンテカルロ・シミュレーションを行い，検定統計量及び変換
統計量のカイ二乗近似の精度などについて数値的評価を行った．また，数値例を用いることにより，本研究の結果を実データに適用できるかを確かめた．
参考文献[1] Fujikoshi, Y. (2000). Transformations with Improved Chi-Squared Approximations. Journal ofMultivariate Analysis, 72, 249-263.
[2] Jinadasa, K. G., and Tracy, D.S. (1992). Maximum Likelihood Estimation for Multivariate Normal
Distribution with Monotone Sample. Communications in Statistics - Theory and Methods, 21, 4150.
[3] Yagi, A., Seo, T., and Hanusz, Z. (2019). Improved Simplified T 2 Test Statistics for a Mean Vector
with Monotone Missing Data. Communications in Statistics - Simulation and Computation, 48,
1606-1622.

Modal Interval Regression for Visualizing
Data-frequent Areas
Sai YAO(Tohoku University)∗1
Yuko ARAKI (Tohoku University)∗2
January 26, 2024

1 Introduction
Modal interval [1, 2] is a statistic representing the spread of data, which can show
the area where observations appear frequently and robust to distribution. Regression analysis is often used to analyze the continuous relationship between two randomvariables. A modal interval regression method for visualization of bivariate data is
proposed in reference [2]. In this presentation, we improved the method in reference [2]by using kernel density estimation instead of bin segmentation. Numerical experiments
show that, compare with the method in [2], proposed method improves the accuracyof the results.

2 Preliminaries
2.1 Modal Interval
For every α ∈ (0, 1), the 100α% modal interval [1, 2] is defined by
Z xR
 L

R
MIα := mα , mα :=argminxR − xL s.t.
fX (x) dx = α.
[xL ,xR ]⊆cl(xinf ,xsup )xL

2.2 Quantile Regression
Asymmetric absolute value function is defined by(pt
if t ≥ 0,
Jp (t) :=−(1 − p)t if t < 0,with p ∈ (0, 1). The quantile regression is formulated as the following optimization
problem [3]NX
minimizeJp (yi − rθ (xi )),
(1)
θ

i=1

where (xi , yi ) is observation, function rθ (x) is regression model, and vector θ denotes
adjustable parameters to be optimized.
2.3 Spline Smoothing
The spline function is a piecewise polynomial and possesses certain-times continuous
differentiability [4]. Spline smoothing is to solve the following optimization problem
ZNX
wi L(yi − s(xi )) + λ |s′′ (x)|2 dx,(2)minimizeρs∈Sd (⊔b )

∗1

i=1

I

e-mail: yaosai@dc.tohoku.ac.jp
e-mail: yaraki@tohoku.ac.jp
This work was partly by JST SPRING, Grant Number JPMJSP2114.
Keywords: Modal Interval, Regression Analysis, Spline Smoothing, Data Visualization.
∗2

where Sdρ (⊔b ) is the set of all univariate spline functions of degree d at most andsmoothness ρ with b polynomials, wi > 0 is weight, L : R → [0, ∞) is loss function, λ >0 is the regularization parameter to balance bias and variance to prevent overfitting.
2.4 Generalization of Nadaraya-Watson Estimator
The generalization of Nadaraya-Waston estimator can be used to estimate conditional modal interval [1, 2]. Conditional cumulative distribution function of Y given
X = x can be estimated byPN
Kh (x − xi )I(yi ≤ y),(3)
F̂Y |x (y) := i=1PN
i=1 Kh (x − xi )where Kh (·) is a kernel function with bandwidth h, and I(yi ≤ y) = 1 when yi ≤ y,otherwise I(yi ≤ y) = 0. The 100α% conditional modal interval for given X = x is
estimated asargmin yup − ylow s.t. F̂Y |x (yup ) − F̂Y |x (ylow ) ≥ α.
(4)[ylow ,yup ]3 Spline Modal Interval Regression
Our proposed method is formulated as the following optimization problem
Z
NX
minimizewi Jp̂⟨i⟩ (yi − s(xi )) + λ |s′′ (x)|2 dxlow¯
s,s̄∈Sdρ (⊔b )
I ¯
i=1
¯
Z
N
X
+
wi Jp̂⟨i⟩(yi − s̄(xi )) + λ |s̄′′ (x)|2 dx
up
i=1

(5)

I

subject to ∀x s(x) ≤ s̄(x),
¯
⟨i⟩
⟨i⟩
where p̂low = F̂Y |xi (ylow ) and p̂up = F̂Y |xi (yup ), whichqare estimated with the general5
ization of Nadaraya-Waston estimator in (4). wi = fˆ(xi ) and fˆ(xi ) is estimated bykernel density estimation method. s(x) and s̄(x) are lower curve and upper curve of¯
conditional modal interval, respectively. The condition ∀x s(x) ≤ s̄(x) is to ensure¯
that the two curves will not cross each other.

4 Numerical Experiments
In numerical experiments, we estimate the 50% conditional modal interval of several
simulation datasets with different distributions and sample sizes. Numerical experiments show that the results of proposed method in (5) have higher accuracy than the
results of method in [2], since proposed method does not rely on bin segmentation. We
also used proposed method to analyze the relationship between temperature and water
vapor pressure in meteorological data [5]. Compared with general regression methods,
the propsed method can not only visualize the relationship between two variables, butalso visualize the changing trend of the spread of conditional distribution.

References[1] Jan G.De Gooijer and Ali Gannoun, Nonparametric conditional predictive regions for time series, Computational Statistics & Data
Analysis, 33, No.3, (2000), 259-275.
[2] S. Yao, D. Kitahara, H. Kuroda and A. Hirabayashi, Modal Interval Regression Based on Spline Quantile Regression, IEICE
Transactions on Fundamentals of Electronics, Communications and Computer Sciences, E106.A, No.2, (2023), 106-123.
[3] R. Koenker and G. Bassett Jr., Regression quantiles, Econometrica, 46, No.1, (1978), 33–50.
[4] S. Wold, Spline functions in data analysis, Technomet., 16, No.1, (1974), 87–111.
[5] Meteorological data of Beijing, China, from August 1 to 7, 2021, China Meteorological Data Service Centre, http://data.cma.cn/en.

Characterization of Asymptotic Distribution in Martingale
Convergence of the Supercritical Branching Processes
with Poissonian Offsprings
Ankhbayar Bat-Erdene1 and Shuji Kawasaki2
1 Graduate School of Science and Engineering, Iwate University
2 Faculty of Science and Engineering, Iwate University

Abstract. A branching process (BP) {Zn | n ∈ N0 } is a stochastic model of population increaseby the proliferation of individuals. A normalized BP forms a martingale and converges, as n → ∞,almost surely. The asymptotic distribution on (0, ∞) of the martingale convergence is suggested to
be used in various applications such as the potential theory of BP. Nevertheless, the explicit expression of the asymptotic distribution is unknown. We propose a characterization of the distribution
based on the analysis of a functional equation for the Laplace transform (LT) of the distribution.
Particularly, we present the Laguerre expansion of the asymptotic distribution. It turns out that,by numerical validation, our proposed model agrees well with the true distribution.
1. Background and Problem
Let X be a random variable (r. v.) taking values in N0 = {0, 1, 2, · · · }, the set of nonnegative
integers, and let p = {pk | k ∈ N0 } be its distribution. Then, a branching process (BP) {Zn | n ∈ N0 }is defined by
Zn−1
X
Z0 = 1,Zn =Xn,i , n ∈ N,
i=1

where {Xn,i | i ∈ N} is an independently identically distributed sequence of r.v.s with Xn,i ∼ X foreach n ∈ N0 . That is, each of the i-th individuals in (n − 1)-th generation produces Xn,i offspring
according to the same distribution p independently of each other to form the n-th generationindividuals. Because Z0 = 1, Z1 ∼ X especially. For a BP, the extinction and survival probabilities
are, respectively, defined by(1.1)

η = P(∃n : Zn = 0)and

ζ = P(Zn > 0, ∀n ≥ 0) = 1 − η.

If λ = E[Z1 ] ≤ 1, the population dies out with extinction probability η = 1. If µ > 1, the process iscalled supercritical. In this latter case, the extinction probability satisfies η < 1. Let Wn = µ−n Zn ,n ∈ N0 . Then, {Wn } is known [1, Theorem I.6.1] to form a martingale and converges to an r. v.
W∞ on R+ = [0, ∞) almost surely as n → ∞.
W∞ is known to have a probability density function (p. d. f.) on (0, ∞) [1, Corollary I.12.1],
which we denote by w : (0, ∞) 7→ R+ , while W∞ has a point mass at the origin, P (W∞ = 0) = η[1, Theorem I.6.2]. Thus, we may write W∞ ∼ ηδ(x) + w(x) on x ∈ R+ , where δ(x) denotes the
Dirac delta [2, Theorem I.8.3].
However, not much is known about the p. d. f. w(x) so far [3, p. 97], despite its importance in
BP. Because the expression of w(x) is unknown, further theory of BP has been developed withthe unknown w(x) as it is (see e. g., [1]). The objective of this study is to characterize the p. d. f.
Particularly, we identify w(x) in the form of Laguerre expansion.

2. Characterization and Main Results
We assume p is Poissonian and let f : [0, 1] 7→ [0, 1] be the probability generating function
for p . Then,f (s) = eλ(s−1) , s ∈ [0, 1]. Furthermore, let φ : R+ 7→ R+ be the LT of W∞ :R
φ(u) = (0,∞) e−ux w(x) dx. For φ and f , the following functional equation, called Abel’s equation,holds [1]:  h u iu(2.1)
φ(u) = f φ= exp λ φ,u ∈ R+ .
λλEspecially, taking u = 0 or ∞ in (2.1) corresponds in general to the two solutions of the equations = f (s): s = 1 = φ0 (0) or s = η = φ0 (∞) [1, Corollary I.10.3], respectively. Thus, McLau(n)P
rin’s expansion of the LT φ(u) = n∈N0 φ n!(0) un . We consider McLaurin’s coefficient φ(n) (0) andcharacterize the distribution w(x).
The following are the main results of the study.
Proposition 1. φ(n) (0) has the following form:
φ(n) (0) = −Kn (−ρ)n−1 ,

(2.2)

for

n ∈ N0 ,λ, and Kn = Kn (λ) is a positive-valued function defined recursively by
λ−1

n−1 
1 X n − 1 −(n−1−l)
λKn−l (λ)Kl (λ),
for n ≧ 2;
K1 (λ) ≡ 1,
Kn (λ) =lρn−1

where ρ ≜
(2.3)

l=1

The function Kn is irreducible of ρ, i.e., Kn is not a multiple of ρ.
Theorem 2. The p. d. f. w(x) has the Laguerre expansion
X(2.4)
w(x) = e−x
qr Lr (x),x ∈ R+r∈N0

in L2 (R
(2.5)

2+ ), for a real-valued sequence {qn } ∈ ℓ (N0 ) given byr  Xr Km+1(−ρ)m .
qr =
m!
m
m=0

The constructed w(x) is indeed well-defined in L2 (R+ ) as can be seen from the following proposition. In the following, let q n = [ q1 , q2 , · · · , qn ]⊤ .
Proposition 3.
(2.6)

sup ∥ q n ∥ℓ2 < ∞.

n∈N

References[1] Athreya, K. B., and Ney, P. E., Branching Processes, Dover Publ., 2004.
[2] Harris, T. E. The Theory of Branching Processes, Dover Publ., 1963.
[3] Hofstad, R., Random Graphs and Complex Networks, Cambridge Univ Press, 2016.

Application of the Topological Data Analysis to Cosmology(位相的データ解析の宇宙論への応用)

竹内 努

1,2

3

, 本武陽一 , 西澤 淳 4,5,6, 村上 広椰 1, 加納 龍生 1, 馬 海霞 1, 内田 舜也 1,
Suchetha COORAY7, 河野 海 1

1. 名古屋大学理学研究科素粒子宇宙物理学専攻, 2. 統計数理研究所統計的機械学習研究センター,
3. 一橋大学ソーシャル・データサイエンス研究科, 4. 岐阜聖徳学園大学 DX 推進センター 5. 名古屋大
学高等研究院, 6. 名古屋大学素粒子宇宙起源研究所, 7. 国立天文台科学研究部

銀河とは, 星と星間物質(ガスと星間塵), 暗黒物質からなる巨大な天体である. 銀河は現在観測可能な宇宙に数千億個存在しており, 我々の目に見える波長(可視光線)での宇宙の姿を形作っている. しかし宇宙誕生当時, 物質分布はほぼ一様であった. 銀河は平均からわずかに密度が高い領域が重力によって成長し, 合体成長を経て現在の姿へと進化してきたのである(図 1). 銀河進化は周囲の銀河の密度やガス密度など, 銀河の置かれた環境にも強く依存する極めて複雑な過程である. そして銀河の 図 1: 一様な物質分布から銀河が形成さ形成と進化は 138 億年の宇宙進化の歴史の中でも非 れ、進化する様子.
常に重要な現象のひとつであると考えられている.
我々は 130 数億年にわたる銀河の形成・進化という複雑な物理現象を, これまでの宇宙物理学の方法とはまったく異なる, データ科学の最新手法である位相的データ解析(topological data analysis: TDA)に基づくアプローチによって新たな角度から解明することを試みている. 本研究は特に, 銀河形成以前の通常物質(宇宙物理学ではバリオンと呼ぶ)が
光 と の 相 互 作 用 に よ っ て 生 じ る 粗 密 波 で あ る 「 バ リ オ ン 音 響 振 動 」 (baryon acoustic
oscillation: BAO)が銀河分布に残す痕跡について扱う. バリオンと光の相互作用が切れる宇
宙の晴れ上がり(脱結合)時に物質が放射した熱放射は宇宙膨張によって波長が引き伸ばされ, 現在は電波(ミリ波)の波長で宇宙マイクロ波
背景放射(cosmic microwave background: CMB)として観測される. BAO は CMB のゆらぎに特徴的なスケールとして現れ, ゆらぎをフーリエ変換
したパワースペクトルにおいて最も顕著なピー
クとして観測される(図 2). このスケールは共動距離で観測した際に, どの赤方偏移においてもおよそ 148 Mpc (1 Mpc = 3.08×1024 cm)の一定値
図 1: Planck 衛星が観測した CMB のパワースをとる. この性質から, BAO シグナルの観測的ペクトル. 最も顕著なピークが BAO である.

スケールは宇宙論パラメータの推定の最も信頼できる指標のひとつとなっている.
従来, BAO の検出には銀河分布の 2 点相関関数やパワースペクトルが用いられてきたが,BAO シグナルは銀河の重力的クラスタリングによって生じるシグナルに対して 1/100 程
度と微弱であるため, 統計的に有意な結果を得るためには大規模かつサンプリングレートの高い銀河探査が必要であり, 現存する最大の銀河探査であるスローンデジタルスカイサーベイ(Sloan Digital Sky Survey: SDSS)によって初めてそのシグナルが検出された. この規模の探査は決して容易ではなく, 独立な, そしてより有効な検出方法が望まれている.
我々は, 銀河の空間分布がもつ位相幾何学的情報を抽出することで BAO シグナルの検出や宇宙大規模構造を解析することを目的として, パーシステントホモロジー(PH) の計算を行った. PH はトポロジカルデータ解析(TDA) の中心的方法論であり, ホモロジー群に対して関心のあるパラメータに関するフィルトレーションによって定義される. さらに,穴の位置や形状などの宇宙論的に重要な情報を得ることができる点が特徴である. 本研究では, 共動距離(宇宙膨張の影響を取り除いた本質的な距離)をパラメータにとることにより, 空間に対する解像度を変化させた際に生じる構造を次元ごとに解析することができる.
我々はまず TDA の BAO 検出
に対する有用性の検証のため, バ
リオン有り・無しの二通りの N 体
シミュレーションデータに対してPH の計算を行った. その結果, バ
リオン有の場合にのみ有意な 2 次ホモロジー(H2), すなわち中空のシェル構造を検出した(図 3). また,

図 3: 宇宙論的シミュレーションデータから求めたパーシステン 検出された H2 構造の半径の平均
ト図. 左がバリオン有, 右がなしの場合を示す.

値 は 𝑟𝑟̅ = 150.16± 8.46 [Mpc] で あり, BAO の理論と整合的であった. さらに, SDSS eBOSS サーベイによって得られた実際のクエーサー観測データに対しても同様の 解 析 を 行 い , H2 に 対 し て 𝑟𝑟̅ =146.6± 2.0 [Mpc]という推定値を得た.
また本研究では PH の逆解析も行い,得られた H2 に対応する実際の宇宙の構造を可視化することに成功した(図4). 銀河分布に対しては世界初の例である. これらの計算結果に H1, H0 についての考察を加味し, 銀河の空間分布 図 4: SDSS eBOSS データから求めた PH (左), およびそに関して得られた結果を報告する.

の逆解析によって得られた銀河分布の対応構造.

代表的な接合関数の裾従属係数と裾次数*
吉羽 要直†

東京都立大学 大学院経営学研究科

金融機関のリスク管理においては，代表的な資産価格変動や信用状態を𝑑変量のリスクファクターとして，その𝑑変量の従属関係を正規接合関数や t 接合関数といった楕円接合関数で捉えることが多い．しかしながら，正規接合関数は変量間の裾従属性が弱く，t 接合関数はリスクファ
クターの上下の変動に対して対称な従属関係しか表現できない．そのため，正規接合関数や t 接合関数に非対称性を入れたモデリングも発展してきている．本研究では，実務的に用いられることが多い代表的な接合関数の裾従属性について，𝑑 = 2変量に焦点を当てて整理する．具体的には伝統的に用いられている裾従属係数とともに，弱い裾従属関係を捉える裾次数について，Joe(2014)，吉羽(2020, 2021)，Yoshiba, Koike and Kato (2023)の研究の成果に沿って整理する．
一般に，2 変量確率ベクトル𝑿 = (𝑋1 , 𝑋2 )⊤の同時分布関数𝐹(𝑥1 , 𝑥2 )は，各変量が連続確率変数
であれば，各周辺分布関数𝐹1 (𝑥1 ), 𝐹2 (𝑥2 )と接合関数𝐶(𝑢1 , 𝑢2 )を用いて𝐹(𝑥1 , 𝑥2 ) = 𝐶(𝐹1 (𝑥1 ), 𝐹2 (𝑥2 ))(1)と一意に表現できる．確率ベクトル𝑿が相関𝜌の 2 変量標準正規分布に従う場合，内包される接合関数が相関𝜌の正規接合関数𝐶𝑁 (𝑢1 , 𝑢2 ; 𝜌)である．𝑿が相関𝜌，自由度𝜈の 2 変量 t 分布に従う場
合は，相関𝜌，自由度𝜈の t 接合関数𝐶𝑡 (𝑢1 , 𝑢2 ; 𝜌, 𝜈)が内包される接合関数である．確率ベクトル𝑿に非対称性を入れると，その同時分布に内包される非対称接合関数を定義できる．
一方，2 変量の従属関係については，アルキメデス型の接合関数が用いられることもあり，代表的なものに 1 パラメータ𝜃により以下で表せる Gumbel，Clayton，Frank の各接合関数がある．
1⁄𝜃

Gumbel𝐶𝐺 (𝑢1 , 𝑢2 ; 𝜃) = exp [−{(− ln 𝑢1 )𝜃 + (− ln 𝑢2 )𝜃 }Clayton𝐶𝐶 (𝑢1 , 𝑢2 ; 𝜃) = (𝑢1−𝜃 + 𝑢2−𝜃 − 1)

Frank

1⁄𝜃

,

],

𝜃>1

(3)

𝜃>0

1
(e−𝜃𝑢1 − 1)(e−𝜃𝑢2 − 1)
𝐶𝐹 (𝑢1 , 𝑢2 ; 𝜃) = − ln {1 +},𝜃e−𝜃 − 1

(2)

𝜃∈ℝ

(4)

裾従属性は裾次数𝜅𝐿 あるいは𝜅𝑈 が 1 の場合は裾従属係数𝜆𝐿 あるいは𝜆𝑈 で従属性の強さを捉えられ，裾次数が(1,2)の範囲であれば，裾次数で従属性の強さを捉えられる．代表的な 2 変量接合関数について，下側・上側それぞれの裾従属係数𝜆𝐿 , 𝜆𝑈 あるいは裾次数𝜅𝐿 , 𝜅𝑈 は，Joe (2014)の Table
4.1 と Clayton（Mardia–Takahasi–Clayton–Joe）接合関数に関する 4.6.1 節の議論より，自由度𝜈の1 変量 t 分布関数を𝑇𝜈 (⋅)として表 1 のようにまとめられる．
̃ = (𝑍̃1 , 𝑍̃2 )⊤ とそれとは独立な標準正規確率変数𝑍0 ，非相関𝜌̃の 2 変量標準正規確率ベクトル𝒁
対称性を表すパラメータ𝜹⊤ = (𝛿1 , 𝛿2 ) ∈ (−1,1)2を用いて，2 変量確率ベクトル𝒀 = (𝑌1 , 𝑌2 )⊤を
𝑌𝑗 = 𝛿𝑗 |𝑍0 | + √1 − 𝛿𝑗2 𝑍̃𝑗 ,*†𝑗 = 1,2

本研究は JSPS 科研費 JP21K01581 からの助成を受けている．
E-mail: tyoshiba@tmu.ac.jp

1(5)とすると，𝑌1と𝑌2の相関は𝜌となり，𝒀は非対称性𝜹の 2 変量非対称正規分布𝑆𝑁(𝜌, 𝜹)に従い，内包される接合関数として非対称正規接合関数𝐶𝑆𝑁 (𝑢1 , 𝑢2 ; 𝜌, 𝛿1 , 𝛿2 )が定義される．非対称正規接合
関数は，正規接合関数と同様に，上下ともに裾従属係数はゼロで漸近独立であり，裾次数については，Yoshiba, Koike and Kato (2023)で以下のように整理されている．
𝜅𝐿 (𝐶𝑆𝑁 ; 𝛿1 , 𝛿2 ≤ 0) = 𝜅𝑈 (𝐶𝑆𝑁 ; 𝛿1 , 𝛿2 ≥ 0) =

2
1+𝜌

(6)

𝜅𝐿 (𝐶𝑆𝑁 ; 𝛿1 , 𝛿2 ≥ 0) = 𝜅𝑈 (𝐶𝑆𝑁 ; 𝛿1 , 𝛿2 ≤ 0) =2
1 + 𝜌̃

(7)

表 1 代表的な 2 変量接合関数の裾従属係数と裾次数
下側 𝜆𝐿 /𝜅𝐿

接合関数

𝜅𝐿 =

正規

上側 𝜆𝑈 /𝜅𝑈

2
1+𝜌𝜅𝑈 = 𝜅𝐿 =

2
1+𝜌

t

(1 − 𝜌)(𝜈 + 1)
𝜆𝐿 = 2𝑇𝜈+1 (−√)(1 + 𝜌)𝜆𝑈 = 𝜆𝐿

Gumbel𝜅𝐿 = 21/𝜃

𝜆𝑈 = 2 − 21/𝜃Clayton

𝜆𝐿 = 2−1/𝜃𝜅𝑈 = 2Frank𝜅𝐿 = 2

𝜅𝑈 = 𝜅𝐿 = 2

2 変量確率ベクトル𝒀 ∼ 𝑆𝑁(𝜌, 𝜹)と独立にガンマ分布に従う確率変数𝑉 ∼ 𝐺(𝜈/2, 𝜈/2)を用いて，𝑿 = 𝒀⁄√𝑉で Azzalini and Capitanio (2003)の非対称 t 分布に従う𝑑変量確率ベクトル𝑿が構成さ
れ，内包される接合関数として AC 非対称 t 接合関数𝐶𝐴𝐶𝑆𝑡 (𝑢1 , 𝑢2 ; 𝜌, 𝜈, 𝛿1 , 𝛿2 )が定義される．吉羽(2020, 2021)で示されているように，𝐶𝐴𝐶𝑆𝑡 (𝑢1 , 𝑢2 ; 𝜌, 𝜈, 𝛿1 , 𝛿2 )は漸近従属で，裾従属係数𝜆𝐿 , 𝜆𝑈 は非対称性𝛿1 , 𝛿2 によって異なる．
一方，相関𝜌の 2 変量標準正規確率ベクトル𝒁と𝑉 ∼ 𝐺(𝜈/2, 𝜈/2)を用いて，非対称性パラメータを𝜸 = (𝛾1 , 𝛾2 ) ∈ ℝ2 として，𝑿 = 𝜸𝑉 −1 + 𝒁⁄√𝑉により，一般化双曲型（generalized hyperbolic）の
非対称 t 分布に従う𝑑変量確率ベクトル𝑿が構成され，内包される接合関数として GH 非対称 t 接合関数𝐶𝐺𝐻𝑆𝑡 (𝑢1 , 𝑢2 ; 𝜌, 𝜈, 𝛾1 , 𝛾2 )が定義される．この裾従属係数は，𝛾1 = 𝛾2 = 𝛾での特殊ケースで𝛾 < 0ならば𝜆𝐿 = 1, 𝜆𝑈 = 0となり，𝛾 > 0ならば𝜆𝐿 = 0, 𝜆𝑈 = 1で与えられる．裾従属係数がゼロで漸近独立の場合の裾次数の分析は，今後の課題である．
参考文献
Azzalini, A. and Capitanio, A. (2003), “Distributions generated by perturbation of symmetry with emphasis
on a multivariate skew t-distribution,” Journal of the Royal Statistical Society Series B, 65(2), 367–389.
Joe, H. (2014), Dependence Modeling with Copulas, Boca Raton, Chapman & Hall/CRC.
吉羽 要直 (2020),「非対称 t 接合関数の性質と統計的推定方法」『統計数理』68(1), 45–63.
吉羽 要直 (2021),「極値での従属性および非対称性と信用ポートフォリオリスク」『日本統計学
会誌』51(1), 157–178.
Yoshiba, T., Koike, T. and Kato, S. (2023), “On a measure of tail asymmetry for the bivariate skew-normal
copula,” Symmetry, 15(7), 1410.

2

コピュラ・グラフィック推定量- 理論と応用 –

江村剛志*
概 要
生存時間解析の多くの手法は, 打ち切り時間と生存時間が互いに独立であるという仮定（独立打ち切りの仮定）下で妥当である. 一方, この仮定が成立しない「従属打ち切り」に起因するバイアスの問題が生存時間解析の文脈において活発に議論されている. 例えば, 従属打ち切
り下における Kaplan-Meier 推定量のバイアスがよく議論される. 本研究では, Kaplan-Meier
推定量に代わる推定量である, Copula-Graphic 推定量を紹介し, その理論と応用例を示す.
キーワード: Kaplan-Meier 推定量 ▪ コピュラ ▪ Copula-Graphic 推定量 ▪ 生存関数 ▪ ２
標本検定 ▪ ノンパラメトリック検定 ▪ 右打ち切り

1. はじめに「従属打ち切り」とは, 生存時間変数と打ち切り変数が互いに独立でない場合を指す
が, 通常の生存時間解析は, 独立性の仮定を置いて行われる. この「独立打ち切り」の仮定が満たされないような状況は, 生存時間解析の応用において散見される. これに関連する諸問題は生存時間データ解析の文脈において広く議論されてきた（ Braekers and Veraverbeke(2005), Emura and Michimae(2017), Staplin etal.(2015), Xu et al.(2018), Moradian et al.(2019), Schneider et al.(2020,2022), Deresa and Van Keilegom(2020a,b), Huang et al.(2021), Li et al.(2021),
Lo and Wilke(2010), Deresa et al.(2022)). 独立打ち切りの仮定が満足されない場
合, それを基にした Kaplan-Meier 推定量などは一致性を持たない（推定量が漸近バイアスを持つ）ため, 解析結果はバイアスを持つ.
従属打ち切りによる漸近バイアスの研究結果がいくつか報告されている. 従属打
ち切り下で Rivest and Wells(2001)は Kaplan-Meier 推定量の漸近バイアスを解析的に導出した. また Emura and Chen(2016, 2018)は Cox 回帰法の漸近バイアスを解析
的・数値的に調べた. Emura and Hsu(2020)は, 2 標本問題における治療効果の漸近バ
イアスを導出した. 道前(2023)は, 競合リスクによる従属打ち切りが尤度推定量に与えるバイアスをシミュレーションで検証した. Emura et al.(2024)は因子実験計画における治療効果の推定バイアスをシミュレーションで検証した. 上記の全ての研究結
果では「コピュラ」(Nelsen 2006)で生存時間と打ち切り時間の従属構造をモデリングし, バイアスの構造を解析している.
本稿では, 従属打ち切りの影響をコピュラで調節するために導出された CopulaGraphic 推定量を紹介する. この推定量は, 生存時間と打ち切り時間の相関構造をコピュラで正しく特定することにより, バイアスの無い生存関数の推定を行うものであ*統計数理研究所,

Email:takeshiemura@gmail.com

1

る. さらに, それを応用した２標本検定法（Emura and Hsu, 202)）と因子実験計画法（Emura et al.2024）を紹介する.

2. Copula-Graphic 推定量とその応用
本節では, 独立打ち切りのある生存データ解析の基本的な解説と, Copula-Graphic 推
定量の定義を与える. またその２標本問題への応用を示す.

2.1 Copula-Graphic (CG) 推定量
以下で生存関数や打ち切りの基本的な定義を与え, Copula-Graphic 推定量を定義する.
２標本問題を考える. ２つの群のいずれか(𝑖 = 1, 2)に属する個体（𝑗 = 1, . . . , 𝑛 ）の生存時間を 𝑇 で表し, その生存関数を 𝑆 (𝑡) ∶= 𝑃(𝑇 > 𝑡) とする. 𝑛 は群 𝑖 における個体数である. 合計 𝑛 ∶= 𝑛 + 𝑛 個の個体は全て独立であると仮定する. 本
稿では 𝑇 がハザード関数を 𝜆 (𝑡) ∶= 𝑓 (𝑡)/𝑆 (𝑡) を持つとする.
打ち切り時間を 𝑈
データを, その生存関数を 𝐺 (𝑢) ∶= 𝑃( 𝑈 > 𝑢) とおく. ２群の生存; 𝑖 = 1, 2; 𝑗 = 1, . . . , 𝑛 と定義し, ここで 𝑋 ∶= min( 𝑇 , 𝑈 ),𝑋 , 𝛿𝛿 ∶= 𝑰{𝑇 ≤ 𝑈 }, 𝑰{. } は指標関数とする.
また, 𝑌 (𝑠) ∶= 𝑰{𝑋 ≥ 𝑠}, 𝑌 (𝑠) ∶= ∑𝑌 (𝑠) を定義する.

従属打ち切り下で 𝑆 (. ) を推定することを考える. 以下では, 打ち切りの従属構造
をコピュラでモデリングすることにより 𝑆 (. ) を一致推定することを考える.
まず 𝑇 と 𝑈 の同時分布が次のアルキメデスコピュラで特定されるとする:
𝑃 𝑇 > 𝑡, 𝑈 > 𝑢 = 𝜙
ここで 𝜙

𝜙

{𝑆 (𝑡)} + 𝜙{𝐺 (𝑢)} ,𝑖 = 1, 2.

(. ) は生成素と呼ばれ 𝜃 , 𝑖 = 1, 2, はパラメータである. 定義より, 生

成素は連続な減少関数であり, 𝜙(0) = ∞, 𝜙

クレイトン型コピュラで生成素は 𝜙(1) = 0 を満足する.

(𝑡) = (𝑡− 1)/𝜃 , 𝜃 > 0 であり, このコ

ピュラでの相関の強さを表すケンドール相関は 𝜏= 𝜃 /(𝜃 + 2) である. グンベル型

コピュラで生成素は 𝜙(𝑡) = { − log(𝑡) }, 𝜃 ≥ 0 であり, ケンドール相関は= 𝜃 /(𝜃 + 1) である.

𝜏

Rivest and Wells(2001)は次式で Copula-Graphic 推定量を提案した:𝑆(𝑡) ∶= 𝜙𝜙:

生成素を 𝜙
し, 𝑆

,

𝑌 𝑋

−1
𝑛

−𝜙

𝑌 𝑋
𝑛

,

0 ≤ 𝑡 ≤ max 𝑋.

(𝑡) = − log (𝑡) とおくと, アルキメデスコピュラは独立コピュラを生成(𝑡) = 𝑆 (𝑡) = ∏ ;,{1 − 1/𝑌 (𝑋 )} となる. よって, Copula-Graphic 推

定量は Kaplan-Meier 推定量を一般化した推定量である.

2

クレイトン型コピュラの下での Copula-Graphic 推定量は𝑆

𝑌 𝑋(𝑡) = 1 −
:

−1
𝑛

,

−

𝑌 𝑋
𝑛

,

グンベル型コピュラの下での Copula-Graphic 推定量は⎡⎤
𝑌 (𝑋 ) − 1
𝑌 (𝑋 )⎢⎥
𝑆 (𝑡) = exp ⎢−− log− − log
⎥
𝑛
𝑛
⎢
⎥
:
,
⎣
⎦
で与えられる. これら推定量は R パッケージ compound.Cox (Emura et al. 2019a) のR 関数 CG.Gumbel(.), CG.Clayton(.), CG.Frank(.) で計算することが出来る.

この Copula-Graphic 推定量は, モデルが正しく特定されていれば生存関数を一致推定できる. これを応用し, ２標本検定などを行うことができる. また, ３群以上の場合や, ２つ以上の因子を持つ複雑な群構造を構成する分散分析も行うことができる. そのほか, 生存木を構築したり, 予後予測を行う際にも, Copula-Graphic 推定量は有用となる. 詳細は, 講演にて説明する.

参考文献
Braekers, R. and Veraverbeke, N. (2005). A copula‐graphic estimator for the conditional survival
function under dependent censoring, Canadian Journal of Statistics, 33, 429-447.
Deresa, N. W. and Van Keilegom, I. (2020a). Flexible parametric model for survival data subject
to dependent censoring, Biometrical Journal, 62, 136-156.
Deresa, N. W. and Van Keilegom, I. (2020b). A multivariate normal regression model for survival
data subject to different types of dependent censoring, Computational Statistics & Data
Analysis, 144, 106879.
Deresa, N. W., Van Keilegom, I. and Antonio, K. (2022). Copula-based inference for bivariate
survival data with left truncation and dependent censoring, Insurance: Mathematics andEconomics, 107, 1-21.
Emura, T., and Chen, Y.-H. (2016). Gene selection for survival data under dependent censoring,a copula-based approach, Statistical Methods in Medical Research, 25, 2840–2857.
Emura, T., and Chen, Y.-H. (2018). Analysis of Survival Data with Dependent Censoring:
Copula-Based Approaches, Springer.
Emura, T., and Hsu, J.-H. (2020). Estimation of the Mann-Whitney effect in the two-sample
problem under dependent censoring, Computational Statistics and Data Analysis, 150,
106990.
Emura, T., Matsui, S., and Chen, H.-Y. (2019a). compound. Cox: univariate feature selection andcompound covariate for predicting survival, Computer Methods and Programs in
Biomedicine, 168, 21-37.

3

Emura, T., and Michimae, H. (2017). A copula-based inference to piecewise exponential models
under dependent censoring, with application to time to metamorphosis of salamander
larvae, Environmental and Ecological Statistics, 24, 151–173.
Emura, T., Dizhaus M., Dobler, M., Murotani, K. (2024). Factorial survival analysis for treatment
effects under dependent censoring, Statistical Methods in Medical Research, in press.
Huang, X., Wang, W., and Emura, T. (2021). A copula-based Markov chain model for seriallydependent event times with a dependent terminal event, Japanese Journal of Statistics
and Data Science, 4, 917-51.
Li, D., Hu, X. J., and Wang, R. (2021). Evaluating association between two event times withobservations subject to informative censoring, Journal of the American Statistical
Association, doi.org/10.1080/01621459.2021.1990766.
Lo, S. M., and Wilke, R. A. (2010). A copula model for dependent competing risks, Journal of
the Royal Statistical Society: Series C (Applied Statistics), 59, 359-376.
道前洋史 (2023). 「競合リスクを伴う左側切断・右側打ち切りデータの解析～現状と今後の展望について～」『日本統計学会誌』, in press.
Moradian, H., Larocque, D., and Bellavance, F. (2019). Survival forests for data with dependent
censoring, Statistical Methods in Medical Research, 28, 445-461.
Nelsen, R. B. (2006). An Introduction to Copulas, Springer.
Rivest, L.-P., and Wells, M. T. (2001). A martingale approach to the copula-graphic estimator forthe survival function under dependent censoring, Journal of Multivariate Analysis, 79,
138-155.
Schneider, S., Demarqui, F. N., Colosimo, E. A., and Mayrink, V. D. (2020). An approach tomodel clustered survival data with dependent censoring, Biometrical Journal, 62, 157174.
Schneider, S., Demarqui, F., and de Freitas Costa, E. (2022). Free-ranging dogs’ lifetime
estimated by an approach for long-term survival data with dependent censoring.
Environmental and Ecological Statistics, 29, 869–911.
Staplin, N., Kimber, A., Collett, D., and Roderick, P. (2015). Dependent censoring in piecewise
exponential survival models, Statistical Methods in Medical Research, 24, 325-341.

4医学研究における不完全データの感度解析
大阪大学大学院医学系研究科医学統計学教室 服部 聡

E-mail: hattoris@biostat.med.osaka-u.ac.jp
医学研究においては、治療法などの有用性を評価するために、多くの場合ヒトを対象にして研究が行われる。倫理面に配慮してデータを取得する必要があることから、意図通りにデータを取得する事は必ずしも容易ではなく、寧ろ不完全データを相手にせねばならないことの方が通常である。
複数の治療法を比較する研究において、被験者に治療法を無作為に割付ける無作
為化研究は、一般的にもっとも信頼できる研究デザインとして広く受け入れられている。研究開始前に試験実施計画書を作成し、どのようなタイミングで評価項目の測定を行うかなど、詳細が事前に決定される。事前に周到に計画したとしても、副作用の発現や被験者からの申し出などの倫理的な理由から、試験途中で当該被験者の主要評
価項目のフォローアップを中断せざるを得ない中途脱落の問題が生じえる。
無作為化研究は、組み入れられた被験者内での治療法の比較の妥当性、内的妥当性
を極小の仮定で保証することが可能で、最も根拠レベルの高いゴールドスタンダードともいうべき研究デザインである。一方で、費用が高いなど実施上の困難が伴う。
対象とする被験者が限定的であり、実際の医療で対象とする症例への有効性すなわち
外的妥当性の担保には疑問が残る。このような背景から、電子カルテデータなどの既存データベースの解析により、効率的に無作為化研究を補足あるいは代替えする気運
が、リアルワールドエビデンスの名のもとに高まっている。このようなリアルワール
ドデータの解析は観察研究の解析に従うこととになるが、治療が無作為に割付けられていないことから、交絡による治療効果推定の偏りが問題となる。解析においては、交絡因子の影響を排除する必要があり、傾向スコア解析が頻繁に利用されている。
臨床試験における中途脱落に伴う欠測の問題、観察研究における交絡調整の問題などに対しては、傾向スコアによる逆確率重み付け (inverse probability weighting;IPW) 法を用いることが可能であり、実際の臨床研究での応用例も多い。通常、Missing At Random(MAR) あるいは強い意味での無視可能 (Strong Ignorabile Treatment
Assignment) 性のもとで、以下のステップで実行される。
(1) 傾向スコアにロジスティック回帰等のパラメトリックモデルを想定する。
(2) 未知パラメータを適当な推定方程式を解くことで推定する。多くの場合最尤法に基づくスコア方程式を解く。

(3) 各被験者の傾向スコアを決定し、IPW 推定を行う。
傾向スコアのモデルが正しく特定されている場合には、IPW 推定量は平均因果効果の一致推定量を与えるが、傾向スコアを誤特定している場合や、SITA 条件が成立しない場合には、交絡の影響を調整できない。傾向スコアにより調整できるのは観測されている交絡因子に限定され、調整され損ねている未測定の交絡因子の存在、残差交絡の問題が、観察研究の解析では常に問題となる。SITA 条件を前提としない解析として、アウトカム依存傾向スコアが提案され、適当なパラメトリックモデルと推定方
程式を用いることで、上述の (1)-(3) のステップで SITA 条件なしに IPW 推定を行う
ことができる (Miao et al. 2016) が、識別性の問題があり、十分に柔軟な解析ができない場合がある。未観測交絡に対する感度解析法が広く研究され、実際の臨床研究での適用例も見られるが、何らかの検証できない仮定に基づいているなど十分に使いやすいものとなっておらず、簡明かつ有効な方法の確立が望まれる。
本研究では、傾向スコアの IPW 推定に基づいて提案された、ある簡明な感度解析法の異議を検討する。アウトカム依存傾向スコアにおける (1)-(3) のステップにおいて、パラメトリックモデルの仮定 (1) を排除した上で、真の傾向スコアが取りえる範囲上での IPW 推定量の上限および下限を構成することで、平均因果効果の範囲を特定する方法を提案する。具体的には、(1) のモデルを要請しないため、(2) の推定のステップが必要がないものの、推定方程式を真の傾向スコアが満たすべき制約として利用して、線形計画法の簡単な適用で範囲を設定する。Tang et al. (2023) では、連続型あるいは二値データの場合の挙動を調べている。提案法のアイディアは簡明である
ことから、統一的に他の問題へ拡張できる可能性がある。生存時間アウトカムの場合への拡張は応用上重要である。メタアナリシスは複数の臨床試験の結果を併合して、より強固な根拠を得るための統計的方法である。公表バイアスの問題はメタアナリシ
スにおける不完全データの問題であり、特有の扱いがなされる。講演者らは、臨床試
験登録を用いることで、公表バイアスの問題が通常の不完全データと同様のデータ構造を有する見なすことができることを指摘し (Huang et al., 2021)、IPW 法での調整法を提案した (Huang et al., 2022)。本研究での感度解析法は公表バイアスの最悪評価に用いることができる。時間が許せば、これらの進展中の研究についても触れる。

引用文献
Tang C, Zhou Y, Huang A, Hattori S (2023). A simple sensitivity analysis method for unmeasured confounders via linear programming with estimating equation constraints.
arXiv:2312.04150
Miao W, Ding P, Geng Z. (2016). Identiﬁability of normal and normal mixture models
with nonignorable missing data. J Am Stat Assoc. 111, 1673–1683.
Huang A, Komukai S, Friede T, Hattori S. (2021). Using clinical trial registries to inform Copas selection model for publication bias in meta-analysis. Research Synthesis
Methods 12, 658–673.
Huang A, Morikawa K, Friede T, Hattori S. (2022). Adjusting for publication bias in metaanalysis via inverse probability weighting using clinical trial registries. Biometrics 79,
2089-2102.

野菜の価格と気象との関連性について
⽯井 ⼀夫(公⽴諏訪東京理科⼤学), ⾦打 聖菜 (公⽴諏訪東京理科⼤学)【背景と⽬的】地球温暖化は年々進⾏し私達の⽣活を脅かすものとなっている．気温上昇や⾃然災害は，農産物の収穫量や品質に⼤きく影響し，⾷料安全保障上の懸念点となっている．
野菜は，穀物などのように⽣産調整や価格調整が⾏われることが少なく，その価格に受容供給関係が直接に反映される事が多い．そこで，気候変動による経済影響を捉えるため野菜の価格と気象との関連性を統計⼿法により⽰すことを⽬的とした．
【⽅法】野菜の価格変動を捉えるため，株価の価格変動トレンドを分析するのによく⽤いられるローソク⾜チャートや移動平均線，MACD，RSI などのテクニカル分析⼿法で⽤いられる指標を⽤いることを考えた．また，気象データ（気温，降⽔量など）や農産物の作況デー
タとの経時的関連変化を，カイ⼆乗検定による p 値および-logp 値で評価した．
野菜の価格と⽣産量については，農林⽔産省の提供している⻘果物卸売市場調査（⽇別調査）[東京都]から，野菜の価格データおよび⽣産量データを⼊⼿した．気象データについて
は気象庁の提供している過去の気象データ・ダウンロードから各地の気象情報(⽇最⾼気温・
⽇最低気温・⽇平均気温・降⽔量の⽇合計・⽇照時間)を⼊⼿した．分析対象期間としては，
2010 年 10 ⽉ 1 ⽇〜2023 年 4 ⽉ 10 ⽇までの⽇々のデータを使⽤した．データ分析を⾏った．なかでも，令和 3 年産野菜⽣産統計において，収穫量ランキングの 2 位との差が 2 倍以上のもの(⽇本における特定の野菜の⽣産の多くを 1 都道府県で担っているもの)に限定
し，たまねぎ，キャベツ，レタス，トマトで分析を⾏った．
【結果】もっとも典型的な変動トレンドを⽰した北海道のタマネギの価格を⽤いた野菜の
価格変動分析の結果を図１~６に⽰した．今回，検討したテクニカル分析の指標は，MACD(図 2)，RSI(図 3)，移動平均乖離率(図 4)，ボリンジャーバンド(図 5)，ストキャス
ティクス(図 6)である．いずれの指標も価格変動に対応して，指標の変動を検出可能であった．

!" # $ 䇮

図１ ローソク⾜チャート! 俍凪$ ㏺&' ( ) *図２ ローソク⾜チャートと MACD

! 俍凪$ 咧反叒! "#図３ ローソク⾜チャートと RSI
ボリンジャーバンド変動

図 5 ボリンジャーバンド

図４ ローソク⾜チャートと移動平均線乖離率

ストキャステックス変動

図 6 ローソク⾜チャートとストキャスティクス

次に，野菜の価格変動（終値）と，野菜価格と，収穫量，降⽔量の関連性を評価した結果を
図７〜10 に⽰した．各図は上から順に，価格のローソク⾜チャート，テクニカル検定結果(MACD/RSI[1 年])，検定対象(降⽔量/⽇照時間)，カイ⼆乗検定(p 値)，カイ⼆乗検定(-logpの値)の結果である．
MACDと降⽔量でカイ⼆乗検定

図７ MACD と降⽔量の関係
RSI(1年)と降⽔量でカイ⼆乗検定

図 9 RSI(1 年)と降⽔量の関係

MACDと⽇照時間でカイ⼆乗検定

図８ MACD と⽇照時間の関係
RSI(1年)と⽇照時間でカイ⼆乗検定

図 10 RSI(1 年)と⽇照時間の関係【考察】このように，テクニカル分析指標である MACD，RSI，移動平均線乖離率，ボリ
ンジャーバンド，ストキャスティクスをもちいて野菜の価格変動を評価し，気候変動との関連について⽰すことが可能であった．

Sparse quantile regression via ℓ0 -penalty
本田敏雄
一橋大学経済学研究科
本講演では以下の内容について報告した．

1. Best subset selection
まず本研究の動機となった，情報量基準や ℓ0 ペナルティーを用いて総当たりにより変数を選択する best subset selection 手法に関して，Bertsimas et al.(2016)[1], Hastie

et al.(2020)[3], Hazimeh et al.(2023)[4] などの重要な文献について紹介した．[1] と [4]では混合整数計画法を用いている．[3] は，いくつかの手法について数値的な比較を行ったレビュー論文である．その他分位点回帰では，Chen and Lee(2023)[2] がある．

2. スパースな分位点回帰モデル， ℓ0 回帰，一致性
本研究では，（超）高次元の線形モデル，加法モデル，変動係数モデルを統一的に扱っており，特に線形モデル，加法モデルについて報告した． τ 分位点線形モデル
では以下の通りである．

Yi =

p
X

Xi j β∗j + ϵi

(1)

j=1

bは
で，非ゼロ要素の添え字集合 S0 の推定量 SeS ) ≡ Rn (βeS ) + λζn |S|Qn (βlog pnnかつb ≡ argmin Qn (β
eS ),
S

(2)

S

n
X
T
eS ≡ argmin Rn (βS ) かつ Rn (βS ) ≡ 1βρτ (Yi − XiS
βS ),
n
βS
i=1

(3)

で定義される． ρτ は τ 分位点回帰のチェック関数， pn = p ∨ n ， λ はチューニング
パラメータ， ζn は適宜指定される無限大に発散する場合もある定数である．そして本研究では以下の二つの場合を扱った． |S0 | ( ≪ pn ) は無限大に増加するとする．

1. k0 < |S0 | < k1 かつ C|S0 | < k0 がわかっており， ζn = k0 とおく場合．
2. |S0 | に関する情報はなく， ζn = 1 とおく場合．
以上のそれぞれの場合について， ℓ0 回帰モデル選択の一致性について報告した．
本講演で紹介した理論的な結果は，ある意味で Honda et al.(2019)[6] の定理 2 の真のモデルサイズが発散する場合への拡張と言え，証明手法はその論文の証明手法に新1

しいアイデアを加えた上で，より精密化したものである．詳細は Honda(2023)[5] を参照されたい．

3. 計算法[4] と [2] の手法を合わせた，近似法による加法モデルと変動係数モデルの場合の，提案した手法の計算法のアイデアについて報告した．

参考文献[1] D. Bertsimas, A. King, and R. Mazumder. Best subset selection via a modern optimization lens. Ann. Statist., 44:1–32, 2016.
[2] L.-Y. Chen and S. Lee.

Sparse quantile regression.

Journal of Econometrics,

235:2195–2217, 2023.
[3] T. Hastie, R. Tibshirani, and R. Tibshirani. Best subset, forward stepwise or lasso?
analysis and recommendations based on extensive comparisons. Statistical Science,
34:579–592, 2020.
[4] H. Hazimeh, R. Mazumder, and P. Radchenko. Grouped variable selection with discrete optimization: Computational and statistical perspectives. Ann. Statist., 51:1–32,
2023.
[5] T. Honda. Sparse quantile regression via ℓ0 -penalty. Discussion Paper Series No.
2023-03, Graduate School of Economics, Hitotsubashi University, 2023.
[6] T. Honda, C.-K. Ing, and W.-Y. Wu. Adaptively weighted group lasso for semiparametric quantile regression models. Bernoulli, 25:3311–3338, 2019.

2

対数周辺尤度とその陰的正則化のスケーリング挙動
九州大学情報基盤研究開発センター
徳田悟
古くは 17 世紀に発見されたケプラーの法則が象徴するように、単純な関数や方程式を用いて観測データを表す数理モデリングは様々な現象に対する理解を深めてきた。現象の本質を特徴付けるためにはできるだけパラメータの少ない単純なモデルが望ましい。こうした指針はオッカムの剃刀として知られる。統計学の文脈では、赤池情報量規準（AIC）やベイズ情報量規準（BIC）に象徴されるモデル選択の規準がオッカムの剃刀をサポートしている。特に BIC はベイズ統計における対数周辺尤度のラプラス近似として導出される統計量であり[1]、対数周辺尤度が有する陰的正則化と一致性を示している。しかし、この近似では観測データの質（ノイズ分散）が無視できる極限を考えており、対数周辺尤度やその陰的正則化がノイズ分散の大きさにどう依存するかは非自明な問題である。
特異学習理論は対数周辺尤度の漸近挙動の解明に貢献してきた[2]。BIC が正則なモデルに対する対数周辺尤度の漸近挙動を表す一方、特異学習理論では非正則なモデルも含めた一般化が成される。同理論では対数周辺尤度の有する陰的正則化が実対数閾値（RLCT）と呼ばれる双有理不変量によって定量化される。しかし、特異学習理論もまたノイズ分散が無視できる極限を考えており、前述の問題に対する回答を提示しない。
対数周辺尤度やその陰的正則化がノイズ分散に依存することを示すスケーリング挙動を明らかにする上で、以下の関連問題が挙げられる。
⚫

RLCT は陰的正則化の極限値であり、観測可能な統計量でない⚫

特異学習理論ではノイズ分散𝛽 −1 = 𝑂(1)の下でサンプルサイズの極限𝑛 → ∞に着目⚫

対数周辺尤度は𝛽 = 𝑂(1)の下で一致性を示すが、𝑛𝛽 = 𝑂(1)での挙動は不明

これらの問題に対し、我々は次のように貢献した[3]（図 1）。
✓

RLCT の観測可能な対応物（ベイズ比熱）を導入✓

𝑛𝛽 = 𝑂(1)の下での極限𝑛 → ∞に着目し、𝑛𝛽に依存するスケーリング関数を導出✓

スケーリング関数を元に対数周辺尤度の𝑛𝛽への依存性や陰的正則化の寄与を実証

特に、3 つ目の貢献を通じて、ベイズ推定が観測データの質や量に応じた複数の「状態」を取り、状態によって異なるモデルを最適とみなす性質を発見した。観測の精細度𝑛𝛽が高い（データの量が多いまたは質が高い）状態であるほど、より多くのパラメータを持つ複雑なモデルを最適とみなす性質である。こうしたスケーリング挙動は BIC では説明がつかず、導出したスケーリング関数によって初めて明らかになったものである。
以上の知見は研究者の洞察に頼ってきたモデルの単純化をデータに根ざして客観化し、様々な現象の理解に役立つことが期待される。モデルの妥当性を実証するために「どの程の量や質のデータが必要か？」に答えるものでもあり、計測の効率化にも示唆を与える。

図 1：先行研究に対する我々の貢献の位置付け
参考文献：
1.

G. Schwarz, “Estimating the dimension of a model”, Ann. Statist. 6(2), 461-464 (1978).

2.

S. Watanabe, Mathematical Theory of Bayesian Statistics (CRC Press, Boca Raton, 2018).

3.

S. Tokuda, K. Nagata, and M. Okada, “Intrinsic regularization effect in Bayesian nonlinear
regression scaled by observed data”, Phys. Rev. Research 4, 043165 (2022).

所得不平等度の分解の推定
千葉大学

川久保 友超名古屋市立大学

各務 和彦

一般化エントロピーとその分解可能性

1xi (i = 1, . . . , N ) を，母集団における i 番目の世帯の所得とし，添字集合 I = {1, . . . , N } によって有限母集団を表すとする．ここで，I の所得不平等の尺度として，以下で定義される一般化エントロピー尺度（generalized entropy measures, GE）を考える．

1
1 X
GEθ =
θ(θ − 1) N

(

i∈I

ただし，µ = N −1

xi
µ

)

θ

−1 ,

θ ̸= 0, 1.

(1)

PN

i=1 xi は母集団 I における平均所得を表す．

次に，母集団 I （全国）を互いに排反な部分母集団（地域）I1 , . . . , IJ に分解し，Ij の GE を考える．ここで，Nj = #(Ij ) を地域 Ij のサイズ（世帯数），λj = Nj /N を Ij の世帯数シェア，µj = Nj−1

Pi∈Ij xi を Ijの平均所得とする．このとき，地域 Ij の GE は，

1
1 X
GEθ,j =
θ(θ − 1) Nj

(

i∈Ij

xi
µj

)

θ

−1 ,

θ ̸= 0, 1.

と定義される．

GE には，以下の分解可能性という望ましい性質がある．
GEθ =

J
X

λ1−θ
sθj GEθ,j +
j

j=1

{z

|

地域内不平等

}


J
X

1
λj
θ(θ − 1)  j=1
{z
|



µj
µ

θ



−1 .

}

地域間不平等

ただし，sj = λj µj /µ は Ij の所得シェアを表す．分解の第一項は地域内不平等（within-region inequality）を表し，各地域の GE の加重平均と理解できる．一方で第二項は地域間不平等（between-region inequality）を表し，「各地域の全世帯の所得が，その地域の平均所得 µj である場合の母集団の GE」と理解できる．

2

GE の分解の推定2.1

グループデータにもとづいた GE の推定

公的統計における所得データの多くは，ある標本調査によって収集された所得の情報が，グループデータの形式で公開されている．S ⊆ I は，ある標本調査において標本として選ばれた世帯 i をあつめた添字集合とする．ここで，予め決められた境界値 0 = c0 < c1 < · · · < cG = +∞ に対して，それぞれの所得階級内の標本における世帯数yg =X

I(cg−1 ≤ xi < cg ),(g = 1, . . . , G),i∈S

がカウントされ，度数分布 y = (y1 , . . . , yG )⊤ が公表されているとする．

グループ数 G = 10 程度のグループデータ y にもとづいて GE を推定するために，パラメトリックな所得分
布を仮定する．F (· | ω) および f (· | ω) を，パラメータ ω を持つ確率分布の累積分布関数および確率密度関数とする．このとき，グループデータ y にもとづいた尤度は，L(ω) ∝GYy{F (cg | ω) − F (cg−1 | ω)} g .

g=1という多項分布の尤度となる．密度関数 f (· | ω) を持つ確率分布の GE は，(1) 式の有限母集団の GE の定義における有限母集団平均（算術平均）をとる操作を，確率分布の平均（期待値）におきかえた以下の式で定義される．

1
GEθ (ω) =θ(θ − 1)(Z  )θ∞x
f (x | ω)dx − 1 ,µ0

所得分布として当てはまりの良いとされる多くのパラメトリックな確率分布は，GEθ (ω)，T0 (ω) および T1 (ω)
が ω の陽な関数として書ける．よって，尤度 L(ω) にもとづいて ω が推定できれば，GE も容易に推定できる．

2.2

GE の分解の推定

次に，全国の GE と各地域の GE を同時に推定する問題を考える．各地域の標本添字集合を Sj = S ∩ Ij と定義し，各地域の世帯所得の度数分布 yR,j = (yj1 , . . . , yjG )⊤ が観測されているとする．ただし，

yjg =

X

I(cg−1 ≤ xi < cg ),(g = 1, . . . , G).

i∈Sjである．y と yj (j = 1, . . . , J) が観測されているもとで，以下の推定手法を提案する．
まず，全国のパラメトリック分布 f (· | ω) と事前分布 ω を仮定して全国の GE を推定し，その推定値をdGEθ とする．サンプルサイズとグループ数 G に応じて適切なパラメトリック分布を選択することで，全国のGE を効率的に推定することができる．次に，各地域 Ij のパラメトリック分布 fj (· | ωj ) と事前分布 πj (ωj )PJを仮定し，Ij の平均所得推定値 µ̂j を求める。また，µ̂ = j=1 λj µ̂j として全国の平均所得推定値 µ̂ を求める．なお，世帯数シェア λj は既知とする．µ̂j と µ̂ にもとづいて，地域間不平等を

b=
B


J
X

1
λj
θ(θ − 1)  j=1



µ̂j
µ̂

θ



−1 .


(2)

と推定する．最後に，分解の等式を制約として，各地域の GE を制約付きベイズ推定量によって推定する．制
約付きベイズ推定量は，以下の事後リスクの制約付き最小化問題

mind1 ,...,dJの解であり，

JX

hi
2
ϕj E {dj − GEθ,j (ωj )} | yRdθ =s.t. GE

j=1



CB
B
dθ − Bb − GEw ,d = GEd + rj GE
GE
θ
θ,j
θ,j
q

J
X

wj dj + B̂,

(3)

j=1

j = 1, . . . , J,

と解ける．

3

日本の所得データを用いた実証分析
住宅・土地統計調査の所得データを用いて，全国の GE を都道府県間不平等と都道府県内不平等に分解し，さらに各都道府県の GE を市区町村間不平等と市区町村内不平等に分解する実証分析を行った．

Bayesian analysis of stochastic conditional duration models withintraday and intra-deferred future seasonalities in high-frequency
commodity market
Makoto Nakakita1 , Tomoki Toyabe2 , and Teruo Nakatsuma2
1

2

AIP, RIKEN, Japan
Faculty of Economics, Keio University, Japan

1 Introduction
We propose an extension of the stochastic conditional duration (SCD) model to capture intraday
seasonality patterns in stock trading intervals and an improved Markov chain Monte Carlo (MCMC)
method. The trading time intervals of financial data are known to have very slow-decaying autocorrelations, and the Ancillarity-Suﬃciency Interweaving Strategy (ASIS) overcomes the disadvantage of
instability in estimation by incorporating a centered parameterization of the model. The estimation
results using gold futures tick data from JPX show that trading time intervals for each forward month ofthe contract month has a trend of a certain pattern and that the use of ASIS leads to stable estimation
results.

2 Stochastic Conditional Duration Model
Considering the time of the ith transaction τi and the duration yi = τi+1 − τi as a product of a latentvaluable αi and a positive random variable εi , the SCD model with intraday seasonality is written as:yi = exp (x′i β + αi ) εi ,αi = φαi−1 + ηi ,2εi > 0,2

ηi ∼ N (0, σ ),2α1 ∼ N (0, σ /(1 − φ )),L
xi = (xBi ; xi )B|φ| < 1,i ∈ {1, . . . , N },(1)i ∈ {2, . . . , N },(2)L

β = (β ; β ).

Note that the transaction time is normalized to [0, 1] (i.e., the market starts at time 0 and ends attime 1) and both εi and ηi are mutually and serially independent in this model. As for the distributionof the random variable εi , the weibull distribution and the gamma distribution have been used in theliterature. Following Toyabe and Nakatsuma (2022), we use Gamma(γ, γ) instead of the commonly usedGamma(γ, 1) for computational stability. Thus the probability density function (pdf) of εi and yi are⎧ γ−1⎨γεi exp(−εγi ), (Weibull distribution),γ γ γ−1p(εi |γ) =(3)⎩εexp(−γεi ), (gamma distribution),Γ(γ) i
and

⎧
γ−1
γ
′
′
⎪
⎪γyi exp (−γ(xi β + αi ) − yi exp(−γ(xi β + αi ))) ,⎪⎪⎪⎨ (Weibull distribution),γ γ γ−1p(yi |αi , β, γ) =⎪yexp (−γ (x′i β + αi + yi exp(−x′i β − αi ))) ,⎪⎪ Γ(γ) i⎪⎪⎩ (gamma distribution).

(4)

To capture the structure of the intraday seasonality of trading time intervals, this model considers thatB
xB′
is a Bernstein polynomial
i β
B
xB′
i β =

N
%

βkB xk,i =

k=0

& '
m k
bk,m (v) =v (1 − v)m−k ,kN%βkB bk,m (τi+1 )

(5)

k=0

k = 0, . . . , m,v ∈ [0, 1].

(6)Note that the τi+1 , not the time interval yi , is in the equation here. This allows us to estimate thepattern of intraday seasonality with time as an argument.
By the same method as intraday seasonality, we added the order book information asxL′
i = (di1 , di2 , di3 , di4 , di5 )(7)

where di1 is the price return, di2 is the absolute value of price return, di3 is the trade volume, di4 = di1 ∗di3
and di5 = di2 ∗ di3 .

3 Empirical results
Setting of hyperparameters is as follows.
(aφ , bφ ) = (3, 2),(aτ , bτ ) = (3, 10),(γa , γb ) = (3, 10)

In the proposed model, it is not feasible to estimate the degree of the Bernstein polynomial using
MCMC. We use the widely applicable information criterion (WAIC) proposed by Watanabe (2013) as acriterion for model selection. The detailed results are omitted here for reasons of space limitation.

4 Conclusions
From the curve depicting the pattern and seasonality of trading time intervals over a two-month period
estimated from the proposed model, it can be observed that, consistently, trading time intervals exhibit aninverted U-shaped form. This represents an approach that captures the continuous behavior of the mostactively traded futures, previously considered merely for the sake of removing strong autocorrelation. Italigns with the results of previous research Toyabe and Nakatsuma (2022), showing an inverted U-shaped
intraday seasonality in stock trading.
2微小ノイズをもつ空間 2 次元線形放物型確率偏微分方程式
モデルの係数パラメータの推定
大阪大学大学院基礎工学研究科

渡名喜庸蔵

1 モデルと設定
D = (0, 1)2 上の線形放物型確率偏微分方程式
dXt = (θ2 ∆ + (θ1 , η1 ) · ∇ + θ0 )Xt dt + ϵdWt ,t ∈ [0, 1]を考える. ただし, 初期条件 X0 ∈ L2 (D), ディリクレ境界条件 Xt (y, z) = 0, (t, y, z) ∈ [0, 1]×∂D を満たすとする. ここで {Wt }t≥0 は D 上のソボレフ空間上の Q-ウィーナー過程とし, (θ0 , θ1 , η1 , θ2 ) ∈ R3 ×

(0, ∞) を未知パラメータ, パラメータ空間 Θ ⊂ R3 × (0, ∞) はコンパクト凸集合で, (θ0∗ , θ1∗ , η1∗ , θ2∗ ) ∈Int (Θ) を真値とする. また ϵ ∈ (0, 1] を既知の微小撹乱パラメータとする.
−Aθ = θ2 ∆ + (θ1 , η1 ) · ∇ + θ0 とおく. l1 , l2 ∈ N, (y, z) ∈ D に対して
λl1 ,l2 = π 2 θ2 (l12 + l22 ) +θ12 + η12
− θ0 ,
4θ2u, v ∈ L2 (D) に対して ⟨u, v⟩ =

R1R1
0

0

θ1

η1

el1 ,l2 (y, z) = 2 sin(πl1 y) sin(πl2 z)e− 2θ2 y e− 2θ2 z ,θ1η1u(y, z)v(y, z)e θ2 y e θ2 z dydz, ∥u∥ =p⟨u, u⟩ とおく. α ∈ (0, 1)は既知で, µl1 ,l2 = π 2 (l12 + l22 ) + µ0 とする. ただし µ0 ∈ (−2π 2 , ∞) は未知パラメータでパラメータ空間は (−2π 2 , ∞) のコンパクト凸集合で, 真値 µ∗0 はその内部に属するとする. Q-ウィーナー過程を次で
定める.

Wt =

X

−α/2

µl1 ,l2 wl1 ,l2 (t)el1 ,l2 ,t ≥ 0.

l1 ,l2 ≥1

ただし {wl1 ,l2 (·)}l1 ,l2 ∈N は独立な 1 次元標準ウィーナー過程である.
観測データとして高頻度時空間データ XN,M = {Xti (yj , zk )}0≤i≤N,0≤j≤M1 ,0≤k≤M2 を考える. た
だし ti = i∆ = i/N , yj = j/M1 , zk = k/M2 , M = M1 M2 とする. また, 初期値 X0 ∈ L2 (D) は⟨X0 , e1,1 ⟩ ̸= 0, ∥Aθ X0 ∥ < ∞ を満たし, λ∗1,1 = 2π 2 θ2∗ +(θ1∗ )2 +(η1∗ )2− θ0∗ > 0 とする.
4θ2∗

2 θ1 , η1 , θ2 の推定
b ∈ (0, 1/2) を固定し, b ≤ ye0 < ye1 < · · · < yem1 ≤ 1 − b, b ≤ ze0 < ze1 < · · · < zem2 ≤ 1 − b を満たす
(1)

観測データ XN,M の空間間引きデータ XN,m = {Xti (eyj , zek )}0≤i≤N,0≤j≤m1 ,0≤k≤m2 を考える.

∆i X(y, z) = Xti (y, z) − Xti−1 (y, z) に対して
ZN (y, z) =N
1 X
∆Xi2 (y, z),N ∆α i=1fα (y, z : θ1 , η1 , θ2 ) =Γ(1 − α) − θθ1 y − ηθ1 z
e 2 e 24παθ21−α(1)(Γ(·) はガンマ関数) とおき, XN,m を用いて (θ1 , η1 , θ2 ) の最小コントラスト推定量を次で定める.
(θ̂1 , η̂1 , θ̂2 ) = arginf
θ1 ,η1 ,θ2

m1 n
m2 X
X
k=1 j=1

ϵ−2 ZN (eyj , zek ) − fα (eyj , zek : θ1 , η1 , θ2 )

o2
.

2

2−α

N
RN,ϵ = N 1/2∧(1−α) ∧ ϵN 5/4−α ∧ ϵ log
N とすると N, m1 , m2 → ∞, ϵ → 0 のとき正則条件の下,RN,ϵ (θ̂1 − θ1∗ , η̂1 − η1∗ , θ̂2 − θ2∗ ) = OP (1)が成り立つ.

3 θ0 の推定
in ≤ N とする. eti = ⌊ Nn ⌋ N , i = 0, . . . , n に対して観測データ XN,M の時間間引きデータ(2)

Xn,M = {Xeti (yj , zk )}0≤i≤n,0≤j≤M1 ,0≤k≤M2 を考える.
x1,1 (t) := ⟨Xt , e1,1 ⟩ は以下を満たす 1 次元オルンシュタイン・ウーレンベック過程となる.
−α/2

dx1,1 (t) = −λ1,1 x1,1 (t)dt + ϵµ1,1 dw1,1 (t),x1,1 (0) = ⟨X0 , e1,1 ⟩.

(2)

Xn,M と上で構成した推定量 θ̂1 , η̂1 , θ̂2 を用いた x1,1 (·) の近似過程
2xbl1 ,l2 (eti ) =
M

M2 X
M1
X

θ̂1

λ,µ

η̂1

zk

k=1 j=1

による (λ1,1 , µ1,1 ) の最小コントラスト推定量(λ̂1,1 , µ̂1,1 ) = arginf

yj

Xeti (yj , zk ) sin(πl1 yj ) sin(πl2 zk )e 2θ̂2 e 2θ̂2

( n
X (b
x1,1 (eti ) − e−λ/n xb1,1 (eti−1 ))2ϵ2 (1−e−2λ/n )
2µα λ

i=1

)
1 − e−2λ/n+ n log.
2µα λ/nを考え, θ0 , µ0 の推定量をθ̂0 = −λ̂1,1 +−2λで定める. G1 (λ, µ) = 1−e2λ

ϵ−1 (θ̂0 − θ0∗ ),θ̂12 + η̂124θ̂2+ 2π 2 θ̂2 ,µ̂0 = µ̂1,1 − 2π 22αµα x1,1 (0)2 , G2 (µ) = 2µ2 とすると n → ∞, ϵ → 0 のとき正則条件の下,

√

 d

n(µ̂0 − µ∗0 ) → N 0, diag {G1 (λ∗1,1 , µ∗1,1 ), G2 (µ∗1,1 )}−1が成り立つ.

4 数値シミュレーション
高頻度時空間データ XN,M , N = 103 , M1 = M2 = 200 による大規模数値シミュレーションを実行
し, 推定量 θ̂1 , η̂1 , θ̂2 , θ̂0 , µ̂0 の漸近的挙動を検証する.

参考文献
M. Bibinger and M. Trabs. (2020). Volatility estimation for stochastic PDEs using high-frequency observations. Stochastic Processes and their Applications, 130(5), 3005–3052.
Y. Kaino and M. Uchida. (2021). Adaptive estimator for a parabolic linear SPDE with a small noise.
Japanese Journal of Statistics and Data Science, 4, 513–541.
Y. Tonaki, Y. Kaino, and M. Uchida. (2023). Parameter estimation for linear parabolic SPDEs in two spacedimensions based on high frequency data. Scandinavian Journal of Statistics, 50(4), 1568–1589.
Y. Tonaki, Y. Kaino, and M. Uchida. (2023). Parameter estimation for a linear parabolic SPDE in two
space dimensions with a small noise. Statistical Inference for Stochastic Processes. https://doi.
org/10.1007/s11203-023-09301-2.
M. Sørensen and M. Uchida. (2003). Small-diﬀusion asymptotics for discretely sampled stochastic diﬀerential equations. Bernoulli, 9(6), 1051–1069.

VAR models
for surface time series
fPCA analysis of nighttime population in Tokyo
during Jan. 2016-Dec. 2023
Yasumasa Matsuda and Daisuke Kurisu
Tohoku University and University of Tokyo

1

Data(東京圏夜間人口の分析)• 東京圏2400メッシュ(1kmx1km)における午前4時の人口• 2016/1-2018/12の月平均を収集• 96か月 x 2400メッシュの時空間データ

右図の1メッシュを
10x10に分けたメッシュを用いる

2

PC function for 対数差分系列
都心回帰

郊外化

立川
江戸川区

都心回避

東京駅

流山

3

主成分スコア for 原系列
1st PC score

2nd PC score

3rd PC score
2020/1月
コロナによる影響がみられるのは、第３主成分の都心回避成分。現在(2023/12)では元の水準に回復
郊外化、都心回帰成分はコロナの前後でもほぼ変わりなく、都心回帰成分が増加傾向にある

4