高次元複雑データの統計モデリング報告書
A 15H01678 大：

複雑データのの大学：
日時：8 月 29 日 〜 8 月 30 日
場

：九州大学伊都キャンパス ウエスト 1 号館 C501 大講義室
：

年

高次元

データ

複雑

データ
複雑

計
の

リ
の

ン

の

統計モデリング
高次元複雑データの統計モデリングのウの高次元複雑データの統計モデリングの統計

科研費シンポジウム「高次元複雑データの統計モデリング」
講演プログラム

8月29日（木）
時間

名前

所属

10：00〜10：20

タイトル

オープニング (青嶋誠，廣瀬慧)

10：20〜11：00 矢田 和善 筑波大学
石井青嶋
11：00〜11：40

晶誠

東京理科大学
筑波大学

新村 秀一 成蹊大学

11：40〜13：00

高次元遺伝子解析の呪いからの解放1 -統

廣瀬 慧

計が1970年からこの問題を解決できなか
った理由-

朴 喜楥

佐賀大学

脳画像における高次元データ解析法

広島大学

Sparse Tensor subspace method forbiological modulators selection

14：20〜14：40
14：40〜15：20

Geometrical quadratic discriminant
analysis for high-dimension, strongly
spiked eigenvalue models

休憩

13：00〜13：40 川口 淳
田尻 涼
13：40〜14：20

座長

二宮 嘉行

休憩
二宮 嘉行 統計数理研究所

因果推論におけるセミパラメトリックアプ
ローチのための情報量規準

15：20〜16：00 作村 建紀 法政大学

方向データに対する von Mises 分布の母

16：00〜16：30

A Practicable Estimation of Mean Squared
Prediction Error in Small Area Estimation

柳本 武美 統計数理研究所 数推定について
廣瀬 雅代 九州大学

16：30〜17：00

総合討論

18：30〜20：30

懇親会@大名つつじ庵

矢田 和善

科研費シンポジウム「高次元複雑データの統計モデリング」
講演プログラム

8月30日（金）
名前

所属

タイトル

10：20〜11：00 野村 俊一

統計数理研究所 前震の統計的判別に基づく地震予測

11：00〜11：40 鈴木 拓海

東京大学

11：40〜13：00

Penalized least squares approximation
methods

高井 啓二

休憩13：00〜13：40 高井 啓二

関西大学

欠測データを用いたフィッシャースコア
リング法

13：40〜14：20 森川 耕輔

大阪大学

Bayesian Fractional Imputation With
Nonignorable Nonresponse Data

14：20〜14：40

廣瀬 雅代

休憩

14：40〜15：20 和田 裕一郎 名古屋大学

スペクトラル埋め込みを利用した深層ク

15：20〜16：00 小林 景

距離空間の変形と埋め込みを用いたデー

16：00〜16：20

座長

慶應義塾大学

ラスタリングタ解析

クロージング増田 弘毅

Geometrical quadratic discriminant analysis for high-dimension,strongly spiked eigenvalue models
Kazuyoshi Yataa , Aki Ishiib and Makoto Aoshimaaab

1Institute of Mathematics, University of Tsukuba
Department of Information Sciences, Tokyo University of Science

Introduction

We consider a classifier for high-dimensional data under the strongly spiked eigenvalue (SSE) model. We
create a new classification procedure on the basis of the high-dimensional eigenstructure. We propose aquadratic classification procedure by using a data transformation. Suppose we have two classes πi , i = 1, 2,
and define independent p × ni data matrices, X i = [xi1 , ..., xini ], i = 1, 2, from πi , i = 1, 2, where
xij , j = 1, ..., ni , are independent and identically distributed (i.i.d.) as a p-dimensional distribution witha mean vector µi and covariance matrix Σi (≥ O). We assume that lim supp→∞ ||µi ||2 /p < ∞ fori = 1, 2, where || · || denotes the Euclidean norm. Also, we assume that tr(Σi )/p ∈ (0, ∞) as p → ∞for i = 1, 2. Here, for a function, f (·), “f (p) ∈ (0, ∞) as p → ∞” implies lim inf p→∞ f (p) > 0
and lim supp→∞ f (p) < ∞. We assume ni ≥ 3, i = 1, 2. The eigen-decomposition of Σi is given by∑
Σi = H i Λi H Ti = ps=1 λs(i) hs(i) hTs(i) , where Λi = diag(λ1(i) , ..., λp(i) ) having λ1(i) ≥ · · · ≥ λp(i) (≥ 0)
and H i = [h1(i) , ..., hp(i) ] is an orthogonal matrix of the corresponding eigenvectors. Let X i −[µi , ..., µi ] =1/2

H i Λi Z i for i = 1, 2. Then, Z i is a p × ni sphered data matrix from a distribution with the zero mean
and identity covariance matrix. Let Z i = [z 1(i) , ..., z p(i) ]T and z j(i) = (zj1(i) , ..., zjni (i) )T , j = 1, ..., p, fori = 1, 2. We assume that the fourth moments of each variable in Z i are uniformly∑ni bounded for i = 1, 2. Letz oj(i) = z j(i) −(z̄j(i) , ..., z̄j(i) )T , j = 1, ..., p; i = 1, 2, where z̄j(i) = n−1i
k=1 zjk(i) . We also assume that()
P lim inf p→∞ ||z o1(i) || ̸= 0 = 1 for i = 1, 2. Let x0 be an observation vector of an individual belonging∑ ixij /nito πi (i = 1, 2). We assume x0 and xij s are independent. We estimate µi and Σi by x̄ini = nj=1
∑ i(xij − x̄ini )(xij − x̄ini )T /(ni −1). As for the HDLSS data, Aoshima and Yata [1] considand S ini = nj=1
ered distance-based classifiers and Aoshima and Yata [1] gave the misclassification rate adjusted classifier
for multiclass, high-dimensional data whose misclassification rates are no more than specified thresholds
under the following condition for eigenvalues:
λ21(i)tr(Σ2i )→ 0 as d → ∞ for i = 1, 2.

(1.1)

Aoshima and Yata [3] called (1.1) the “non-strongly spiked eigenvalue (NSSE) model”. Also, Aoshima andYata [3] gave the following eigenvalue model called the “strongly spiked eigenvalue (SSE) model”:{ λ2
lim infd→∞1(i)tr(Σ2i )}> 0 for i=1 and 2.

(1.2)

Aoshima and Yata [4] considered a liner classifier under the SSE model 1.2 by using a data transformationfrom the SSE model to the NSSE model. They gave a consistency property of the classifier and discussed
the asymptotic normality when p → ∞ and ni → ∞ (i = 1, 2). On the other hand, Ishii [5] proposed alinear classification procedure which has the consistency property even when ni ’s are very small under thefollowing eigenvalue condition:∑p2
s=2 λs(i)= o(1) as p → ∞ for i = 1, 2.
(1.3)
λ21(i)
Note that (1.3) implies the conditions that λ2(i) /λ1(i) → 0 and λ21(i) /tr(Σ2i ) → 1 as p → ∞.

2

New geometrical quadratic discriminant analysis under the SSE model

Aoshima and Yata ([2]) gave a geometrical quadratic discriminant analysis (GQDA) under the NSSE model(1.1). We consider a new GQDA for (1.3). We assume the following condition:(A-i) h1(1) = h1(2) (= h1 , say) and λ1(1) /λ1(2) ∈ (0, ∞) as p → ∞.
Note that (A-i) is much milder than Σ1 = Σ2 . Aoshima and Yata [4] considered a distance-based classifierby using a data transformation from the SSE model to the NSSE model. They gave the consistency property
for the classifier and discussed the asymptotic normality when p → ∞ and ni → ∞ (i = 1, 2). On the other
hand, Ishii [5] gave a distance-based classifier by using the data transformation when p → ∞ while ni ’s are
fixed. In this talk, we create a new quadratic classifier by using the data transformation. We construct thefollowing new GQDA:e DT (x0 ) =pG||x0 − x̄1n1 ||2 − {(x0 − x̄1n1 )T h̃1(2) }2−ptr(S 1n1 ) − λ̃1(1)||x0 − x̄2n2 ||2 − {(x0 − x̄2n2 )T h̃1(1) }2tr(S 2n2 ) − λ̃1(2){}tr(S 2n2 ) − λ̃1(2)pp−+− p log.
n1 n2
tr(S 1n1 ) − λ̃1(1)(2.1)

e DT (x0 ) < 0 and π2 otherwise. Here, λ̃1(i) ’s and h̃1(i) ’s are the noise
Then, one classifies x0 into π1 if G
reduction (NR) estimators of λ1(i) ’s and h1(i) ’s.
Theorem 2.1. Assume (A-i), (1.3) and some regularity conditions. For the classifier given by (2.1), we havethat as p → ∞e(1) → 0 and e(2) → 0.
(2.2)
Here, e(i) denotes the error of misclassifying an individual from πi into πj for i, j = 1, 2 and i ̸= j.
e DT (x0 ) by simulation studies and real data examples.
We gave the performances of G

References[1] Aoshima, M., Yata, K.: A distance-based, misclassification rate adjusted classifier for multiclass, highdimensional data. Ann. Inst. Statist. Math. 66, 983-1010 (2014).
[2] Aoshima, M., Yata, K.: Geometric classifier for multiclass, high-dimensional data. Sequential Anal.
(Special Issue: Celebrating Seventy Years of Charles Stein’s 1945 Seminal Paper on Two-Stage Sampling) 34, 279-294 (2015).
[3] Aoshima, M., Yata, K., 2018. Two-sample tests for high-dimension, strongly spiked eigenvalue models.
Stat. Sin. 28, 43-62.
[4] Aoshima, M., Yata, K.: Distance-based classifier by data transformation for high-dimension, strongly
spiked eigenvalue models. Ann. Inst. Statist. Math. 71, 473-503 (2019).
[5] Ishii, A.: A classifier under the strongly spiked eigenvalue model in high-dimension, low-sample-size
context. Commun. Stat. Theory Methods. (2019), in press.

高次元遺伝子解析の呪いからの解放
成蹊大学名誉教授

新村秀一「高次元 Microarray データを用いて癌遺伝子の特定と癌の亜種を見つける研究」が 1970 年頃から行われてきた(Golub, 1999)。これらの研究で用いられたデータが公開されているので、統計に限らず機械学習（AI）、パター
ン認識、Bio 工学の新テーマとして研究されてきたが、いずれの研究も成功していない(判別分析の Problem5)。
しかし、仮に症例数が n=100 で遺伝子数が p=10,000 の発現量データとすれば、2 群判別が最も適した手法である。筆者は 2015 年 10 月 28 日から 12 月 20 日の僅か 54 日間で、簡単にこの問題を解決した。用いたデータは、1999 年から 2004 年の間に米国の 6 医学研究プロジェクトが論文を発表し、研究に用いた公開データである。これらは癌と健常、あるいは 4 種の異なった癌の 2 クラスである。
結果は非常に単純である。6 種のデータは線形分離可能なデータ(Linearly Separable Data, LSD)である（判別分
析の Fact3）。この重要な事実であり信号が、これまでの研究で誰も指摘していない。さらに、筆者が開発したMatryoshka Feature Selection Method (Method2)で簡単に線形分離可能な含まれる遺伝子数が n 個以下の遺伝子のｋ組の部分空間(SMall Matryoshka, SM)と最小誤分類数(Minimum Number of MisclassificationSMNM)が 1 以上の遺伝子の雑音部分空間に分割できた(Fact4)。各 SM は統計分析が容易な小標本であるが、ロジスティック回帰以外の統計手法（一元配置の分散分析、ｔ検定、相関分析、クラスター分析、PCA）で線形分離可能な事実が示されなかった(Problem6)。そこで、MNM 基準による改定 IP-OLDF（RIP）の判別スコア(RIP Discriminant Score, RipDS)を変数とし、n 症例*k 次元の信号データ(k<=n)を作成した。これを上記の統計手法で分析し「癌の遺伝子診断の統計分析法を世界で初めて提案」できた。今年 5 月に Springer から Springer2 を出版した。
以上の研究が簡単にできたのは、大学卒業以来行ってきた判別分析の新理論(Springer1）が 2015 年に完成し、新理論がその応用問題として 1970 年頃から未解決の Microarray を用いた癌の遺伝子解析（Problem5）を簡単に解決できた。本来であれば癌の遺伝子研究の専門家でない筆者が「癌の遺伝子診断」までを行うことは適していない。しかし、癌は遺伝子の病気であり、高次元の Microarray 空間で 2 群が完全に分かれていて、さらに MNM＝０である k 組の SM に分割できる。そして RipDS で信号データを作ることで、上記の統計手法で線形分離可能であり有効と考えられる Malignancy Indexes が数多く発見できた。しかし、これ等のどの Malignancy Indexes が医学的に役に立つか否かは医学専門家の検証が必要である。
残念ながら Golub らの研究後に、「NIH が乳がん以外の癌に関して Microarray による研究は成果が出ないと判断し、医学研究が終わったようである」。このため、いかに医学専門家の検証につなげるかを 2016 年から模索している。しかし統計や工学研究者は、NIH の報告を知らずに研究を続けているのは一般的に問題であろう。また、データが LSD であるのに、そのデータを学習標本に用いた AI 研究が LSD の事実を指摘しない点だけが、まだ説明できていない。
大学卒業以来の研究テーマである判別分析の新理論を確立し、その応用として「高次元 Microarray の癌の遺伝
子解析と診断」にはじめて成功した。そこでこれまでの研究を見直した結果、LSD である高次元データは、ケー
ス数 n 個以下の遺伝子のｋ組の小標本に必ず分割できるという事実が統計にとって一番重要と考えた。
すなわち、
我々は「高次元データの呪いから数理計画法(MP)の LP と IP で定式化した LDF で解放される（2 次計画法 QPで定式化した SVM ではできない）。そして、分割された SM を統計分析するとこれまで見えてこなかった新しい癌の遺伝子診断の世界が広がる」。
以上の重要なテーマを以下の 4 回のシンポジュームで報告したい。
1) 九州大学：「高次元遺伝子解析の呪いからの解放 1 -統計が 1970 年からこの問題を解決できなかった理由-」
12) 新潟大学：「高次元遺伝子解析の呪いからの解放 2 -癌の遺伝子診断-」3) 東京工業大学：「高次元遺伝子解析の呪いからの解放 3 -機械学習などの工学研究の問題点-」
4) 秋田大学：「高次元遺伝子解析の呪いからの解放 4 -高次元データの分割法の最新結果-」
予稿集は、以下の 6 章からなり、1 章と 2 章は共通の基礎知識、3 章から 6 章は 4 大学での発表に対応していて、各大学の予稿はその章しか含んでいないで、必要であれば他大学の予稿を参考にしてほしい。
１章では、Springer1 と「新村(2010). 最適線形判別関数．日科議連」の中から、統計的判別関数がなぜ癌の遺伝子解析に役に立たなかったかの理由を報告する。すなわち LSD 判別は、MNM 基準による改定 IP‐OLDF(RIP)とハードマージン最大化 SVM（H-SVM）でしか理論的に正しくできないことが原因である。回帰係数や判別係
数を０にすることで、Problem5 に対応できると考える LASSO 研究の間違いを指摘する
２章では、1970 年頃から解決できなかった Problem5 の結果を紹介する。
３章では、高次元の Microarray がなぜ n 個以下の MNM＝０であるｋ組の SM に簡単に分割できるかを数理計画法(MP)の基礎知識（新村(2010).

数理計画法による問題解決法．日科議連」）と連立方程式の解の基礎知識の簡単な組み合わせで説明する。そして、統計的判別関数と 2 次計画法 QP で定式化された H

SVM が高次元のLSD である Microarray をｋ組の SM に分割できない理由を説明する。
4 章では、本研究における 8 種類の LDF の役割を概観する。Microarray が LSD である Fact3 は、RIP、H-SVM、Revised LP-OLDF と SVM4(C=104)で発見できる。そして RIP と Revised LP-OLDF が高次元の呪いから研究者を解放し、SVM ができない理由を示す。高次元 Microarray はｋ組の n 個以下の遺伝子の部分空間の SM に分割できる。これらは小標本であり統計手法で簡単に分析し癌の遺伝子診断が行えると考えた。しかし、ロジスティック
回帰だけが全ての SM がＮＭ＝０で、LSD であることが分かる。しかし他の統計手法で LSD の事実が得られなかった（Problem6）。試行錯誤の末、信号データを作成してこれを解決し、癌の遺伝子診断が可能になった。しかし LSD である Microarray データを学習に用いているのに、なぜ AI 研究は LSD の事実を発見できないかを参加者と議論したい。6 月開催の IEEE の機械学習の国際会議で良い情報が得られればそれも報告する。
5 章では、世界で初めて成功した癌の遺伝子診断の結果の概略を Springer２から説明する。
6 章では、Method２が Microarray データだけでなく、6 変数の普通車と小型車の 2 群判別にも適用できることを示す。即ち本研究は、LSD は必ずより小さい SM や最小次元の Basic Gene Sets(BGS)に分割できることが今後の統計にとって重要なことを示す。そして Method2 が求めた SM をさらに BGS（iPS 研究の山中 4 因子と同じ概念）に分割すると、多くの場合に 2 個の BGS が含まれることが分かった。
筆者の方法によって、今後高次元の他の LSD であっても、容易に統計分析の研究対象になる。本研究は、質が高く、2 群が LSD であるという検証しやすいデータを用いたことで、LINGO[2]と JMP[1]の組み合わせで初めて役に立つ研究を退官後に完成できたことは医学データを研究対象としたことが幸運であったと考える。
1

新村秀一 (2004).『JMP活用

2

新村秀一 (2010). 『最適線形判別関数』. 日科技連出版.

3

新村秀一 (2011). 『数理計画法による問題解決法』. 日科技連出版.

4

Shinmura S (2000b). Optimal Linear Discriminant Function using Mathematical Programming. Disertation,

統計学とっておき勉強法』．講談社．

Okayama Univ.
5

Springer1: Shinmura S (2016). The New Theory of Discriminant Analysis after R Fisher, Springer. DOI:
10.1007/978-981-10-2164-0

6 Springer2: Shinmura S (2019a) High Dimensional Microarray Data Analysis – Cancer Gene Diagnosis andMalignancy Indexes by Microarray. Springer.
2

脳画像における高次元データ解析法
佐賀大学 医学部
川口淳 田尻涼・概要
近年の医学ビッグデータ解析は，これまで個々に解析されていた複数のデータセットを統合的に解析することによって，疾患の特徴づけなどを多様な角度から同時に行っている．脳画像解析ではマルチモダルと呼ばれ，脳の形態と機能などの側面からの脳病態が評価される．
本研究では多種モダリティ脳画像間の関連性を明らかにする方法を提案した．この方法ではモダリティの情報のみならず，アウトカムの情報を組み込んで解釈のしやすい結果を導く
ように工夫した．またスパース推定を用いアウトカムに関連する変数の特定を行った．
・適用データ
US-ADNI2 にて PET と MRI、DTI を撮影した 42 名より正常群、AD 群それぞれ 10 名、計20 名をランダムにサンプリングした。この 20 名の MRI と PET、DTI 画像に対して教師付
きマルチブロック主成分法を適用し、性能の評価を行った．本データは申請を通せば使える
オープンデータである．
・前処理
sMRI 画像は SPM12 を用いてセグメンテーションを行い、標準脳に対して解剖学的標準化を行った．これにより灰白質画像(GM: Gray matter)、白質画像(WM: White matter)、脳脊髄液
画像(CSF: Cerebrospinal fluid)の 3 種画像を得る．DTI 画像は fsl を用いて FA(Fractional
Anisotropy)画像を作成、これに対して解剖学的標準化を行った．PET(使用薬剤: FDG)画像
は SPM12 を用いて、検査中に複数回撮影される画像を 1 画像にするために画像の加算を行い、解剖学的標準化を行った．以上 3 モダリティ、5 種の画像を用いて解析を行う．
・結果
教師付きマルチブロック主成分法にて、μ = 0.5とし、5 成分を求めた．各モダリティの各成分に対する寄与を見るため以下に各成分のスーパーウェイトを記す．
GM

WM

CSF

FA

PET

成分 1

0.43

0.24

0.43

0.71

0.27

成分 2

0.39

0.26

0.28

0.83

0.14

成分 3

0.37

0.21

0.37

0.80

0.22

成分 4

0.29

0.23

0.30

0.55

0.69

成分 5

0.21

0.10

0.30

0.31

0.87

表１ 各成分のスーパーウェイト
成分１～３では FA 画像が非常に大きな重みをもっており、それに付随する形で構造情報として灰白質、脳脊髄液画像が比較的大きな重みとなっている．成分 4 では FA に代わって PETの割合が大きくなっており、神経情報と機能情報をバランスよく含んでいると解釈できる．
成分 5 では PET の割合が極めて大きくなり、機能情報を主に含んだ成分になっている． また、以下にブロックウェイトとして第 1 成分の灰白質と FA の重みを元画像に表示したもの
を示す．

←萎縮傾向

図 1 成分 1, GM のブロックウェイト

図 2 成分 1, FA のブロックウェイト

得られたスーパースコアの妥当性を検証のためロジスティック回帰を行った．変数選択
は stepwise 法を用い、Null モデルを初期モデルとし、AIC 最小となるモデルを求めた
結果、次のモデルが選択された．
logit(Pr) = 𝛽0 + 𝛽1 𝑐𝑜𝑚𝑝1 + 𝛽2 𝑐𝑜𝑚𝑝2 + 𝛽3 𝑐𝑜𝑚𝑝4 + 𝛽4 𝑐𝑜𝑚𝑝5
パラメータ
推定値

𝛽0

𝛽1

𝛽2

𝛽3

𝛽4

-1.67

0.212

0.0875

0.199

0.549

SE

3.49

0.248

0.870

0.259

0.713

p値

0.633

0.395

0.314

0.444

0.441

表2

ロジスティック回帰の各パラメータ値

AIC=17.3

各成分のパラメータにおいて５％を切ったものはなかった．これは 20 例での検証のため、例数不足の可能性がある．しかし、AIC の観点では得られた 5 つスコアのうち 4 つのパラメータが AD に対して良い予測スコアとなっていることが示された．従って、教師データに対応する部位の抽出に成功したと考えられる．

Sparse Tensor Subspace method forbiological modulator selection
Heewon Park（広島大学情報科学部）
小西貞則（中央大学理工学部）
We are often interested in analysis of multiple datasets comprising the same variables asmeasured in different groups, and we call this type of data as tensor data. In this study, we consider
that “the multiple datasets have a common eigenvector structure but with different sets ofeigenvalues” (Pepler, 2014). Flury (1984) focused that “the covariance matrices of different
groups have a common basic structure, even though the underlying covariance matrices are notexactly identical across all groups”, and generalized a technique for PCA to identify a common
structure of multiple groups. Wang et al. (2011) proposed a method for common component
analysis (CCA) to extract common low-dimensional subspace, that accurately describes alldatasets of multiple groups. However, the existing methods for CCA provide fully dense common
loadings, thus tensor subspaces are constructed by all variables. It implies that tensor subspaceconstruction procedures can be disturbed by noisy features, because noisy features are inevitably
included in datasets. Moreover, the fully dense loadings lead to difficulty in interpretation of CCA
results, like in PCA (Al-Kandari and Jolliffe, 2005).
To resolve these issues, we incorporate sparsity into CCA, and propose sparse CCA to construct
sparse common subspace (i.e., sparse tensor subspace) of datasets from multiple groups. In orderto perform sparse CCA, we first focus that “ordinary CCA can be implemented by eigen analysis
of a matrix (i.e., common loadings of multiple datasets are estimated as eigenvectors of a matrix
𝑴”. We then consider that the common loadings of a matrix can be also obtained by SVD ofsquare root of the matrix 𝑴. In order to incorporate sparsity into CCA, we introduce the use ofsparse PCA (Zou et al.,2006) with non-centered matrix, i.e., the common loadings of multiple
datasets can be estimated by sparse PCA of the square root of the Gram matrix, which is a noncentered data matrix. Incorporation of sparsity into CCA leads to effective analysis of multiple
high-dimensional datasets, because the high-dimensional data inevitably contain noisy features.
Furthermore, our method provides efficient interpretation results for the constructed common
subspace (i.e., we can select crucial common-features of multiple datasets).
We apply the proposed method to construct tensor subspace of drug sensitivity-specific gene
regulatory networks, i.e., drug sensitive (resistance) subspace. It can be expected that thereconstructed drug sensitive and resistant networks on the common subspace can effectively
describe characteristics of drug sensitive and resistant networks, respectively. Thus, we consider
candidate drug identification based on the similarity test between the reconstructed drug sensitiveand resistance tensors on common subspace. For biological modulator selection, we also propose
a permutation-based statistical test to evaluate whether the connectivity of genetic network varies
with modulators. The similarity test between two network tensors is conducted based on the
proposed similarity test on not original data space but the constructed drug sensitivity-specificsubspace. We can see through the simulation studies and real data analysis that the proposed
method can effectively perform for common structure identification from multiple datasets, andthe effectiveness leads to reliable results to biological modulator selection (i.e., candidate drug
discovery).

ҼՌਪ࿦ʹ͓͚ΔηϛύϥϝτϦοΫΞϓϩʔνͷͨΊͷ৘ใྔ‫ن‬४ೋ‫ ٶ‬Յߦ ʢ౷‫਺ܭ‬ཧ‫ॴڀݚ‬ʣ

1४උपลߏ଄ฏ‫ۉ‬Ϟσϧ͸ɼજࡏ݁Ռม਺ͷपล‫ظ‬଴஋ʹର͢ΔϞσϧͱͯ͠ఏҊ͞Ε͍ͯΔɽࠓɼॲஔ͕
H छྨ͋Γɼh ൪໨ͷॲஔΛड͚ͨͱ͖ͷજࡏ݁Ռม਺Λ y (h) (∈ Rm )ɼͦͷॲஔΛड͚Δͱ 1 Ͱड͚ͳ!H (h)
͍ͱ 0 ʹͳΔϥϯμϜͳׂΓ౰ͯม਺Λ t(h) ͱ͢Δ (h ∈ {1, 2, . . . , H},= 1)ɽͦͯ͠ɼ֤જࡏ
h=1 t݁Ռม਺ʹઢ‫ܗ‬ճ‫ؼ‬ϞσϧΛ૝ఆͨ͠पลߏ଄ฏ‫ۉ‬Ϟσϧ
H
H
#
$
"
"
y=
t(h) y (h) =t(h) X (h) β + ε
h=1

h=1

Λߟ͑Δɽӈลͷ X(∈ R) ͸ಠཱม਺ߦྻɼε (∈ Rm ) ͸ฏ‫ ͕ۉ‬0m Ͱ෼ࢄ‫ڞ‬෼ࢄߦྻ͕ σ 2 Im ͷ‫ࠩޡ‬ม਺ϕΫτϧͰ͋Δɽ͜͜Ͱɼ0m ͸ m ࣍‫ྵݩ‬ϕΫτϧɼIm ͸ m ࣍‫ݩ‬୯ҐߦྻͰ͋Δɽࠨลͷ y ͸‫؍‬ଌ͞ΕΔ݁Ռม਺ͱͳΔ͜ͱʹ஫ҙ͢Δɽ͜ͷϞσϧͰ͸ɼt(h) = 0 ͳΔ H − 1 ‫ݸ‬ͷજࡏ݁Ռม਺ y (h)͕ܽଌ͍ͯ͠Δͱߟ͑ΒΕΔɽͦͯ͠ɼҰൠʹ͸ E[y (h) ] ̸= E[y (h) | t(h) = 1] Ͱ͋ΔͨΊɼ‫؍‬ଌ஋ͷΈ͔Β୯७ʹ X (h) β Λਪఆ͢ΔͱόΠΞε͕ੜͯ͡͠·͏ɽຊൃදͰ͸ɼ͜ͷόΠΞεͷআ‫͕ڈ‬ՄೳͱͳΔΑ͏ʹɼy (h) ͱ t(h) ͷަབྷม਺ϕΫτϧ z (∈ Rs ) ͕‫؍‬ଌ͞Ε͍ͯΔ͜ͱΛ૝ఆ͢Δɽ
͜͏͍ͬͨϞσϧͱม਺ʹର͠ɼ௨ৗஔ͔ΕΔ৚݅ΛԾఆ͢Δɽ·ͣ X (h) ʹ͍͕ͭͯͩɼҰൠతͳઃఆΛѻ͏ͨΊʹަབྷม਺ z Λ‫ؚ‬ΜͰ΋Α͍ (ϥϯμϜͰ͋Δ͜ͱΛ‫͢ʹͱ͜ )͢ڐ‬Δɽ·ͨɼࣜʹ͓͚Δද
!H‫ݱ‬ͷ൥ࡶ͞Λܰ‫͢ݮ‬ΔͨΊɼຊ࣭తͰ͸ͳ͍͕͜ͷಠཱม਺͸ E[ h=1 X (h)T X (h) ] = Ip ͱͳΔΑ͏ʹ‫ج‬४Խ͞Ε͍ͯΔ΋ͷͱ͢Δɽ͜ͷ‫ج‬४Խͷ༗ແʹΑΒͣɼҎ߱Ͱ݁Ռͱͯ͠ಋ͔ΕΔ‫ج‬४͸ಉ͡΋ͷͱͳΔɽ࣍ʹɼ্ड़ͷόΠΞεͷআ‫ڈ‬ΛՄೳͱ͢ΔͨΊͷ΋ͷͰ͋Δɼऑ͘ແࢹͰ͖ΔׂΓ౰ͯ৚݅
(h)

m×p

y (h) ⊥ t(h) | z(h ∈ {1, 2, . . . , H})

ΛԾఆ͢Δɽ͜ͷ৚݅ʹ͓͚Δ y (h) ͸ ε ʹม͑ΒΕΔ͜ͱʹ஫ҙ͢Δɽͯ͞ɼ͜ͷϞσϧʹ͕ͨ͠͏N ‫ݸ‬ͷαϯϓϧ͕͋Δͱ͠ɼୈ i αϯϓϧͷม਺ʹ͸ఴ͑ࣈ i Λ෇͚Δ͜ͱʹ͢Δɽͦͯ͠ɼỹ (h) =
(h)T
(h)T
(h)T
(h)
(h)T
(h)T
(h)T
(y1 , y2 , . . . , yN )T ɼT (h) = diag{ti Im }ɼX̃ (h) = (X1 , X2 , . . . , XN )T ɼε̃ = (εT1 , εT2 , . . . ,T TεN ) ͱ͠ɼ·ͱΊͯ

ỹ =

H
"

T (h) ỹ (h) =

h=1

H
"

h=1

ͱද͢͜ͱʹ͢Δɽͨͩ͠ɼαϯϓϧ͸ಠཱɼͭ·Γ
(h)

(h)

(h)

(h)

#
$
T (h) X̃ (h) β + ε̃(ti , Xi , εi , zi ) ⊥ (tj , Xj , εj , zj )(i ̸= j, h ∈ {1, 2, . . . , H})(h)Ͱ͋Δ͜ͱΛԾఆ͢Δɽ͜ΕΑΓ౰વ yi ⊥ yj (i ̸= j) Ͱ͋Δɽ·ͨɼ௨ৗͷճ‫ؼ‬ϞσϧͷΑ͏ʹ Xi ͱεi ͸ಠཱͰ͋ΔͱԾఆ͓ͯ͘͠ɽ
જࡏ݁Ռม਺ y (h) ͱަབྷม਺ z ͷؔ܎Λਖ਼͘͠ϞσϦϯάͰ͖͍ͯΕ͹ɼແࢹͰ͖ΔׂΓ౰ͯ৚݅ͷ΋ͱͰ‫؍‬ଌ஋ͷΈ͔Β y (h) ͷपล‫ظ‬଴஋ΛҰகਪఆ͢Δ͜ͱ͕Ͱ͖Δ͕ɼҰൠʹ͜ͷϞσϦϯά͸೉͍͠ɽ(h)(h)ͦ͜Ͱۙ೥͸ɼ͜ͷਖ਼͍͠ϞσϦϯάΛඞͣ͠΋ඞཁͱ͠ͳ͍ɼ܏޲είΞ ei ≡ P(ti = 1 | zi ) Λ༻͍ͨηϛύϥϝτϦοΫΞϓϩʔνΛ༻͍ΒΕΔ͜ͱ͕ଟ͍ɽຊൃදͰ͸ɼ͜ͷΞϓϩʔνʹ͓͍ͯ୅දతͳ
IPW ਪఆΛѻ͏ɽ͜ͷਪఆ๏Ͱ͸ɼ‫؍‬ଌ஋ʹ܏޲είΞͷ‫਺ٯ‬ΛॏΈͱ͔͚ͯ͠Δ͜ͱͰٙࣅతʹܽଌ஋(h)(h)Λ෮‫͠ݩ‬ɼͦͷ‫ʹޙ‬௨ৗͷਪఆΛ͓͜ͳ͏ɽ۩ମతʹ͸ɼॏΈߦྻ W (h) ≡ diag(ti Ir /ei ) Λ༻͍ͯॏΈ෇͖ೋ৐ଛࣦؔ਺Λ
H #
"

h=1

ỹ − X̃ (h) β$T#$W (h) ỹ − X̃ (h) βͱఆٛ͠ɼ͜ΕΛ β ʹ͍ͭͯ࠷খԽ͢Δ͜ͱͰ IPW ਪఆྔ
%H
&−1 H
"
"
IPW
(h)T
(h)
(h)
β̂
≡
X̃
W X̃
X̃ (h)T W (h) ỹ
h=1

h=1

Λ༩͑ΔɽIPW ਪఆྔ͸ɼऑ͘ແࢹͰ͖ΔׂΓ౰ͯ৚݅ͷ΋ͱͰҰகੑΛ΋ͭɽ

1

(1)

2

Ϟσϧબ୒‫ج‬४σʔλʹܽଌ͕͋Δ৔߹ͷ MSE ͱͯࣗ͠વʹߟ͑ΒΕΔ΋ͷΛೋछྨߟ͑Δɽ·ͣҰछྨ໨ͱͯ͠
mean weighted squared error Λ
H"

MwSE =

E

h=1
H
"

=

E

h=1

−

'#
'#

H
"

H
"

ỹ − X̃ (h) β̂

'#

E

h=1

+2

'
($T
#
'
($(
X̃ (h) β̂ − E ỹ (h) | X̃ (h)W (h) X̃ (h) β̂ − E ỹ (h) | X̃ (h)

E

h=1

$T

#
$(
W (h) ỹ − X̃ (h) β̂

'
($T
#
'
($(
ỹ − E ỹ (h) | X̃ (h)W (h) ỹ − E ỹ (h) | X̃ (h)

'#

'
($T
#
'
($(
ỹ − E ỹ (h) | X̃ (h)W (h) X̃ (h) β̂ − E ỹ (h) | X̃ (h)ͱఆٛ͢Δɽ͜͜Ͱɼ௨ৗͷ Cp ‫ج‬४ͷಋग़ʹͳΒ͍ɼఆٛ‫ʹޙ‬͸ࡾͭͷ߲ʹ෼ղ͍ͯ͠Δɽ͜͜ͰͷॏΈ෇͖ೋ৐࿨͸ɼ‫؍‬ଌ஋ʹॏΈΛ෇͚ͯ෮‫ͨ͠ݩ‬σʔλͷ‫ظ‬଴஋ͱͦͷਪఆ஋ͷࠩͷೋ৐࿨ͱΈͳͤɼͦͷॏΈ͸ IPW ਪఆͰ΋༻͍ΒΕ͍ͯΔ΋ͷͰ͋Δɽ࣮ࡍɼ෼ղ‫ޙ‬ͷୈҰ߲໨͸ (1) ͷ‫ظ‬଴஋Λͱͬͨ΋ͷ
Ͱ͋Δɽͭ·Γɼਪఆྔͷಋग़ͱਪఆྔͷ‫ࠩޡ‬ධՁʹಉछͷଛࣦؔ਺Ͱߟ͓͑ͯΓɼͦͷ఺ͰࣗવͰ͋Δɽ࣍ʹɼೋछྨ໨ͱͯ͠ mean unweighted squared error Λ

MuSE =

H
"

E

h=1

=

H
"

h=1

−

E

'#
'#

H
"

E

h=1

+2

H
"

h=1

'
($T
#
'
($(
X̃ (h) β̂ − E ỹ (h) | X̃ (h)T (h) X̃ (h) β̂ − E ỹ (h) | X̃ (h)ỹ − X̃ (h) β̂

'#

E

$T

#
$(
T (h) ỹ − X̃ (h) β̂

'
($T
#
'
($(
ỹ − E ỹ (h) | X̃ (h)T (h) ỹ − E ỹ (h) | X̃ (h)

'#

'
($T
#
'
($(
ỹ − E ỹ (h) | X̃ (h)T (h) X̃ (h) β̂ − E ỹ (h) | X̃ (h)ͱఆٛ͢Δɽ͜͜Ͱ͸‫؍‬ଌ஋ͦͷ΋ͷͷσʔλͷ‫ظ‬଴஋ͱͦͷਪఆ஋ͷࠩͷೋ৐࿨Λߟ͍͑ͯΔɽ‫؍‬ଌ͞Εͨ΋ͷʹର͢Δਪఆਫ਼౓ͷ޲্ʹΑΓ஫ྗ͍ͨ͠ͱ͍͏ߟ͑ํͷ΋ͱͰ͸ɼͪ͜Βͷଛࣦؔ਺ͷํ͕ࣗવ
ͱ͍͑Δɽ͍͏·Ͱ΋ͳ͘ɼW (h) ΍ T (h) ͸ඇෛఆ஋ରশߦྻͰ͋Γɼ྆ MSE ͱ΋ϊϧϜʹͳ͍ͬͯΔɽ
௨ৗͷ Cp ‫ج‬४ͷಋग़ʹͳΒ͍ɼ෼ղͨ͠΋ͷͷୈҰ߲໨ͷ‫ظ‬଴஋Λ֎͠ɼϞσϧʹґଘ͠ͳ͍ୈೋ߲໨
Λແࢹ͠ɼE[ỹ (h) | X̃ (h) ] = X̃ (h) β ͱͯ͠ୈࡾ߲໨Λ઴ۙධՁͨ͠΋ͷΛҼՌਪ࿦ʹ͓͚Δ Cp ‫ج‬४ͱ͢
Δɽ઴ۙධՁͰ͸ɼ‫ظ‬଴஋ͷத਎ʹରͯ͠઴ۙతʹओ߲ͱͳΔ΋ͷΛऔΓग़͠ɼͦͷ‫ظ‬଴஋Λཅʹ‫ٻ‬ΊΔɽ
ͦͯ͠ɼMwSE ͋Δ͍͸ MuSE ͔Βಋ͔ΕΔ‫ج‬४ΛͦΕͧΕ wCp ͋Δ͍͸ uCp ͱද͢͜ͱʹ͢Δɽ

Theorem 1. ܏޲είΞ͕‫ط‬஌Ͱ͋ΔέʔεͰ IPW ਪఆྔΛ༻͍ͨͱ͖ͷ Cp ‫ج‬४͸
wCp =

H #
"

h=1

uCp =

H #
"

h=1

ỹ − X̃ (h) β̂ IPWỹ − X̃ (h) β̂ IPW

$T

$T

)
*
H
#
$
"
1
W (h) ỹ − X̃ (h) β̂ IPW + 2
E (h) εT X (h) X (h)T ε ,
e
h=1

#
$
T (h) ỹ − X̃ (h) β̂ IPW + 2σ 2 pͰ༩͑ΒΕΔɽ!N !H (h) (h)
(h)
(h)
(h)T
wCp ͷϖφϧςΟ߲͸ཅʹ‫͖Ͱࢉܭ‬Δͱ͸‫ݶ‬Βͳ͍͕ɼ i=1 h=1 ti (yi − Xi β̂ IPW )T Xi Xi
(h)
(h)
(h)2
(yi − Xi β̂ IPW )/(N ei ) ͷΑ͏ʹҰகਪఆྔΛ༰қʹ༩͑Δ͜ͱ͕Ͱ͖Δɽ

2

ํ޲ぶがの〠ର『぀von Mises ෼෍〣฼਺ਪఆ〠〙⿶〛࡞ଜ ‫∗ لݐ‬1༄ຊ෢ඒ †〝「〛ఏҊ『぀ɽvon Mises ෼෍〣ࣗવ฼਺〤ɼ θ !

〤」〶〠(θ1 , θ2 ) ! (τ cos η, τ sin η) 〜⿴〿ɼରԠ『぀े෼౷‫ܭ‬ํ޲〠ؔ『぀෼෍〝「〛ɼ࠷〷༗໊〟〷〣〠 vonྔ〤 x ! (cos x, sin x) 〜⿴぀ɽࣄલ෼෍〠〤ࢀরࣄલMises ෼෍⿿⿴぀ɽ〈〣෼෍〤ɼࢦ਺෼෍଒〠ଐ「〛⿼ີ౓ぇ༻⿶぀ (Garvan & Ghosh, 1997; Robert, 2007)ɽ

%〿ɼํ޲ぶがの〠ؔ『぀୅දత〟෼෍〠Ґஔ『぀ɽ〈〣෼෍〤 2 〙〣む゘ゐがのぇ࣋〘〛⿼〿ɼ〒ぁ〓ぁํ޲む゘ゐがのɼूத౓む゘ゐがの〝「〛஌〾ぁ〛⿶぀ɽํ޲ぶがの〠ؔ『぀‫ڀݚ‬〤ɼྫ⿺〥ɼ෩޲਺ θ1 〣ࣄ‫ޙ‬ฏ‫ۉ‬〤ɼ぀む゘ゐがのਪఆ〣‫〛「〝ڀݚ‬〤ɼSchou (1978) 〠〽

∫θ̂1 ! cos x぀पล࠷໬ਪఆྔ〣ఏҊ〟〞⿿‫ڍ‬〆〾ぁ぀⿿ɼ〒〣ଞ〣〝〟぀ɽಉ༷〠ɼྫ〤গ〟⿶ɽ
Ұํɼࢦ਺‫ܕ‬෼෍଒〠⿼⿶〛ࣗવ฼਺〣ࣄ‫ޙ‬ฏ‫ۉ‬〤θ̂2 ! sin x

࠷దੑぇอ〙 (Yanagimoto & Ohnishi, 2009)ɽ〈ぁぇ〶぀ఏҊਪఆྔ〤ɼ『〈〝⿿‫ظ‬଴《ぁ぀ɽη̂ ! tanvon Mises ෼෍ඪຊ X1 , . . . , X n 〤ɼ࣍〣֬཰ີ౓ؔ਺ぇ࣋〙 vonτ̂ !

Mises ෼෍ VM(η, τ) (η ∈ (−π, π), τ > 0) ⿾〾〣ແ࡞〝〟぀ɽҝඪຊ〝『぀ɽ

!
"1
exp τ cos(x − η) , x ∈ (−π, π)2πI0 (τ)
i!1 cos(x i ), S(x) !

3

(2)

! x,

θ̂12 + θ̂22 !

∫τπ m (τ|x)A(τR(x))dτ,(3)〈〈〜ɼ τ 〠ؔ『぀ MLE 〤ҎԼ〜ಘ〾ぁ぀ɽ(1)τ̂ML ! A−1x̄ ! tan−1 (S(x)/C(x)) 〜⿴぀ɽ

(

R(x)
n

)

〳〔ɼSchou (1978) 〤ɼ τ 〠ؔ「〛ɼ࣍ࣜぇຬ〔『MMLE, τ̂MML ぇఏҊ「〛⿶぀ɽ√nA( τ̂MML ) ! R(x)A( τ̂MML R(x)), R(x) > n√〳〔ɼ0 ≤ R(x) ≤ n 〣〝 ɼ τ̂MML ! 0 〜⿴぀ɽ《Yanagimoto and Ohnishi (2009) 〠〽〿ɼ࠷దੑ⿿
อূ《ぁ〔ਪఆྔ〜⿴぀ࣗવ฼਺〣ࣄ‫ޙ‬ฏ‫ۉ‬ぇɼvon–Mises ෼෍〣฼਺ (η, τ) 〣ਪఆ〠ର『぀༗ޮ〟ਪఆྔ
†

θ̂1

)

ぶ゚む゘ゐがの〣ࣄ‫ޙ‬ฏ‫ۉ‬ʢPMʣぇରর〝『぀ɽ

ఏҊਪఆྔ

∗

θ̂2

ʢMLEʣ
ɼपล࠷໬ਪఆྔʢMMLEʣ (Schou, 1978)ɼゑぷ゚ R(x) 〝ඪຊฏ‫ ޲ํۉ‬x̄ ぇߟ⿺぀〝ɼ〜⿴〿ɼR(x) ! 0 〣〝(〣ൺֱ࣮‫ݧ‬ぇߦ⿸ɽ‫ط‬ଘਪఆྔ〝「〛〤ɼ࠷໬ਪఆྔ

i!1 sin(x i ) 〝「ɼ

C 2 (x) + S2 (x),*−1ఏҊਪఆྔ〣ੑೳぇௐࠪ『぀〔〶〠ɼ‫ط‬ଘਪఆྔ〝#nC(x) ! R(x) cos x̄, S(x) ! R(x) sin x̄ ぇຬ〔『߹੒よぜ$τπ m (τ|x)A(τR(x))dτ,

4 ਪఆྔ〣ൺֱ〈〈〜ɼؔ਺ I0 ( · ) 〤ୈ 1 छม‫ܗ‬よひな゚ؔ਺〜⿴぀ɽ

R(x) !

∫〜⿴぀⿾〾ɼ θ ! (θ1 , θ2 ) ! (τ cos η, τ sin η) 〽〿ɼ‫ٻ‬⿺぀〝ɼ࠷దੑ〣‫〾⿾఺؍‬ɼఏҊਪఆྔ〤ྑ⿶ੑೳぇࣔ〈〈〜ɼC(x) !

τπ m (τ|x)A(τR(x))dτ,⿿ಘ〾ぁ぀ɽ〈ぁ〾〤ɼむ゘ゐがの η, τ 〠ؔ「〛ෆมvon Mises ෼෍〠ద༻「ɼఏҊਪఆྔ〝『぀〈〝ぇߟ#n.

η 〠ґଘ「〟⿶〈〝⿾〾ɼҎ߱ɼπR (τ) 〝ද『ɽࣗવ฼〜Ԡ༻⿿‫ظ‬଴《ぁ〛⿶぀ɽ〳〔ɼvon Mises ෼෍〠ؔ『p(x;η, τ) !

& 1/2〈〈〜ɼA(x) ! I1 (x)/I0 (x) 〜⿴぀ɽ〈ぁ〤む゘ゐがの〣෼ੳ (Mardia & Jupp, 1972) 〟〞⿿⿴〿ɼۙ೥ɼ޿⿶෼໺2

A(τ)− A2 (τ)τπR (η, τ) ∝ 1 −〾〠ɼ τ 〣 PM 〤࣍ࣜ〜༩⿺〾ぁ぀ɽ๏੓େֶཧ޻ֶ෦, ˟184-8584 ౦‫౎ژ‬খۚҪࢢֿ໺ொ 3–7–2
౷‫਺ܭ‬ཧ‫ॴڀݚ‬, ˟190-0014 ౦‫ཱ౎ژ‬઒ࢢ྘ொ 10–3

τ̂PM !
1∫

τπ m (τ | x)dτ

5

Table1: ଛࣦ KL e 〣〷〝〜〣゙とぜൺֱ

Prop.
MLE
MMLE
PM

4

τ^

3

n

τ

Prop.

MMLE

MLE

PM

10

0
0.3
1
3
10

0.55860.62980.98581.104
1.055

0.51650.6686
1.169
1.065
1.052

1.043
1.055
1.065
1.100
1.087

0.90110.89680.9540
1.075
1.049

30

0
0.3
1
3
10

0.51750.6925
1.064
1.033
1.021

0.49230.8573
1.045
1.023
1.019

1.009
1.014
1.000
1.028
1.032

0.87870.8517
1.018
1.024
1.019

100

0
0.3
1
3
10

0.50650.9277
1.031
1.012
1.022

0.4822
1.122
1.018
1.009
1.020

1.004
1.019
1.000
0.9991
1.025

0.87430.9126
1.020
1.009
1.020

2
1
0
0

1

2

3

R (x )

Figure1: つアゆ゚つぐど n ! 3 〠⿼々぀ɼR(x) 〠ର『぀ τ 〣ਪఆྔ〣ৼ぀෣⿶: ࣮ઢ⿿ఏҊਪఆྔʢProp.ʣɼ఺ઢ⿿ MLEɼഁઢ⿿ PMɼ௕ഁઢ⿿ MMLE ぇࢦ『ɽ4.1ਪఆྔ〣ৼ぀෣⿶

ूத౓む゘ゐがの τ 〣〒ぁ〓ぁ〣ਪఆྔ〤ɼR(x) 〠ରԠ「〛ਪఆ《ぁ぀ɽ〒〈〜ɼ〳』ɼR(x) 〠Ԡ」〔ਪఆྔ τ̌ 〣ৼ぀෣⿶〠〙⿶〛‫ݕ‬౼『぀ɽつアゆ゚つぐど〰〱ಉ౳〜⿴぀ɽ〳〔ɼn ! 30, 100 〣৔߹〠⿼⿶〛〤ɼ

ぇ n ! 3 〝‫ݻ‬ఆ「ɼR(x) ぇมԽ《【〔〝〣ਪఆྔぇτ ! 1 〣〝 〣݁Ռ〤‫ٯ‬స「〛぀〷〣〣ɼ〰〱ಉ」܏޲

exact 〠‫݁〔〶ٻ‬Ռぇ Figure 1〠ࣔ『ɽਤ〽〿ɼҎԼ〣√〈〝⿿‫⿺ݴ‬぀ɽMMLE 〤 τ ! 3 〜ંぁ఺⿿⿴぀ɽ〳〠⿴぀ɽ〳〔ɼMLE 〝ൺ〮぀〝ɼఏҊਪఆྔ〤〰〝え〔ɼPM 〤 R(x) ! 0 〣〝 〠 τ̂ ! 0 〜⿴぀ɽR(x) ⿿େ〠 τ ! 0 〠ۙ⿶〰〞ɼఏҊਪఆྔ〣〰⿸⿿ྑ⿶܏޲〠〞〣৚݅〣〷〝〜༏ӽ「〛⿶぀ɽPM 〝ൺֱ『぀〝ɼಛ。〟぀〝ɼఏҊਪఆྔ〝 MMLE 〤ඇৗ〠ۙ⿶஋ぇऔ⿴぀ɽ〈ぁ〾〤ɼvon Mises Fisher ෼෍〟〞〣ؔ࿈෼〿ɼPM 〷〳〔ۙ〚⿶〛⿶぀ɽMLE 〤ఏҊਪఆྔ〽〿෍〠⿼々぀ਪఆྔ〣ߏ੒〠ؔ「〛ɼ〒〣‫〟〝ૅج‬぀݁Ռ〷ৗ〠େ〜⿴぀ɽ《〾〠ɼ〈〈〜〤ແ৘ใࣄલ෼෍ぇԾఆ「〛゙

4.2

⿶஋ぇऔ぀ɽ

゙とぜൺֱ

とぜൺֱぇߦ〘〔⿿ɼ৘ใ〣⿴぀ࣄલ෼෍〜〣બ୒〹〒ఏҊਪఆྔ〣ੑೳぇ゙とぜൺֱ〠〽〘〛‫『ূݕ‬぀ɽ〣ਪఆྔ〣ߏ੒〠〙⿶〛〷༗ӹ〟݁Ռ〜⿴぀〝‫ظ‬଴〜゙とぜ〤ස౓࿦〣‫『ࢉܭ〾⿾఺؍‬぀ɽ〙〳〿ɼR( τ̌) !

∫぀ɽ

ReferencesL( τ̌, τ)p(x|τ)dx

Garvan, C. W. & Ghosh, M. (1997, Dec.). Noninformative
priors for dispersion models. Biometrika, 84(4), 976-982.
Mardia, K. V. & Jupp, P. E. (1972). Directional statistics.
John Wiley & Sons, LTD.
Robert, C. (2007). The bayesian choice: from decision-theoretic
foundations to computational implementation. Springer Science & Business Media.
Schou, G. (1978). Estimation of the concentration parameter in von mises-fisher distributions. Biometrika,
65(2), 369-377.
Yanagimoto, T. & Ohnishi, T. (2009). Bayesian predictionof a density function in terms of e-mixture. Journal of
Statistical Planning and Inference, 39, 3064-3075.

〈〈〜ɼ p(x|τ) 〤 η ! x̄ 〝「〔 von Mises ෼෍〣ඪຊີ౓〜⿴〿ɼL( τ̌, τ) 〤ଛࣦؔ਺ぇද「ɼ〈〈〜〤ɼe てとふわ〠⿼々぀ KL はぐみがでこアと, KL e 〝ɼm てとふわ〠⿼々぀ KL はぐみがでこアと, KL m ɼೋ৐ଛࣦ,MSE ぇ༻⿶぀ɽ
‫ط‬ଘਪఆྔぇରর〝「〔ఏҊਪఆྔ〝〣゙とぜൺ
ֱぇɼてゎゔ゛がてゖア〠〽〘〛ߦ⿸ɽてゎゔ゛がてゖア〠⿼⿶〛ɼํ޲む゘ゐがの η ! 0 〠‫ݻ‬ఆ「ɼつアゆ゚つぐど n ! {10, 30, 100}, ूத౓む゘ゐがの

τ ! {0, .3, 1, 3, 10} 〣ߏ੒〜ɼ〒ぁ〓ぁ〠⿼⿶〛 10,000

つアゆ゚ぇੜ੒「〛ߦ⿸ɽ݁Ռぇ Table 1 〠ࣔ『ɽTable 1 〽〿ɼఏҊਪఆྔ〤ɼMMLE ぇରর〝『぀〝ɼ n ! 10 〠⿼⿶〛ɼ τ ! 0 〜〤 MMLE 〠ྼ぀〷〣〣ɼ τ ! 0.3 ⿼〽〨 τ ! 1 〜〤༏ӽ「〛⿼〿ɼ τ ≥ 3 〜

2

A Practicable Estimation of Mean Squared
Prediction Error in Small Area Estimation
Masayo Y. Hirose
The Institute of Mathematics for Industry, Kyushu University

Abstract
The empirical best linear unbiased prediction (EBLUP) estimator
is utilized for eﬃcient inference in various research areas, especially forsmall-area estimation. In order to measure its uncertainty, we generally
need to estimate its mean squared prediction error (MSPE). In thisstudy, we seek a MSPE estimation method achieving several desired
properties under aggregated level model.

1

Introduction

There has been high demand for reliable statistics on smaller geographic
areas and sub-populations where large samples are not available. Even insuch a situation, an explicit model-based approach can achieve more accurate estimates by borrowing strength from related areas.
The Fay–Herriot model (Fay and Herriot, 1979), in particular, is widelyused as an aggregated level model for small-area inference as follows:
For i = 1, . . . , m,ind

Level 1 :yi |θi ∼ N (θi , Di );ind

Level 2 :θi ∼ N (x′i β, A).

(1)

The level-1 model takes into account the sampling distribution of the direct
estimator yi for the ith small area. The true small-area mean for the ith
area, denoted by θi , is linked to the area-specific auxiliary variables xi =(xi1 , · · · , xip )′ in the level-2 model. In practice, the coeﬃcient vector β in
Rp and the model variance parameter A in this linking model are unknown.
The sampling variance Di is often assumed to be known.
1Since A and β are unknown in practice, the empirical best linear unbiased
predictor (EBLUP) of θi is generally used for small-area inference. It would
also be quite important to measure the MSPE of EBLUP as its uncertainty.
For small-area inference, its MSPE needs to be estimated with high accuracy.
In this study, we consider what a desirable properties of practicable
MSPE estimation. We also compare the performance of several MSPE estimators in our investigated class and others through a Monte Carlo simulation study. This study has been published in Hirose (2019).

Acknowledgement
This research was supported by Grant-in-Aid for Research Activity Start-up,
JSPS Grant Number 26880011.

References[1] Fay, R.E. and Herriot, R.A. (1979). Estimates of income for small
places: an application of James-Stein procedures to census data, Journal of the American Statistical Association 74, 269-277[2] Hirose, M.Y. (2019). A Class of General Adjusted Maximum likelihood
Methods for Desirable Mean Squared Error Estimation of EBLUP under the Fay-Herriot Small Area Model. Journal of Statistical Planningand Inference, 199, 302-310.

2

前震の統計的判別に基づく地震予測
統計数理研究所

野村俊一

任意の地震群は，群内の地震のマグニチュード差や時空間的距離の分布形に見られる特徴及び地域性に基づき，前震群である確率を本震が起こる前に評価することができる．図 1 に示す Single-link 法[1]に従い地震群を構成し，構成した地震群ごとに群内の最大地震を本震，本震の前と後の群内地震をそれぞれ前震，余震と定義する．そして，活動中の地震群が 30 日以内に本震を起こす前震群である確率を評価するための統計モデルを構築し，その識別性能を検証する．
地震群のマグニチュード差や時空間的距離に関する特徴は，地震群内の地震数 N と最大マグニチュー
ド M1 に依存して変化するものの，その 2 つが同条件の下では前震群と他の地震群の間に相対的な傾向の違いが現れる[2]．そこで，群内の最大マグニチュード M1 と二番目に大きいマグニチュード M2 の差
ΔM，地震群の継続期間 T（日），群内の平均震央間距離 D（km），地震群の中心経度 X，緯度 Y の位置
情報（地域性）を前震識別のための特徴量として，ある地震群について，群内の最新地震から 30 日以内の本震発生確率を本稿での前震確率と定義し，次式に基づいて評価する:
logit p = g(X, Y) + f1(N, M1, ΔM) + f2(N, M1, T) + f3(N, M1, D) + α0 + αs
ここで，logit p = log{p/(1－p)}は p の取る値を 0 < p < 1 の範囲へ限定するロジット関数であり，このとき右辺の各項は特徴量による前震確率の利得あるいは前震群と他の地震群の対数オッズ比と解釈される．
右辺第 1 項 g(X, Y) は地域性のみに依存したベースラインであり薄板スプライン関数により推定される．
また，f1(N, M1, ΔM), f2(N, M1, T), f3(N, M1, D)は各特徴の組合せに基づくオッズの自然対数を，それぞれ 3
次スプライン関数により表現した項である．また α0 は定数項であり，さらに，特異な地震群の影響による推定の不安定化を抑えるために，以上の項に含まれない地震群固有の変量効果 αs を加えた．
1926 年 1 月 1 日から 1999 年 12 月 31 日までの気象庁震源カタログによる上の評価式の推定結果を
図 2 に示した．図 2(d)の地域差の推定結果からは，主に伊豆半島沖と東北沖において前震確率が高くなっている．また，図 2(a)からはマグニチュード差 ΔM が小さいほど前震確率が高く，図 2(b)からは継続
期間 T が短いほど前震確率が高く，図 2(c)からは平均震央間距離が長いほど若干前震確率が高くなっており，この傾向は地震数 N によって変化していることがわかる．なお，図 2 に共通する特徴として，最大マグニチュード M1 が小さいほど前震確率が高くなっているが，これは，M1 が小さいほど対象とする
本震マグニチュードの下限が下がるためである．
以上で推定された前震確率の評価式を，2000 年 1 月 1 日から 2017 年 10 月 31 日までの同カタログに適用し，群内地震数が N = 2, 4, 8 に達した各地震群について前震確率を評価した．地震群ごとの前震確率の評価値を 10%区切りの階級に分けて集計し，各階級中の実際の前震割合すなわち適中率と比較した結果を表 1 に示す．適中率は，該当の地震群数が少ない階級を除いて，各階級の前震確率評価値と概ね整合的である．また，60%超の高い前震確率を得たのは伊豆半島沖の地震群であり，実際に M6 前後の本震を起こしている．ただし，たとえば 2016 年熊本地震の前震系列のように群内地震数が 10 を超えてくると，事例数の低下により前震確率の評価は不安定なものとなる．今後は下限マグニチュードの引
き下げや海外のカタログの活用を検討して事例数の充足を図りたい．

参考文献[1] Frohlich, C., and Davis. S. D. (1990). Single-link cluster analysis as a method to evaluate spatial
and temporal properties of earthquake catalogues, Geophys. J. Int., 100, pp.19-32.
[2] Ogata, Y., Utsu ,T. and Katsura, K. (1995). Statistical features of foreshocks in comparison withother earthquake clusters, Geophys. J. Int., 121, pp.233-254.

図 1.

図 2.

Single-link 法による前震群候補の構成方法．

地震数 N = 4 の前震群（赤点）と他の地震群（青点）の分布の違い．両者の分布の違いから推定された確率利得（対数オッズ比）をカラースケールで示している．

表 1.

2000 年以降の前震確率評価値の度数分布と適中率（＝前震群数／全地震群数）．

前震確率の評価値 0-10% 10-20% 20-30% 30-40% 40-50% 50-60% 60-70% 70-80% 合計
全地震群数
463
483
227
64
14
2
0
0 1253
N = 2 前震群数
32
74
47
20
4
2
0
0
179
適中率
7% 15% 21% 31% 29% 100% N/A N/A 14%
全地震群数
151
98
54
22
14
5
1
2
347
N = 4 前震群数
14
16
10
8
2
1
1
2
54
適中率
9% 16% 19% 36% 14% 20% 100% 100% 16%
全地震群数
82
17
10
2
6
5
2
0
124
N = 8 前震群数
7
3
4
1
3
2
2
0
22
適中率
9% 18% 40% 50% 50% 40% 100% N/A 18%

penalized LSA estimator

P-O estimatorsimulation

penalized LSA estimator

P-O estimatorsimulation

setting

B = (Ω, F , F, P ), F = (Ft )t∈R+ ɿa stochastic basisΘ ⊂ Rp ɿύϥϝʔλۭؒΘ ͸༗ք։ू߹θ∗ ɿਅ஋ɼθ∗ ∈ Θ ͔ͭ 0 ͷ੒෼͕ p0 ‫ݸ‬Penalized least squares approximation methods

θ1∗ ̸= 0, ..., θp∗0 ̸= 0ɼθp∗0 +1 = ... = θp∗ = 0

rT ɿT ͷؔ਺ɽଟ͘ͷ৔߹͸ rT−1 =
ླ໦ɹ୓ւ

√

T Λߟ͑Δɽ00ϕΫτϧ v ∈ Rp ΛαϒϕΫτϧ vJ 1 ∈ Rp ͓Αͼ vJ 0 ∈ Rp−p Λ༻͍ͯ࣍ͷΑ͏ʹද͢͜ͱʹ͢Δɽ
!
"
v 1
v= J
vJ 0

౦‫ژ‬େֶେֶӃ਺ཧՊֶ‫ڀݚ‬Պ

2019 ೥ 8 ݄ 30 ೔ಉ༷ʹ p × p ߦྻ M Λ

M=

!

MJ 11
MJ 01

MJ 10
MJ 00

"

ͷΑ͏ʹද͢͜ͱʹ͢Δɽ

ླ໦ɹ୓ւ (౦‫ژ‬େֶେֶӃ਺ཧՊֶ‫ڀݚ‬Պ)γϯϙδ΢Ϝʮߴ࣍‫ݩ‬ෳࡶσʔλͷ౷‫ܭ‬ϞσϦϯάʯ

penalized LSA estimator

2019 ೥ 8 ݄ 30 ೔

P-O estimator

1 / 10simulation

definitionླ໦ɹ୓ւ (౦‫ژ‬େֶେֶӃ਺ཧՊֶ‫ڀݚ‬Պ)γϯϙδ΢Ϝʮߴ࣍‫ݩ‬ෳࡶσʔλͷ౷‫ܭ‬ϞσϦϯάʯ

penalized LSA estimator

P-O estimator

2019 ೥ 8 ݄ 30 ೔

2 / 10simulation

main result

Objective function(q)

QT (θ) = (θ − θ̃)′ Ĝ(θ − θ̃) +

p
#
j=1

κjT |θj |q ,Ծఆ 1

rT−1 (θ̃ − θ∗ ) = Op (1)ɽθ̃ɿॳ‫ظ‬ਪఆྔĜ →p Gྫ͑͹͋Δϩεؔ਺ LT (θ) ʹରͯ͠ɼθ̃ ∈ argminθ LT (θ) ͳΔ΋ͷG ∈ Rp×p ɿ͋ΔʢϥϯμϜͳʣਖ਼ఆ஋ରশߦྻ

p×p

Ĝ ∈ R ɿʢϥϯμϜͳʣਖ਼ఆ஋ରশߦྻ
κjT = αT |θ̃j |−γTheorem 1 (rT−1 -consistency and selection consistency)γ ɿఆ਺ɼγ > −(1 − q)Ծఆ 1 ͷ΋ͱͰɼྫ͑͹ɼγ = 1ɽ−(2−q+γ)(αT )T ɿdeterministic ͳ਺ྻɼrT2ྫ͑͹ɼαT = rT
ɽ

αT → ∞,−1rT
αT = o(1)rT−1 (θ̂(q) − θ∗ ) = Op (1),(q)

P [θ̂J 0 = 0] → 1͕੒Γཱͭɽ

penalized LSA (Least Squares Approximation) ਪఆྔ θ̂(q)(q)θ̂(q) ∈ argmin QT (θ)θ∈Θ̄ླ໦ɹ୓ւ (౦‫ژ‬େֶେֶӃ਺ཧՊֶ‫ڀݚ‬Պ)γϯϙδ΢Ϝʮߴ࣍‫ݩ‬ෳࡶσʔλͷ౷‫ܭ‬ϞσϦϯάʯ

penalized LSA estimator

2019 ೥ 8 ݄ 30 ೔

P-O estimator

3 / 10simulation

main resultླ໦ɹ୓ւ (౦‫ژ‬େֶେֶӃ਺ཧՊֶ‫ڀݚ‬Պ)γϯϙδ΢Ϝʮߴ࣍‫ݩ‬ෳࡶσʔλͷ౷‫ܭ‬ϞσϦϯάʯ

penalized LSA estimator

P-O estimator

2019 ೥ 8 ݄ 30 ೔

4 / 10simulation

main result

Ծఆ 2Definition 3֬཰աఔ X = {XT }T ͕ L∞− -༗քͱ͸ɼsupT E[|XT |p ] < ∞, ∀p ≥ 1 ͷͱ͖ʹ͍͏ɽ1

rT−1 (θ̃ − θ∗ ) →ds G− 2 ζζ ɿp มྔඪ४ਖ਼‫཰֬ن‬ม਺ɼG ͱಠཱ

Ծఆ 3→ds ɿG-stable convergenceG ɿσ-fieldɼσ(G) ⊂ G ⊂ F{Ĝ}T ɼ{Ĝ−1 }T ͓Αͼ {rT−1 (θ̃ − θ∗ )}T ͸ L∞− -༗քɽTheorem 2 (Asymptotic normality)$%
G = Ip0 (GJ 11 )−1 GJ 10 ͱ͓͘ɽ͜ͷͱ͖ɼԾఆ 1 ͷ΋ͱͰɼ

Theorem 4
1. Ծఆ 3 ͷ΋ͱͰ {rT−1 (θ̂(q) − θ∗ )}T ͸ L∞− -༗քɽ2. ೚ҙͷ L > 0 ʹର͠ɼఆ਺ CL ͕ଘࡏ͠ɼ

rT−1 (θ̂(q) − θ∗ )J 1 − G{rT−1 (θ̃ − θ∗ )} →p 0.

(q)

P [θ̂J 0 = 0] ≥ 1 − CL rT2L

ಛʹɼԾఆ 2 ͷ΋ͱͰɼ1͕೚ҙͷ T > 0 ʹରͯ͠੒Γཱͭɽ

rT−1 (θ̂(q) − θ∗ )J 1 →ds GG− 2 ζ ∼ MNp0 (0, (GJ 11 )−1 ).

ླ໦ɹ୓ւ (౦‫ژ‬େֶେֶӃ਺ཧՊֶ‫ڀݚ‬Պ)γϯϙδ΢Ϝʮߴ࣍‫ݩ‬ෳࡶσʔλͷ౷‫ܭ‬ϞσϦϯάʯ2019 ೥ 8 ݄ 30 ೔

5 / 10ླ໦ɹ୓ւ (౦‫ژ‬େֶେֶӃ਺ཧՊֶ‫ڀݚ‬Պ)γϯϙδ΢Ϝʮߴ࣍‫ݩ‬ෳࡶσʔλͷ౷‫ܭ‬ϞσϦϯάʯ2019 ೥ 8 ݄ 30 ೔

6 / 10penalized LSA estimator

P-O estimatorsimulation

definition

penalized LSA estimator

P-O estimatorsimulation

cox model

P-O ਪఆྔ
1. rT−1 -ҰகੑΛ࣋ͭॳ‫ظ‬ਪఆྔ θ̃ Λ༻ҙ͢Δɽ(q)

2. θ̃ Λ༻͍ͯɼ୯Ґߦྻ Ip Λ܎਺ߦྻͱ͢Δ penalized LSA ਪఆྔ θ̂Ip Λ‫ٻ‬ΊΔɿ(q)θ̂Ip ∈ argmin
θ∈Θ̄

p &
#
j=1

'
(θj − θ̃j )2 + κjT |θj |q .

(2.1)

࣍ͷΑ͏ͳ intensity λ(t, θ) Λ΋ͭ point process N Λߟ͑Δɽ
&#'λ(t, θ) = exp
θj Xtj(3.1)

j∈J

͜͜ͰɼX = (Xtj )t∈[0,T ],j=1,...,p ͸ࠨ࿈ଓͳ p-࣍‫཰֬ݩ‬աఔͱ͢ΔɽϞσϧ (3.1) ʹରͯ͠ɼٖࣅ໬౓ղੳΛߟ͑Δɽٖࣅର਺໬౓ؔ਺Λ#( T#( T αℓT (θ) =log (λα (t, θ))dNtα −λ (t, θ)dt.

3. ৽͍͠ϩεؔ਺ LT (θ) Λ༻͍ͯɼP-O ਪఆྔ θ̌ ΛҎԼͰఆٛ͢Δɽθ̌ ∈ argmin LT (θ)
θ∈Θ̂

α∈I

(q)
͜͜ͰɼΘ̂ = {θ ∈ Θ; θj = 0, j ∈ Jˆ0 }ɼJˆ0 = {j = 1, ..., p; θ̂Ip ,j = 0} ͱͨ͠ɽ

0

α∈I

(3.2)

0

ͱ͓͘ͱɼLT (θ) = −ℓT (θ) ͸ϩεؔ਺ͱΈͳͤΔɽ·ͨɼٖࣅ࠷໬ਪఆྔ θ̃ Λ ℓT (θ) Λ࠷େԽ͢Δ΋ͷͱͯ͠༩͑Δͱɼθ̃ ͸Ծఆ 1-3 Λຬͨ͢ɽP-O ਪఆྔ͸ɼม਺બ୒ʹ͍ͭͯ͸ penalized LSA ਪఆྔͱಉ౳ͷੑ࣭Λ࣋ͭҰํɼ༗ޮਪఆྔͱͳ͍ͬͯΔɽ͞Βʹɼ܎਺ߦྻΛ୯Ґߦྻʹͨ͜͠ͱʹΑΓɼ(2.1) Λ‫ͯݟ‬෼͔ΔΑ͏ʹɼ໨తؔ਺ͷ࠷దԽ໰୊͕ 1 ࣍‫ݩ‬ͷ࠷దԽ໰୊ʹ‫ؼ‬ண͞Ε͍ͯΔɽҰൠʹม਺બ୒ͷ໰୊͸ߴ࣍‫ͳʹݩ‬Δ͜ͱ͕ଟ͍ͨΊɼ͜Ε͸ඇৗʹ༗༻Ͱ͋Δɽ
ླ໦ɹ୓ւ (౦‫ژ‬େֶେֶӃ਺ཧՊֶ‫ڀݚ‬Պ)γϯϙδ΢Ϝʮߴ࣍‫ݩ‬ෳࡶσʔλͷ౷‫ܭ‬ϞσϦϯάʯ

penalized LSA estimator

2019 ೥ 8 ݄ 30 ೔

P-O estimator

7 / 10simulation

resultླ໦ɹ୓ւ (౦‫ژ‬େֶେֶӃ਺ཧՊֶ‫ڀݚ‬Պ)γϯϙδ΢Ϝʮߴ࣍‫ݩ‬ෳࡶσʔλͷ౷‫ܭ‬ϞσϦϯάʯ

penalized LSA estimator

2019 ೥ 8 ݄ 30 ೔

P-O estimator

8 / 10simulation

resultTable: Results of the variable selection under T = 50, 100, 200, 400.

covariate process X = (Xt )t∈[0,T ] ͸ 20-࣍‫ ݩ‬OU աఔͰҎԼͷ SDE Λຬͨ͢΋ͷΛߟ͑Δɽ
dXti = −ai Xti dt + 0.4dWti , X0 = 0, t ∈ [0, T ].

%(p-LSA)% (unified LASSO)% (Bridge type)͜͜Ͱɼai , i = 1, · · · , 20 ͸

a1 = a6 = a11 = a16 = 0.15,
a3 = a8 = a13 = a18 = 0.25,

true
initiala4 = a9 = a14 = a19 = 0.3,

p-LSA

a5 = a10 = a15 = a20 = 0.35

P-O

ͱ͠ɼW = (W i )i=1,··· ,20 ͸ 20-࣍‫ݩ‬ඪ४ Wiener աఔͱ͢Δɽσʔλ͸ intensity λ(t, θ ∗ )Λ΋ͭ point process N = (Nt )t∈[0,T ] ͷαϯϓϧύεͱ͢Δɽ·ͨɼਅ஋ θ∗ ͸

θ = [2, −1, 1, −0.5, −1.5, 1.5, 0.5, 0.75, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]′

T = 200
96.8
47.1
52.3

T = 400
99.9
79.6
80.3

9 / 10θ2-1-0.9936θ3
10.9941

θ4-0.5-0.5003

θ5-1.5-1.4909θ6
1.5
1.4978θ7
0.5
0.5009θ8
0.75
0.7490θ9
00.0001θ10
0-0.0002(0.0830)(0.0858)(0.0889)(0.0962)(0.0745)(0.0851)(0.0880)(0.0908)(0.0974)

1.9918-0.9872

0.9877-0.4758-1.4877

1.49460.47500.7383-0.0001-0.0008(0.0728)(0.0840)(0.0868)(0.1035)(0.0965)(0.0748)(0.1049)(0.0904)(0.0265)(0.0340)

1.9959-0.9971

0.9964-0.4979-1.49311.49980.49460.7523-0.0003-0.0007(0.0682)(0.0713)(0.0892)(0.0823)(0.0602)(0.0915)(0.0730)(0.0228)(0.0307)

100.0

100.0

100.0

98.8

100.0

100.0

98.4

100.0

99.4

99.2

true
initialθ1100.0014θ12
0-0.0038θ13
00.0007

θ14
00.0026θ15
00.0063θ16
0-0.0015θ17
00.0042θ18
00.0013θ19
0-0.0048θ20
0-0.0014(0.0760)(0.0794)(0.0872)(0.0947)(0.0971)(0.0701)(0.0914)(0.0941)%(p-LSA)2019 ೥ 8 ݄ 30 ೔

θ12
1.9938(0.0722)(0.0565)

P-Oγϯϙδ΢Ϝʮߴ࣍‫ݩ‬ෳࡶσʔλͷ౷‫ܭ‬ϞσϦϯάʯT = 100
70.4
17.4
21.9

%(p-LSA)

p-LSA

ͱ͢Δɽ

ླ໦ɹ୓ւ (౦‫ژ‬େֶେֶӃ਺ཧՊֶ‫ڀݚ‬Պ)T = 50
32.1
5.9
8.9

Table: The summary of results for the simulation under T = 200.

a2 = a7 = a12 = a17 = 0.2,∗(γ, r, q)(1, 1.2, 0.3)
(1, 1.2, 1)(0, 1, 0.3)

(0.0802)(0.0896)0.00000.00030.00070.0002-0.0004

0.00000.0011-0.00030.0013(0.0000)(0.0185)(0.0252)(0.0240)(0.0248)(0.0000)(0.0194)(00239)(0.0247)0.0000-0.00030.00070.0000-0.00030.00000.0007-0.0004

0.00100.0004(0.0000)(0.0152)(0.0221)(0.0188)(0.0259)(0.0000)(0.0134)(0.0214)(0.0199)(0.0213)

100.0

99.7

99.4

99.4

99.5

100.0

99.7

99.5

99.7

99.5

ླ໦ɹ୓ւ (౦‫ژ‬େֶେֶӃ਺ཧՊֶ‫ڀݚ‬Պ)γϯϙδ΢Ϝʮߴ࣍‫ݩ‬ෳࡶσʔλͷ౷‫ܭ‬ϞσϦϯάʯ

0.0004(0.0208)2019 ೥ 8 ݄ 30 ೔

10 / 10ܽଌσʔλΛ༻͍ͨϑΟογϟʔείΞϦϯά๏ؔ੢େֶ ߴҪ‫ܒ‬ೋຊൃදͰ͸ɼܽଌσʔλʹ΋ͱ͍ͮͯ࠷໬ਪఆ஋Λ‫͢ࢉܭ‬Δख๏ͱͯ͠ෆ‫׬‬શσʔλͷϑΟογϟʔείΞϦϯά๏ (Incomplete-data Fisher Scoring; IFS) Λఏএ͢ΔɽIFS ͷಋग़ํ๏ɼऩଋͷੑ࣭ɼ࣮ͦͯ͠ࡍͷద༻ྫΛࣔ͢ɽैདྷɼܽଌσʔλʹ΋ͱ͍ͮͯ࠷໬ਪఆ஋Λ‫͢ࢉܭ‬Δʹ͸ɼEM ΞϧΰϦζϜ͕ඪ४తʹ༻͍ΒΕ͖ͯͨɽ‫Ͱࡏݱ‬͸ɼEM ΞϧΰϦζϜ͸୯ʹܽଌσʔλʹ΋ͱ͍ͮͯύϥϝʔλΛਪఆ͢Δ͚ͩͰͳ͘ɼજࡏม਺Ϟσϧͷύϥϝʔλͷਪఆ΍Ұൠͷਪఆؔ਺ͷύϥϝʔλਪఆʹ΋ར༻͞Ε͖͍ͯͯΔɽ‫ࡏݱ‬ͷͱ͜Ζɼܽଌσʔλʹ΋ͱ͍ͮͯύϥϝʔλΛਪఆ͢Δࡍʹ͸ɼEM ΞϧΰϦζϜ͸ࣄ্࣮ͷୈҰબ୒ީิͱͳ͍ͬͯΔɽEM ΞϧΰϦζϜͷར఺͸ɼͦͷߋ৽ࣜͷಋग़͕͠΍͍͢͜ͱɼͦͯ͠ߋ৽͕ࣜ୯७ͳ
‫ͳʹܗ‬Δ͜ͱͰ͋Δɽ͜Ε͸࣮૷ͷ͠΍͢͞΍ɼߋ৽ࣜͷղऍͷ͠΍͢͞ʹ௚͍݁ͯ͠
Δɽ͔͠͠ɼͦΕ͸ EM ΞϧΰϦζϜ͕࠷΋༏ΕͯͨΞϧΰϦζϜͰ͋Δ͜ͱΛҙຯͯ͠͸͍ͳ͍ɽEM ΞϧΰϦζϜ͸ɼैདྷ͔Βऩଋ͕஗͍͜ͱ͕ࢦఠ͞Ε͖ͯͨɽ͜͜Ͱ‫ݴ‬͏ऩଋͷ஗͞ͱ͸ɼऩଋ·Ͱͷ‫܁‬Γฦ͠ճ਺ͷ͜ͱͰ͋ΔɽEM ΞϧΰϦζϜ͸ɼଞͷΞϧΰϦζϜʹൺ΂ͯऩଋ·Ͱʹཁ͢Δ‫܁‬Γฦ͠ͷ‫ࢉܭ‬ճ਺͕ଟ͘ͳΔ܏޲ʹ͋Δɽ͜ͷΑ͏ͳ໰୊఺Λղܾ͢΂͘ɼ༷ʑͳՃ଎๏͕ࢼΈΒΕ͖ͯͨɽྫ͑͹ɼdoubling step ๏΍ɼ
Aitken Ճ଎ͳͲͰ͋ΔɽҰൠʹ͜͏͍ͬͨํ๏͸ɼຊ࣭తʹ͸໨తؔ਺ͷೋ֊ඍ෼ͷ৘ใΛ࢖͏͜ͱʹͳΔɽͦΕΏ͑ɼ࠷ॳ͔Β໨తؔ਺ͷೋ֊ඍ෼͋Δ͍͸ͦΕʹ૬Ԡ͢ΔΑ
͏ͳ৘ใΛ࢖͏͜ͱ͕߹ཧతͰ͋Ζ͏ɽͦ͏͍ͬͨߟ͑ํʹ΋ͱ͍ͮͯɼEM ΞϧΰϦζ
Ϝͷ৘ใͱͦͷଞͷ৘ใ (ྫ͑͹ޯ഑) Λར༻͢Δ hybrid EM ΞϧΰϦζϜ΍ɼEM ͷΞΠσΟΞΛ༻͍ͳ͕Β΋ EM ΞϧΰϦζϜΛ࢖Θͳ͍ߋ৽ࣜ (ٙχϡʔτϯ๏΍‫ݧܦ‬తͳϑΟογϟʔ৘ใߦྻΛ༻͍ͨϑΟογϟʔείΞϦϯά๏) ͕։ൃ͞Ε͖ͯͨɽຊ‫Ͱڀݚ‬ఏএ͢Δ IFS ͸ɼલऀͱ‫ऀޙ‬ͷதؒʹҐஔ෇͚Δ͜ͱ͕Ͱ͖ΔɽIFS ͸εςοϓ෯Λௐ੔͢Δ͜ͱʹΑΔՃ଎๏Ͱ͋Δɽ۩ମతʹ͸ɼύϥϝʔλ α ͷ‫ࡏݱ‬ͷਪఆ஋ α(t)͔Β৽ͨͳਪఆ஋ α(t+1) ΁ͷߋ৽Λα(t+1) = α(t) + s"−11!
Jcom (α(t) )∇ℓobs (α(t) )nͰߦͳ͏ɽ͜͜Ͱɼn ͸αϯϓϧαΠζɼs ͸εςοϓ෯ɼJcom (α) ͸‫׬‬શσʔλͷϑΟο
γϟʔ৘ใߦྻͰ͋ΔɽIFS ͷಛච͢΂͖ੑ࣭͸ɼྫ͑͹ਖ਼‫ن‬෼෍ͷ৔߹ʹ͸ s = 1 ͱ͢Δͱɼ(ಉҰͷॳ‫ظ‬஋Ͱ͋Ε͹)EM ΞϧΰϦζϜͱશ͘ಉҰͷਪఆ஋ͷྻΛੜ੒͢Δͱ͍
͏͜ͱͰ͋Δɽ͜ͷ͜ͱ͸ɼs > 1 ͱ͢Δ͜ͱʹΑΓɼEM ΞϧΰϦζϜΑΓ΋଎͍ऩଋ͕‫ࠐݟ‬ΊΔͱ͍͏͜ͱͰ΋͋Δɽͦͯ͠ɼεςοϓ෯ s Λద੾ʹௐઅ͢Δ͜ͱʹΑΓɼΑΓߴ଎ͳऩଋΛ࣮‫͖Ͱݱ‬ΔͰ͋Ζ͏ɽ
EM ͕͍࣋ͬͯͨ໨తؔ਺ͷ୯ௐ૿Ճੑʹ͍ͭͯ͸ IFS ΋༗͍ͯ͠ΔɽIFS ʹΑΔ໬౓ͷߋ৽͸ɼ

ℓobs (α(t+1) ) − ℓobs (α(t) ) ≈ ∇ℓobs (α(t) )′ (α(t+1) − α(t) )
!
"−1
s
=
∇ℓobs (α(t) )′ Jcom (α(t) )∇ℓobs (α(t) )
n
1

0.10
0.05
0.00

Difference

AIFS
IFS
EM

0

500

1000

1500

2000

Iteration

ਤ 1: ऩଋͷաఔͰ༩͑ΒΕΔɽ࠷‫ޙ‬ͷ߲ͷ Jcom (α) ͕ਖ਼ఆ஋Ͱ͋Ε͹ɼℓobs (α(t+1) ) − ℓobs (α(t) ) > 0 ͸໌Β͔Ͱ͋ΔɽҰൠతͳ౷‫ܭ‬Ϟσϧʹ͓͍ͯɼϑΟογϟʔ৘ใߦྻ͕ਖ਼ఆͰ͋Δ͜ͱΛԾఆ͢Δ͜ͱ͸੍ܾͯ͠‫ݶ‬తͰ͸ͳ͍Ͱ͋Ζ͏ɽैͬͯɼҰൠతͳ౷‫ܭ‬Ϟσϧʹରͯ͠ (εςοϓ෯ s ͕ద੾Ͱ͋Δ‫ݶ‬Γ͸)ɼ໨తؔ਺ͷ୯ௐ૿ՃੑΛ୲อͰ͖Δɽຊ‫Ͱڀݚ‬͸ɼ͜ͷ௚‫ʹײ‬ཧ࿦తͳཪ͚ͮΛ༩͑ͨɽ۩ମతʹ͸ɼ(1) IFS Λ୯ௐʹ૿Ճͤ͞Δεςοϓ෯ (s > 0) ͕ඞͣଘࡏ͢Δ͜ͱɼ(2) IFS ͷऩଋͷ଎͞͸ ∇2 ℓobs (α(t) )Jcom (α(t) )−1ͷ‫ݻ‬༗஋ʹґଘ͢Δ͜ͱɼͰ͋Δɽ࠷‫ʹޙ‬ɼIFS Λ୯มྔࠞ߹ϙΞιϯ෼෍ͷύϥϝʔλਪఆʹద༻ͨ݁͠ՌΛख୹ʹใࠂ͢Δɽ1910 ೥͔Β 1912 ೥ͷ London Times ࢴͷࢮ๢ਓ਺ͷσʔλ (Hasselblad, 1969;
Titterington, 1984) ʹ͸ɼೋͭͷϙΞιϯ෼෍ͷࠞ߹෼෍͕ྑ͘ద߹͢Δ͜ͱ͕஌ΒΕ͍ͯΔ (Lange, 1995)ɽ͜ͷ෼෍ͷύϥϝʔλਪఆʹ IFS Λ༻͍ͯΈΔɽਤ 1 ͸ɼEM ΞϧΰϦζϜʹΑΔਪఆաఔɼεςοϓ෯Λ 2 ‫ݻ‬ఆͨ͠ IFS ͷਪఆաఔɼ࠷దͳεςοϓ෯ʹ
ͨ͠ IFS ͷਪఆաఔΛ͍ࣔͯ͠Δ (ͦΕͧΕਤதͰ͸ EMɼIFSɼAIFS ͱද‫)ه‬ɽԣ࣠͸‫܁‬Γฦ͠ճ਺Ͱ͋Γɼॎ࣠͸ऩଋઌͷ໬౓ͱͦͷ࣌఺Ͱͷ໬౓ͷࠩͰ͋ΔɽͲͷํ๏΋୯ௐʹऩଋ͍ͯ͠Δ͜ͱ͕෼͔Δɽऩଋ·ͰʹɼEM ΞϧΰϦζϜ͸ 2000 ճҎ্ɼεςοϓ෯Λ‫ݻ‬ఆͨ͠ IFS ͸ 1500 ճఔ౓ɼεςοϓ෯Λௐ੔ͨ͠ IFS ͸ 200 ճఔ౓ͷ‫܁‬Γฦ͠‫ܭ‬ࢉ͕ඞཁͱͳΔɽεςοϓ෯Λௐ੔ͨ͠ IFS ͸ EM ΞϧΰϦζϜͷ 1/10 ΄Ͳͷ‫܁‬Γฦ͠ճ਺ͱͳͬͯɼద੾ͳεςοϓ෯Λ༻͍Δ͜ͱͰՃ଎Ͱ͖Δ͜ͱ͕෼͔Δɽ

2

Spectral Embedded Deep Clustering
Yuichiro WADA1 , Shugo MIYAMOTO3 , Takumi NAKAGAWA2 , Léo Andéol4,5 , Wataru
KUMAGAI5 , and Takafumi KANAMORI2,5
1 Nagoya University, Furocho, Chikusaku, Nagoya 464-8601 Japan

2 Tokyo Institute of Technology, 2-12-1 Ookayama, Meguro-ku, Tokyo 152-8552 Japan
3 The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656 Japan
4 Sorbonne Université, 4 place Jussieu 75005 Paris, France

5 RIKEN AIP, Nihonbashi 1-chome Mitsui Building, 15th floor,1-4-1 Nihonbashi, Chuo-ku,

Tokyo 103-0027 Japan

1

Abstract

We propose a new clustering method based on a deep neural network. Given an unlabeled dataset and the numberof clusters, our method directly groups the dataset into the given number clusters in the original space. We use aconditional discrete probability distribution defined by a deep neural network as a statistical model. Our strategy
is, first to estimate the cluster labels of unlabeled data points selected from high density region, and then to conduct
semi-supervised learning to train the model by using the estimated cluster labels and the remaining unlabeleddata points. Lastly, by using the trained model, we obtain the estimated cluster labels of all given unlabeled data
points. The advantage of our method is that it does not require key condition. Existing clustering methods withdeep neural networks based assumed that the cluster-balance of given dataset is uniform. Moreover, it also can beapplied to various data domains as long as the data is expressed by a feature vector. In addition, it was observedthat our method was robust against outliers. Therefore, the proposed method is expected to averagely performbetter than previous methods. We conducted numerical experiments on five commonly used datasets to confirm
the effectiveness of the proposed method.

2

Proposed Methodn (x 2 R D ) and the number of clusters C, our proposed deep clustering
Given an unlabeled dataset X = {xi }i=1i
named SEDC (Spectral Embedded Deep Clustering) groups X into C clusters. Since this grouping is achieved byobtaining the estimated cluster labels of X, our goal can be replaced by estimating the cluster labels up to permutationof labels. In SEDC, the estimated cluster label of each xi 2 X is defined by argmax j=1,...,C p✓ ⇤ (y = j |xi ), where ✓ ⇤is the trained set of parameters. The classifier p✓ (y|x) is parameterized by a fully connected deep neural network.
The training scheme of this classifier is as follows: we firstly only estimate the cluster labels of selected unlabeleddata points by using only X (this part is done by SGSC (Selective Geodesic Spectral Clustering) algorithm.), andthen conduct semi-supervised learning to train the classifier. Regarding with this semi-supervised learning, we
use the estimated cluster labels of selected unlabeled data points and the remaining unlabeled data points, which
are treated as the given true cluster labels and unlabeled data points respectively. The objective function of thissemi-supervised learning is as follows:

R VAT (✓; X) +

1

h
’

h i=1

⇥
⇤
KL p✓ (y|x(i) )kq(i) +2 H (Y|X) ,(1)

where the first, second and third terms are VAT loss, pseudo empirical loss with estimated cluster probabilities (q(i) ) of selected data points based on KL divergence [1] and conditional entropy loss [1], respectively. 1 and 2
are hyperparameters that range over positive numbers. y and h mean the cluster label and the number of selected
data points, respectively. Eq.(1) will be minimized.
1

Table 1: The mean clustering accuracy (ACC) of Eq.(2) and standard deviation are shown. Five popular clustering
methods and our proposed method were tested on five datasets. MNIST and Reuters are real-world datasets. FC(Four Clusters), TM (Two Monns) and TR (Three Rings) are artificial datasets. For each method, Average means
the averaged ACC over the five datasets. The experiments were conducted seven times on each pair of method anddataset.
MethodMNIST
ReutersFC
TM
TR
Average
10k
k-means [4]
SC [5]IMSAT [2]DEC [7]
SpectralNet [6]
SEDC

3

0.53
0.72
0.98
0.84
0.83
0.89

0.53(0.04)
0.62(0.03)
0.71(0.05)
0.72(0.05)
0.67(0.03)
0.73(0.05)

0.60(0.05)
0.80(0.04)
0.70(0.04)
0.72(0.04)
0.79(0.03)
0.95(0.03)

0.64(0.04)
0.85(0.03)
0.66(0.05)
0.67(0.03)
0.87(0.02)
0.96(0.02)

0.35(0.03)
0.96(0.03)
0.34(0.01)
0.48(0.04)
0.99(0.01)
0.99(0.00)

0.53
0.79
0.68
0.69
0.83
0.90

Numerical Experiment and Conclusion

Table 1 shows the clustering accuracy of each methods. Our method is SEDC. The evaluation metric is as follows:Õn
1[yi = ⌧( ŷi )]ACC = max i=1
,
(2)
⌧
n
n and { ŷ } n be its true cluster label set and estimated cluster label set, respectively. The number ofwhere {yi }i=1
i i=1data points is denoted by n. ⌧ ranges over all permutations of cluster labels, and 1[· · · ] is the indicator function.
The optimal assignment of ⌧ can be computed using the Kuhn-Munkres algorithm [3].
With respect to the conclusion, we propose a deep clustering method named SEDC. Given an unlabeleddataset and the number of clusters, the method groups the dataset into the given number clusters. Regarding
its advantages, it does not require an additional condition except two fundamental assumptions: smoothness andmanifolds assumptions. In addition, SEDC also can be applied to various data domains since it does not havepreferred data domains, as long as raw data is transformed to feature vectors. Furthermore, the performance of
SEDC can be robust against existence of outliers. According to these advantages, our proposed method can be
expected to averagely perform better than previous deep clustering methods.Therefore, we think our method canbe a competitive candidate for users in some practical clustering scenarios where prior knowledge of the given
unlabeled dataset is limited.

References[1] Thomas M Cover and Joy A Thomas. Elements of information theory. John Wiley & Sons, 2012.
[2] Weihua Hu, Takeru Miyato, Seiya Tokui, Eiichi Matsumoto, and Masashi Sugiyama. Learning discrete
representations via information maximizing self-augmented training. In International Conference on Machine
Learning, pages 1558–1567, 2017.
[3] Harold W Kuhn. The hungarian method for the assignment problem. Naval research logistics quarterly,2(1-2):83–97, 1955.
[4] James MacQueen et al. Some methods for classification and analysis of multivariate observations. In
Proceedings of the fifth Berkeley symposium on mathematical statistics and probability, volume 1, pages
281–297. Oakland, CA, USA, 1967.
[5] Andrew Y Ng, Michael I Jordan, and Yair Weiss. On spectral clustering: Analysis and an algorithm. In
Advances in neural information processing systems, pages 849–856, 2002.
[6] Uri Shaham, Kelly Stanton, Henry Li, Boaz Nadler, Ronen Basri, and Yuval Kluger. Spectralnet: Spectralclustering using deep neural networks. arXiv preprint arXiv:1801.01587, 2018.
[7] Junyuan Xie, Ross Girshick, and Ali Farhadi. Unsupervised deep embedding for clustering analysis. In
International conference on machine learning, pages 478–487, 2016.
2Պ‫ݚ‬අγϯϙδ΢Ϝʮߴ࣍‫ݩ‬ෳࡶσʔλͷ౷‫ܭ‬ϞσϦϯάʯ2019 ೥ 8 ݄ 30 ೔ൃදʮ‫ۭؒ཭ڑ‬ͷม‫ͱܗ‬ຒΊࠐΈΛ༻͍ͨσʔλղੳʯใࠂॻ
‫ܚ‬ጯٛक़େֶɹখྛɹ‫ܠ‬

1

ൃද֓ཁ

ຊൃදͰ͸·ͣɼϦʔϚϯଟ༷ମ΍‫্ۭؒ཭ڑ‬ͷ Fréchet ฏ‫ۉ‬ɼextrinsic ฏ‫͓ۉ‬Αͼ Fréchetؔ਺ʹؔ͢Δ‫ط‬ଘ‫ڀݚ‬Λ֓આͨ͠ɽ͜ΕΒ͸ϕΫτϧσʔλͰզʑ͕௨ৗ༻͍͍ͯΔʮฏ‫ۉ‬ʯΛɼઢ‫ͱܗ‬͸‫ݶ‬Βͳ͍‫ۭؒ཭ڑ‬ʢྫ͑͹‫ٿ‬໘ʣ্ʹҰൠԽ͢Δඞཁ͕͋Δͱ͖ʹࣗવʹಋ͔ΕΔ౷‫͋Ͱྔܭ‬Δɽ࣍ʹɼ͜ΕΒͷ‫ڀݚ‬෼໺ͷల։ͷҰͭͱͯ͠ɼ[1] ʹ͓͍ͯஶऀΒ͕ఏҊͨ͠‫ۭؒ཭ڑ‬ͷม‫ܗ‬Λ༻͍ͨσʔλղੳख๏͓Αͼͦͷཧ࿦Λ؆୯ʹ঺հͨ͠ɽ

2Fréchet ฏ‫ۉ‬ʢintrinsic ฏ‫ۉ‬ʣͱ Extrinsic ฏ‫ۉ‬·ͣɼ‫্ۭؒ཭ڑ‬ͷ֬཰෼෍΋͘͠͸ඪຊ෼෍ʹର͢Δ Fréchet ฏ‫͓ۉ‬Αͼ extrinsic ฏ‫͍ͯͭʹۉ‬આ໌ͨ͠ɽ·ͨ઴ۙཧ࿦Λ༻͍Δ͜ͱʹΑΓɼFréchet ฏ‫ʹۉ‬ର͢ΔԾઆ‫ݕ‬ఆ
ʢH0 : µ = µ0 vs H1 : µ ̸= µ0 ʣ΍̎ඪຊͷ Fréchet ฏ‫ۉ‬ͷࠩҟʹؔ͢Δ‫ݕ‬ఆɼ΋͘͠͸
ϕ(µ) ͷ৴པ۠ؒͷߏ੒ͱ͍ͬͨ௨ৗͷϢʔΫϦου্ۭؒͷ౷‫ֶܭ‬తख๏ΛϦʔϚϯଟ༷ମ্ʹҰൠԽ͢Δ͜ͱ͕ՄೳͱͳΔ͜ͱΛઆ໌ͨ͠ɽ
࣍ʹɼFréchet ฏ‫ۉ‬Ҏ֎ͷ‫্ۭؒ཭ڑ‬ͷʮฏ‫ۉ‬ʯͷީิͱͯ͠ɼFréchet ฏ‫͘ͱۉ‬Β΂ͯҰൠʹࢉग़͕༰қͳ extrinsic ฏ‫ۉ‬Λ঺հͨ͠ɽextrinsic ฏ‫ۉ‬͸ɼϢʔΫϦουۭؒʹຒΊࠐ·Εͨ‫ʹۭؒ཭ڑ‬ରͯ͠ɼϢʔΫϦου‫཭ڑ‬Λ༻͍Δ Fréchet ฏ‫ͯ͠ͱۉ‬ఆٛ͞ΕΔɽൃ
දʹ͓͍ͯ͸ extrinsic ฏ‫ۉ‬ͷࣹӨఆཧΛ঺հͨ͠‫ޙ‬ɼͦͷࢉग़͓Αͼ઴ۙతٞ࿦͕ Fréchetฏ‫͘ʹۉ‬Β΂ͯ༰қʹͳΔཧ༝Λઆ໌ͨ͠ɽ
extrinsic ฏ‫ۉ‬͸ࣹӨఆཧΛ༻͍Δͱࢉग़΍઴ۙཧ࿦͕ฏқԽ͞ΕΔͱ͍͏େ͖ͳϝϦοτ͕͋Δ͕ɼͦ΋ͦ΋‫͕ۭؒ཭ڑ‬ϢʔΫϦουۭؒʹຒΊࠐΊͳ͍Α͏ͳ৔߹ʹ͸༻͍Δ͜
ͱ͕Ͱ͖ͳ͍ɽಛʹɼࠓճ঺հ͢Δ‫Ͱڀݚ‬ఏҊ͢ΔΑ͏ʹσʔλۭؒͷ‫཭ڑ‬Λม‫͢׵‬Δࡍʹ͸ɼϢʔΫϦουۭؒ΁ͷຒΊࠐΈ͸͙͢ʹ੒Γཱͨͳ͘ͳΔɽ͜ͷ͜ͱΛ‫ݟ‬ΔͨΊʹɼ‫ڑ‬཭ۭؒͷϢʔΫϦουۭؒ΁ͷຒΊࠐΈՄೳੑʹ͍ͭͯɼSchoenberg ʹΑΔΑ͘஌ΒΕͨ݁ՌΛ঺հͨ͠ɽ

3Fréchet ฏ‫ۉ‬ͷҰҙੑͱ CAT(k) ۭؒҰൠʹ Fréchet ฏ‫߹ूۉ‬΋ extrinsic ฏ‫߹ूۉ‬΋αΠζ͕̎Ҏ্ʹͳΓಘΔͨΊɼฏ‫ۉ‬ͷҰҙੑ͸੒ཱ͠ͳ͍ɽextrinsic ฏ‫ۉ‬ͷ৔߹͸ຒΊࠐ·ΕͨϢʔΫϦουۭؒͰͷ Fréchet ฏ‫͕ۉ‬Ұҙʹఆ·Δ͜ͱ͔ΒɼҰҙੑ͸௚ަࣹӨͷҰҙੑɼͻ͍ͯ͸ϢʔΫϦουۭؒ΁ͷຒΊࠐ·Εํʹґଘ͢ΔɽҰํɼFréchet ฏ‫ۉ‬ͷҰҙੑ͸σʔλۭؒࣗ਎ͷ‫͔཭ڑ‬Βఆ·Δ‫ۂ‬཰ʹΑΓಛ௃͚ͮ͞ΕΔɽൃදʹ͓͍ͯ͸ɼ·ͣσʔλۭ͕ؒϦʔϚϯଟ༷ମͰ͋Δͱ͖ͷ‫ ͱ཰ۂ‬Fréchet ฏ‫ۉ‬ͷҰҙੑͷؔ܎Λ؆୯ʹ঺հͨ͠‫ޙ‬ɼΑΓҰൠʹ Length space ্ͷ‫ۂ‬཰ʢCAT(k) ੑʣͱ Fréchet ฏ‫ۉ‬ͷҰҙੑʹ͍ͭͯ஌ΒΕ͍ͯΔ݁Ռʹ͍ͭͯ΋঺հͨ͠ɽ

Ұํɼ‫཭ڑ‬άϥϑɼଟ໘ମతෳମɼਲ਼ͳͲ͸Ұൠʹඍ෼Մೳͳଟ༷ମͰ͸ͳ͘ɼ͜Ε·Ͱͷٞ࿦Λ౰ͯ͸ΊΔ͜ͱ͕Ͱ͖ͳ͍ɽͦ͜Ͱɼ্‫ه‬ͷΑ͏ͳ‫ز‬Կֶతର৅΍ϦʔϚϯଟ༷ମΛ‫ؚ‬Ή Length space ʹཧ࿦ΛҰൠԽͨ͠‫ط‬ଘ‫ڀݚ‬Λ঺հͨ͠ɽ͞Βʹɼஅ໘‫ ཰ۂ‬k ҎԼͷϦʔϚϯଟ༷ମʹରԠ͢Δ CAT(k) ۭؒͷఆٛͱੑ࣭Λ঺հͨ͠ɽ

4β ‫཭ڑ‬ม‫׵‬࣍ʹɼ[1] ʹ͓͍ͯஶऀΒ͕ఏҊͨ͠σʔλۭؒͷ‫཭ڑ‬ม‫׵‬ख๏ʹ͍ͭͯ؆୯ʹઆ໌ͨ͠ɽ
‫཭ڑ‬ͷม‫׵‬͸ɼ‫཰ۂ‬ɼಛʹ CAT(k) ੑʹண໨ͯ͠ఏҊ͞Εͨ α ‫཭ڑ‬ม‫ ͱ׵‬β ‫཭ڑ‬ม‫ͱ׵‬Α͹ΕΔ̎ஈ֊ͷม‫׵‬ͷ߹੒ͱͯ͠ఆٛ͞ΕΔɽα ‫཭ڑ‬ม‫͍ͯͭʹ׵‬͸؆୯ʹઆ໌ͨ͠‫ޙ‬ɼࠓճͷൃදͷओ୊Ͱ͋Δ β ‫཭ڑ‬ม‫׵‬ͷఆ͓ٛΑͼ‫ྔܭ‬ਲ਼΁ͷຒΊࠐΈͱͷؔ܎Λઆ໌ͨ͠ɽ
·ͣɼҎԼͷΑ͏ʹ‫ ཭ڑ‬d ͱԜؔ਺ gβ (·) Λ߹੒͢Δ͜ͱʹΑΓɼ‫ ཭ڑ‬dβ Λఆٛ͢Δɽdβ (x0 , x1 ) = gβ (d(x0 , x1 )) (β > 0), ɹͨͩ͠!
sin( πz
for 0 ≤ z ≤ β,
2β ),gβ (z) =1,
for z > β"Ͱ͋Δɽβ ͷมԽʹΑΓɼFréchet ؔ਺ F (y) = i dβ (xi , y)γ ͷ‫ۃ‬খ఺ʢKarchar ฏ‫ۉ‬ʣͷ਺Λௐ੔Ͱ͖ΔɽFréchet ؔ਺ͷ‫ۃ‬খ఺ΛΫϥελʔத৺ͷީิͱ͢Δ͜ͱʹΑΓɼΫϥελϦϯάʹԠ༻͢Δ͜ͱ͕Ͱ͖Δɽ͜Ε͸ Fréchet ؔ਺͕ৗʹತؔ਺ͱͳΔϢʔΫϦου‫ڑ‬཭΍ϊϧϜۭؒͰ͸ෆՄೳͳղੳख๏Ͱ͋ΔɽҰํɼβ ‫཭ڑ‬ม‫ʹ׵‬ΑΔ‫཰ۂ‬ͷมԽ͸ɼα ‫཭ڑ‬ม‫ͱ׵‬ҟͳΓ௚઀తʹ͸ٞ࿦Ͱ͖ͳ͍ɽ͜
Ε͸ɼβ ‫཭ڑ‬ม‫׵‬Λ͸͡Ίͱͯ͠ɼLength space(M, d) ΛԜؔ਺ ϕ Λ༻͍ͯม‫ۭ཭ڑͨ͠׵‬ؒ (M, ϕ(d)) ͸Ұൠʹ Length space ͱͳΒͣɼCAT(k) ੑ͕ఆٛͰ͖ͳ͍ͨΊͰ͋Δɽ͔͠͠ β ‫཭ڑ‬ม‫׵‬ಛ༗ͷੑ࣭ͱͯ͠ɼຒΊࠐ·Εͨʮ‫ྔܭ‬ਲ਼ʯͷ‫཰ۂ‬ͷมԽΛධՁ͢Δ͜ͱ͕ՄೳͰ͋ΔɽX ্ͷ‫ ཭ڑ‬dβ ͸ɼຒΊࠐ·Εͨಈ‫ ܘ‬β/2 ͷ‫ྔܭ‬ਲ਼ X̃ Ͱͷ extrinsic ‫ͳͱ཭ڑ‬Δ͜
ͱ͔ΒɼX ্ͷ֬཰෼෍΍ඪຊ෼෍ͷ Fréchet ฏ‫ۉ‬͸ɼ‫ྔܭ‬ਲ਼ʹຒΊࠐΜͩ৔߹ͷ extrinsicฏ‫ͯ͠ͱۉ‬ղऍͰ͖Δɽ·ͨ β ͕খ͘͞ͳΔͱ‫ྔܭ‬ਲ਼ͷ‫͕཰ۂ‬খ͘͞ͳΔ܏޲͕͋Δ͜ͱΛࣔ͢ఆཧΛ঺հͨ͠ɽ
࠷‫ʹޙ‬ɼ·ͱΊͱ՝୊ͱͯ͠ɼϢʔΫϦουۭؒͱ‫ྔܭ‬ਲ਼ͷ extrinsic ฏ‫ۉ‬ͷҧ͍ʹ͍ͭͯ঺հͨ͠ɽ

ࢀߟจ‫ݙ‬[1] Kobayashi, K., Wynn, H.P. (2019), Empirical geodesic graphs and CAT(k) metrics
for data analysis, to appear in Statistics and Computing,https://doi.org/10.1007/s11222-019-09855-3.