学術シンポジウム「統計科学の開拓」
報告書

主催：科学研究費・基盤研究（A)「大規模複雑データの理論と方法論の革新的展開」（課題番号：20H00576）

研究代表者：青嶋 誠（筑波大学）

■学術シンポジウム「統計科学の開拓」実施概要
・開催日時：2022 年 12 月 16 日（金）～12 月 18 日（日）・会場：金沢大学サテライト・プラザ 3 階集会室（金沢市西町三番丁 16 番地）
・開催責任者：星野 伸明（金沢大学）■プログラム12/16(金）
開始

終了
13:00

講演者

演題

13:05 主催者

開催挨拶

江頭健斗(筑波大学理工情報生命学術院)，矢田和善(筑波大学数理

13:10

13:50

13:55

14:35 桃木 光輝（鹿児島大学），吉田 拓真（鹿児島大学）

14:40

15:20 島谷健一郎（統計数理研究所）

15:25

16:05

16:10

16:50

物質系)，青嶋誠(筑波大学数理物質系)

17:35

high-dimensional settings
Mixed effects modeling of clustered extreme
values
スター点過程種市信裕（北海道教育大学札幌校），

スパース性の仮定のもとでの多項分布の適合度検定

関谷祐里（北海道教育大学釧路校）

統計量選択の指針とその応用
Simultaneous estimation of Poisson means in

張元宗・目白大学，
篠崎

信雄

two-way contingency tables under normalized・慶應大学

squarederror loss

星野瞭太（東京理科大学大学院

理工学研究科）
理工学研究科）藤澤健吾（東京理科大学

理工学部）

田畑耕治（東京理科大学

理工学部）

足立 匠 (東京理科大学大学院 理工学研究科 情報科学専攻)、安藤
18:20 宗司(東京理科大学 理工学部 情報科学科)、田畑 耕治(東京理科大

17:40

behaviors in

植物個体群へ応用できる繰り返し可能な非定常クラ

中島亮（東京理科大学大学院
16:55

Hierarchical clustering and its asymptotic

学 理工学部 情報科学科)

正方分割表における点対称性からの隔たりを測る尺
度の推定について
正方分割表におけるモーメントに基づく周辺同等性からの隔たりを測る尺度

12/17(土）
開始

講演者

終了

演題

9:20

10:00 西埜晴久（広島大学大学院人間社会科学研究科）

国民生活基礎調査による日本の所得分布と所得格差

10:05

10:45 国友直人（統計数理研究所特任教授）

操作変数法の理解へ：計量生物と計量経済の邂逅

10:50

11:30 岡野遼（東京大学経済学研究科経済専攻統計学コース博士1年）

多変量分布間の回帰モデル
昼休み

13:00

13:40 永井勇（中京大学）

13:45

14:25

14:30

15:10 田中勝人（一橋大学名誉教授）

GMANOVAモデルでの仮定緩和のさらなる可能性

松内 直輝（神戸大学大学院理学研究科），首藤 信通（神戸大学大 多変量正規母集団における条件付き独立性検定につ学院理学研究科）
Bat-Erdene, A.,

いて
多重積分ブラウン運動の統計学的応用について

Graduate School of Science andEngineering, Iwate University
Kawasaki, S.,
15:15

15:55

Faculty of Science andLi, J.,

School of Computer Science,Supercritical Branching Processes with Poissonian
Offsprings

Chongqing University
Altantsetseg, E.,

On Asymptotic Distribution in Martingale
Convergence ofEngineering, Iwate University

School of the Engineering andApplied Sciences, National University of Mongolia
16:00

16:40 Islam Md Ashraful (大阪大学), 鈴木讓 (大阪大学)

16:45

17:25 山添滉弥(千葉大学融合理工学府)・内藤貫太(千葉大学理学研究院)

17:30

18:10

Forest Construction of Gaussian and Binary
Variables based on WBIC
多次元空間における埋め込み1次元曲線の同時信頼領域

勝又 真（神戸大学大学院理学研究科）

多標本問題における 2-step 単調欠測データの下で

首藤 信通（神戸大学大学院理学研究科）の平均ベクトルの同等性検定の検出力について

12/18(日)
開始

終了

講演者

9:20

10:00 松井宗也（南山大学経営学部）

10:05

10:45 入江薫（東京大学経済学部）

10:50

11:30 前園宜彦（中央大学理工学部）

演題
Subexponentaliy of densities of infinitely divisible
distributions on the whole real line
ガンマ分布の形状パラメータのベイズ推定
カーネル型推定量を利用した推測について
昼休み

13:00

13:40 塚原英敦（成城大学経済学部）

13:45

14:25

On a generalization of Clayton-Oakes model by R.
L. Prentice

助田一晟（東京大学情報理工学系研究科M1），清智也（東京大学 Kendallの順位相関係数を固定した下での最小情報
情報理工学系研究科）

コピュラ

渡邊宏大（千葉大学融合理工学府）・内藤貫太（千葉大学理学研

スピアマンランク行列による主成分分析のロバスト

究院）

ネス

14:30

15:10

15:10

15:15 主催者

閉会挨拶

■報告書目次
講演者

演題頁

江頭健斗(筑波大学理工情報生

Hierarchical

命学術院)，矢田和善(筑波大学

behaviors in high-dimensional settingsclustering

and

its

asymptotic

1

Mixed effects modeling of clustered extreme values

3

植物個体群へ応用できる繰り返し可能な非定常ク

5

数理物質系)，青嶋誠(筑波大学
数理物質系)
桃木光輝（鹿児島大学）, 吉田拓真（鹿児島大学）
島谷健一郎（統計数理研究所）

ラスター点過程
種市信裕（北海道教育大学札幌

スパース性の仮定のもとでの多項分布の適合度検校），関谷祐里（北海道教育大学

定統計量選択の指針とその応用

7釧路校）
張元宗（目白大学），篠崎信雄（慶Simultaneous estimation of Poisson means in two-

應大学）

way contingency tables under normalized squared9

error loss
中島亮（東京理科大学大学院

正方分割表における点対称性からの隔たりを測る

理工学研究科）, 星野瞭太（東京尺度の推定について11

理科大学大学院理工学研究科）,
藤澤健吾（東京理科大学理工学部）, 田畑耕治（東京理科大学理
工学部）
足立匠 (東京理科大学大学院

正方分割表におけるモーメントに基づく周辺同等

理工学研究科), 安藤宗司(東京

性からの隔たりを測る尺度

13

理科大学理工学部), 田畑耕治
(東京理科大学理工学部)
西埜晴久（広島大学大学院人間国民生活基礎調査による日本の所得分布と所得格15

社会科学研究科）

差

国友直人（統計数理研究所特任

操作変数法の理解へ：計量生物と計量経済の邂逅17

多変量分布間の回帰モデル

18

GMANOVA モデルでの仮定緩和のさらなる可能

20

教授）
岡野遼（東京大学経済学研究科
経済専攻）
永井勇（中京大学）

性松内直輝（神戸大学大学院理学

多変量正規母集団における条件付き独立性検定に研究科），首藤信通（神戸大学大

ついて

22

学院理学研究科）
田中勝人（一橋大学名誉教授）

多重積分ブラウン運動の統計学的応用について

24

Bat-Erdene,

A.

On

26

SchoolScienceof(Graduate
andEngineering, Iwate University),

Asymptotic

DistributioninMartingale

Convergence of Supercritical Branching Processes
with Poissonian Offsprings

Kawasaki, S. (Faculty of ScienceandEngineering,

Iwate

University), Li, J. (School ofComputer Science, Chongqing
University),

Altantsetseg,

E.

(School of the Engineering andAppliedSciences,

National

University of Mongolia)
Islam Md Ashraful (大阪大学),

Forest Construction of Gaussian and Binary鈴木讓 (大阪大学)

Variables based on WBIC

山添滉弥(千葉大学融合理工学

多次元空間における埋め込み 1 次元曲線の同時信府), 内藤貫太(千葉大学理学研

頼領域

28
30

究院)
勝又真（神戸大学大学院理学研

多標本問題における 2-step 単調欠測データの下究科）, 首藤信通（神戸大学大学

での平均ベクトルの同等性検定の検出力について

32

院理学研究科）松井宗也（南山大学経営学部）

Subexponentaliy of densities of infinitely divisible

34

distributions on the whole real line
入江薫（東京大学経済学部）

ガンマ分布の形状パラメータのベイズ推定36

前園宜彦（中央大学理工学部）

カーネル型推定量を利用した推測について

37

塚原英敦（成城大学経済学部）

On a generalization of Clayton-Oakes model by R.

39

L. Prentice
助田一晟（東京大学情報理工学

Kendall の順位相関係数を固定した下での最小情系研究科），清智也（東京大学情

報コピュラ

41

報理工学系研究科）
渡邊宏大（千葉大学融合理工学

スピアマンランク行列による主成分分析のロバス

府）, 内藤貫太（千葉大学理学研

トネス

究院）

43

Hierarchical clustering and its asymptotic behaviors
in high-dimensional settings
Kento Egashiraa , Kazuyoshi Yatab , Makoto Aoshimaba

Degree Programs in Pure and Applied Sciences, Graduate School of
Science and Technology, University of TsukubabInstitute of Mathematics, University of Tsukuba

Hierarchical clustering is a methodology to group a set of data by building dendrogram based on a similarity or a dissimilarity between clusters so that data in acluster are similar in the sense of pre-determined linkage function. In hierarchicalclustering, one can observe a process how a cluster is combined or divided throughdendrogram on graphic. Hierarchical clustering has been approved as useful tool
for analysis of gene expression microarray data. In fact, applications of hierarchicalclustering on gene expression microarray data are given by Eisen et al. [4], Perou etal. [8], Bhattacharjee et al. [2], among others. A characteristic of data used in Eisen
et al. [4], Perou et al. [8] and Bhattacharjee et al. [2] is that the number of variables
is much larger than sample size. This type of data represented by gene expression
microarray data is called high-dimension, low-sample-size (HDLSS) data. Substantial work about clustering has been done on HDLSS asymptotics in recent years.
Liu et al. [6] proposed a two-way split clustering called “statistical signiﬁcance ofclustering(SigClust)” especially for HDLSS data. Ahn et al. [1] proposed a hierarchical divisive clustering and considered its high dimensional asymptotics. Huang
et al. [5] developed the SigClust by Liu et al. [6] with soft thresholding approach.
Yata and Aoshima [9] gave consistency properties of sample principal componentscores and applied it to clustering under high dimensional settings. Nakayama et al.
[7] investigated clustering by kernel principal component analysis for HDLSS data.
Borysov et al. [3] studied behaviors of hierarchical clustering under several asymptotic settings from moderate dimension through HDLSS, nevertheless it is consideredthat theoretical assumptions are strict for HDLSS data due to having discussions on
several asymptotic settings at once. Given this background, we focused on HDLSS
settings and considered asymptotic properties of hierarchical clustering with several
linkage functions.
In this talk, we investigated the hierarchical clustering theoretically in the HDLSS
context as dimension goes to inﬁnity while sample size is ﬁxed. We gave asymptotic
properties of hierarchical clustering and showed the threshold to decide the asymptotic behaviors. Finally, we deliberated performances of the hierarchical clusteringin numerical simulations and actual data analyses.

-1-References[1] Ahn, J., Lee, M.H., Yoon, Y.J. (2012). Clustering high dimension, low sample
size data using the maximal data piling distance. Statistica Sinica, 22, 443–464.
[2] Bhattacharjee, A., Richards, W. G., Staunton, J., Li, C., Monti, S., Vasa, P.,
Ladd, C., Beheshti, J., Bueno, R., Gillette, M., Loda, M., Weber, G., Mark,E. J., Lander, E. S., Wong, W., Johnson, B. E., Golub, T. R., Sugarbaker, D.
J., Meyerson, M. (2001). Classiﬁcation of human lung carcinomas by mRNA
expression proﬁling reveals distinct adenocarcinoma subclasses. Proceedings of
the National Academy of Sciences of the United States of America, 98, 13790–
13795.
[3] Borysov, P., Hannig, J., Marron, J.S. (2014). Asymptotics of hierarchical clustering for growing dimension. Journal of Multivariate Analysis, 124, 465–479.
[4] Eisen, M. B., Spellman, P. T., Brown, P. O., Botstein, D. (1998). Cluster
analysis and display of genome-wide expression patterns. Proceedings of the
National Academy of Sciences of the United States of America, 95, 14863–
14868.
[5] Huang, H., Liu, Y., Yuan, M., Marron, J.S. (2015). Statistical Signiﬁcanceof Clustering using Soft Thresholding. Journal of computational and graphicalstatistics, 24, 975–993.
[6] Liu, Y., Hayes, D.N., Nobel, A., Marron, J.S. (2008). Statistical signiﬁcance ofclustering for high-dimension, low-sample size data. Journal of the American
Statistical Association, 103, 1281–1293.
[7] Nakayama, Y., Yata, K., Aoshima, M. (2021). Clustering by principal component analysis with Gaussian kernel in high-dimension, low-sample-size settings.
Journal of Multivariate Analysis, 185, 104779.
[8] Perou, C. M., Sørlie, T., Eisen, M. B., van de Rijn, M., Jeﬀrey, S. S., Rees,
C. A., Pollack, J. R., Ross, D. T., Johnsen, H., Akslen, L. A., Fluge, O.,
Pergamenschikov, A., Williams, C., Zhu, S. X., Lønning, P. E., Børresen-Dale,
A. L., Brown, P. O., Botstein, D. (2000). Molecular portraits of human breast
tumours. Nature, 406, 747–752.
[9] Yata, K., Aoshima, M. (2020). Geometric consistency of principal componentscores for high-dimensional mixture models and its application. Scandinavian
Journal of Statistics, 47, 899–921.

-2-

Mixed eﬀects modeling of clustered extreme values
鹿児島大学大学院

1

桃木 光輝

鹿児島大学大学院

吉田 拓真

はじめに
気象データから災害リスクを定量化することは，災害対策計画や都市開発の観点から重要である．言い換えると，それは災害を引き起こすような極端な値の発生確率を予測することに相当する．しかしながら，稀な事象を対象とするこの分野では，利用可能なデータが決して多くないという問題がある．極値統計学では，数学的に保証されたパラメトリックな分布を活用することで，このような状況下でもより精度の高い予測を実現する．一方，気象デー
タの特徴は複数の気象観測所でデータが蓄積されている点である．本研究で着目する混合効
果モデルは，全ての観測所のデータを一つのモデリングに統合してデータ全体の共通情報を引き出しつつ，各観測所のデータの分布の差異を効率的・効果的に予測することができる（Sugasawa and Kubokawa 2020)．そのため，混合効果モデルは，利用可能なデータが不足しがちな極値統計学との相性が非常に良い．本研究では，混合効果モデルを応用した新たな極値統計モデルを提案し，予測の不確実性を測るための漸近理論を確立した．
極値統計学において，最も重要視されているのは，極値指数（extreme value index）と呼ばれる，分布の裾の重さに関連するパラメータの予測である．先行研究には，ベイズ階層モ
デルを応用した手法があるが，実データ解析において極値指数の予測結果の大きな不確実性が指摘されていた (Dyrrdal et al. 2015)．本研究では，アプリケーションの範囲をリスクの予測が難しいとされるものに限定することで，極値指数のより詳細な解析が可能な手法を実現した．
また，その他の手法には，多変量極値統計学やコピュラなどがある（Davison et al. 2012）．
これらは全ての観測所のデータの同時分布をパラメトリックに解析する手法であり，主に，観測所間のリスクの依存構造に焦点を当てている．しかしながら，多変量極値統計学の場合，多くの観測所のデータを同時にモデリングすることが困難である（Huster and Wadsworth2020）．また，災害は複雑な要因により引き起こされると考えられるが，これらの手法は位置情報以外の共変量情報をモデリングに組み込むことが困難である．本研究では，災害リスクの地理的な依存構造よりも，災害の要因を明らかにするためのモデリングの開発に専念
した．

-3-

2

モデル
本研究はクラスターデータ（clustered data）

(Yij , Xij ) ∈ R+ × Rp , i = 1, 2, · · · , nj , j = 1, 2, · · · , Jを対象としており，気象データはその一例である．J がクラスターの個数，nj が各クラス
ター内のデータ数であり，(Yij , Xij ) は j 番目のクラスターにおける i 番目の観測である．ただし，Yij は目的変数であり，Xij は説明変数である．

Uj を未知の分散 σ02 を持つ正規分布 N (0, σ02 ) に従う未観測の確率変数とする．このとき，
Xij = x と Uj = uj が与えられた下での Yij の条件付き分布関数 F (y|x, uj ) = P (Yij ≤y|Xij = x, Uj = uj ) に，パレート型分布（Pareto-type distribution）
1 − F (y|x, uj ) = y −1/γ(x,uj ) L(y, x, uj )を仮定する．ここで，γ(x, u) > 0 が極値指数であり，L(y, x, u) は任意の s > 0 に対して

limy→∞ L(ys, x, u)/L(y, x, u) → 1 を満足するような関数である．この分布族には t 分布やパレート分布などの裾の重い分布が多く含まれる．分布の裾における挙動を決定する極値指数の予測が本質的に重要である．本研究では，極値指数に混合効果モデル


log γ(x, uj )−1 = α0 + β0 x + uj ,j = 1, 2, · · · , Jを仮定する．ここで，α0 ∈ R と β0 ∈ Rp は未知の回帰係数である．上記モデルは，Uj ≡ 0の古典的な線形モデルについて研究した Wang and Tsai (2009) の拡張である．α0 と β0 は
全クラスター共通のパラメータである．Uj = uj は変量効果と呼ばれ，これにより説明変数

Xij = x では説明のつかない影響，すなわち，クラスター毎の分布の違いを考慮できる．また，その実現値そのものを予測することも可能である．
講演では，研究を通して得られた提案モデルの数学的性質と数値パフォーマンスについて
報告した．

参考文献
[1] Davison, A.C., Padoan, S.A., and Ribatet, M. (2012) Statistical modeling of spatial extremes. Statistical
Science. 27 161–186.
[2] Dyrrdal, A.V., Lenkoski, A., Thorarinsdottir, T.L., and Stordal, F. (2015) Bayesian hierarchical modeling of extreme hourly precipitation in Norway. Environmetrics, 26 89–106.
[3] Sugasawa, S., and Kubokawa, T. (2020) Small area estimation with mixed models: a review. Japanese
Journal of Statistics and Data Science, 3 693–720.
[4] Wang, H., and Tsai, C. L. (2009) Tail index regression. Journal of the American Statistical Association.
104 1233–1240.

-4-ౡཋ̾˳፭ồࣖဇỂẨỦጮụᡉẲӧᏡễ᩼‫ࠝܭ‬ἁἻἋἑὊໜᢅᆉ
ᓥ㇂೺୍㑻࣭⤫ィᩘ⌮◊✲ᡤ 㔝እ⎔ቃୗࡢ᳜≀㞟ᅋ࡟ࡣࠊ✀Ꮚᩓᕸࡀぶಶయ࠿ࡽ㏆㊥㞳࡟㝈ࡽࢀࡿࡇ࡜ࠊ࠶
ࡿ࠸ࡣ㐺ࡋࡓ⎔ቃࡀ㝈ࡽࢀࡿࡇ࡜࠿ࡽࠊ㞟୰ศᕸࢆ♧ࡍࡶࡢࡀከ࠸ࠋ✀Ꮚᩓᕸไ㝈࡟㉳ᅉࡍࡿ㞟୰ศᕸࢆ⾲⌧ࡍࡿࣔࢹࣝࡢᇶᮏࡢࡦ࡜ࡘ࡛࠶ࡿ Thomas 㐣⛬࡜ࡣࠊ㸯㸬ẕぶࡣࣛࣥࢲ࣒࡟ศᕸࡍࡿ㸦ᐃᖖ࣏࢔ࢯࣥ㐣⛬࡟ᚑ࠺㸧ࠋ㸰㸬ྛẕぶࡣ࣏࢔ࢯࣥศᕸ࡟ᚑ࠺ಶᩘࡢፉࢆ⏕⏘ࡍࡿࠋ㸱㸬㸲㸬ፉࡣẕぶ࠿ࡽṇつศᕸ࡟ᚑࡗ࡚㏆ഐ࡟ᩓᕸࡉࢀࡿࠋ
ẕぶࡣṚࢇ࡛ࠊፉࡓࡕࡀ㞟୰ศᕸࢆ♧ࡍࠋ࡜࠸࠺࢔ࣝࢦࣜࢬ࣒࡛Ⅼศᕸࢆ⏕ᡂࡍࡿ✵㛫Ⅼ㐣⛬࡛࠶ࡿࠋ
Thomas 㐣⛬࡟ࡣࠊ኱ࡁࡃ㸰ࡘࡢၥ㢟ࡀ࠶ࡿࠋ㸯㸬ᐃᖖᛶࢆ௬ᐃࡋ࡚࠸ࡿࡀࠊ⌧ᐇ࡟ࡣ⏕ṧ⋡ࡣ⎔ቃ࡟౫Ꮡࡍࡿࠋࡍ࡞ࢃࡕࠊ
㠀ᐃᖖ࡞Ⅼ㐣⛬ࣔࢹࣝ࡟ࡍࡿᚲせࡀ࠶ࡿࠋ㸰㸬᳜≀ࡣ⦾Ṫࢆ⧞ࡾ㏉ࡍࠋࡔ࠿ࡽࠊぶ㞟ᅋࡶ㞟୰ศᕸࡔࡗࡓ࡜௬ᐃࡍࡿ࡯
࠺ࡀ⮬↛࡛࠶ࡿࠋࡑࡢぶ㞟ᅋࡶ⎔ቃ࡟౫Ꮡࡋࡓ㠀ᐃᖖ࡞ศᕸ࡜⪃࠼ࡿ࡯࠺ࡀ
ࡉࡽ࡟⮬↛࡛࠶ࡿࠋ
せࡍࡿ࡟ࠊ⧞ࡾ㏉ࡋྍ⬟࡞㠀ᐃᖖⅬ㐣⛬࡛࠶ࡿᚲせࡀ࠶ࡿࠋ ᐇ㝿ࡢ࡜ࡇࢁࠊ㠀ᐃᖖ࡛㞟୰ศᕸࢆ⏕ᡂࡍࡿ࢔ࣝࢦࣜࢬ࣒ࡣᐜ᫆࡟సࡿࡇ࡜
ࡀ࡛ࡁࡿࠋࡋ࠿ࡋࠊ1 ḟ࣮࣓ࣔࣥࢺࠊ2 ḟ࣮࣓ࣔࣥࢺ࡞࡝ࢆᩘᏛ࡜ࡋ࡚ᑟฟࡋ࡞
࠸࡜ࠊࢹ࣮ࢱ࠿ࡽࡢࣃ࣓࣮ࣛࢱ᥎ᐃ࡞࡝ࡢ⤫ィⓗ᥎ᐃࢆ⾜࠼࡞࠸ࠋ 㠀ᐃᖖᛶࢆຍ࠼ࡓࢿ࢖࣐࣭ࣥࢫࢥࢵࢺ㐣⛬ࡣ 2000 ᖺ௦࡟ධࡾࠊᩘከࡃᥦၐࡉ
ࢀ࡚࠾ࡾࠊࣃ࣓࣮ࣛࢱ᥎ᐃἲࡶᩘከࡃヨࡉࢀ࡚࠸ࡿࠋࡋ࠿ࡋࠊࡑࡢከࡃࡀぶ㞟ᅋ࡜ፉ㞟ᅋࡢ㸰ୡ௦࠿ࡂࡾࡢࣔࢹ࡛ࣝࠊ࠿ࡘࠊፉ㞟ᅋ࡬ࡢ㠀ᐃᖖᛶࡋ࠿⪃៖ࡋ࡚࠸࡞࠸ࠋࡼࡾᗈ࠸㠀ᐃᖖᛶࢆຍ࠼ࡿ࡜ࠊ2 ḟ࣮࣓ࣔࣥࢺࢆึ➼㛵ᩘ࡛᭩ࡁୗࡏ࡞࠸
ࡓࡵࠊఱࡽ࠿ࡢ௦᭰ἲ࡛ࣃ࣓࣮ࣛࢱ᥎ᐃࢆ⾜࠺ࡀࠊィ⟬㔞ࡀ⭾኱࡟࡞ࡾࡀࡕ࡛࠶
ࡿࠋ ࡇࡇ࡛ࡣࠊ㞟୰ศᕸࡍࡿ㠀ᐃᖖ࡞ぶ㞟ᅋ࠿ࡽ㞟୰ศᕸࡍࡿ㠀ᐃᖖ࡞ፉ㞟ᅋࢆ⏕ᡂࡍࡿⅬ㐣⛬ࣔࢹ࡛ࣝࠊ2 ḟ࣮࣓ࣔࣥࢺࡲ࡛ึ➼㛵ᩘ࡛᭩ࡁ⾲ࡏࡿࡶࡢࢆᥦၐ
ࡋࡓࠋፉ㞟ᅋࡢ 1 ḟ࣮࣓ࣔࣥࢺ ρD(1)(x)ࠊ ḟ࣮࣓ࣔࣥࢺ ρD(2)(x, y)ࡣࠊẕぶ㞟ᅋࡢ ḟ࣮࣓ࣔࣥࢺ ρM(1)(x)ࠊ ḟ࣮࣓ࣔࣥࢺ ρM(2)(x, y)ࢆ⏝࠸࡚ࠊ௨ୗࡢ₞໬ᘧ࡛⾲ࡉ
ࢀࡿࠋ

-5-

³ U ( z)E(u)d ( x; z, s )dz  s( x) U D(1) ( x)

A

(1)
M

2

nl

U D( 2) ( x, y )( ³ U M(1) ( z , w)E(u )d nl ( x; z, s 2 )d nl ( y; z, s 2 )dz  ³ ³ U M( 2) ( z, w)E(u (u  1))d nl ( x; z , s 2 )d nl ( y; w, s 2 )dzdw) AA Au s ( x) s ( y )ࡇࡇ࡛ࠊdnl ( x; P,V 2 ) ࡣᖹᆒ μ, ศᩓ σ2 ࡢṇつศᕸࡢ☜⋡ᐦᗘ㛵ᩘ࡛ࠊẕぶ࠿ࡽፉ࡬ࡢᩓᕸ࡟ᑐᛂࡍࡿࠋ➨㸯㡯ࡣẕぶࡀඹ㏻ࡢጜጒ࡟ᑐᛂࡋࠊ➨㸰㡯ࡣẕぶࡀ␗࡞ࡿ㸰ࡘࡢፉಶయ࡟ᑐᛂࡍࡿࠋu ࡣፉᩘ࡜࠸࠺☜⋡ኚᩘࠊE( )ࡣᮇᚅ್ࢆ⾲ࡍࠋs(x)ࡣ x ࡛ࡢ⏕ṧ⋡࡛ random thinning ࡜ࡋ࡚Ⅼศᕸ࡟స⏝ࡍࡿࠋከࡃࡢሙྜࠊm ಶࡢ⎔ቃ᮲௳ li(x)ࡢ⥺ᙧᅇᖐᘧ࡟ఱࡽ࠿ࡢࣜࣥࢡ㛵ᩘࢆ᪋ࡋࡓᙧ࡛⾲ࡉࢀࡿࠋ⎔ቃࢹ࣮ࢱࡣࢧࣥࣉࣜࣥࢢᆅⅬ࡛ィ ࡋࠊ࢞࢘ࢫ࣮࢝ࢿࣝᖹ⁥໬࡞࡝࡛௵ព
ࡢⅬ x ࡟࠾ࡅࡿ⎔ቃ᮲௳ li(x)࡜ࡋ࡚⏝࠸ࡿࡇ࡜ࡀከ࠸ࠋẕぶ㞟ᅋࡢ 1 ḟ࣮࣓ࣔࣥࢺࠊ2 ḟ࣮࣓ࣔࣥࢺࡀ s(x)ࢆྵࡴࡢ࡛ࠊࡑࢀ࡟ṇつศᕸࢆ஌ࡌࡓ₞໬ᘧࡢ✚ศࡣᙜ↛ࠊึ➼㛵ᩘ࡛ࡣ⾲ࡏ࡞࠸ࠋ࡜ࡇࢁࡀࠊ࢞࢘ࢫ࣮࢝ࢿࣝ࡟ࡼࡿᖹ⁥໬ࢆᅇᖐᘧࡢᚋ࡛⾜࠺࡜࠸࠺㏆ఝࢆ⾜࠺࡜ࠊ⏕ṧ⋡ࡣΰྜṇつศᕸࡢᙧ࡛⾲ࢃࡉࢀࠊୖࡢ₞໬ᘧࡣṇつศᕸࡢ
convolution ࡢබᘧ࡛ explicit ࡟ゎࡅ࡚ࡋࡲ࠺ࠋ⤖ᯝ࡜ࡋ࡚ࠊ⦾Ṫᚋ࡟⏕ṧ⋡࡟ࡼࡿ random thinning ࡜࠸࠺㠀ᐃᖖᛶࡀධࡿ⧞ࡾ㏉ࡋྍ⬟࡞㠀ᐃᖖ㐣⛬ࡣࠊ௨ୗࡢ
ࡼ࠺࡞ᶆ‽ᙧࢆ᭷ࡍࡿࡇ࡜ࡀ♧ࡉࢀࡓࠋU D(1) ( x)

¦ C d ( x, q , Q )
k

0
k

nl

k

2
k

U D( 2) ( x, y )U D(1) ( x) U D(1) ( y ) d nl ( x; y, S12 )¦k Ck1 d nl ( x, vk , Tk2 )d nl ( y, vk , Tk2 ) d nl ( x; y, S 22 )¦h Ch2 d nl ( x, uh , Rh2 )d nl ( y, uh , Rh2 )ࣛࣥࢲ࣒࡞ぶୡ௦࠿ࡽጞࡲࡗࡓ k ୡ௦┠࡞ࡽ࡛͐␎ࡉࢀ࡚࠸ࡿ㡯ࡣ඲㒊࡛ k + 1
ಶ࡟࡞ࡿࠋQ, S, T, R, …ࡣ✀Ꮚᩓᕸࡸ࣮࢝ࢿࣝᖹ⁥໬࡛⏝࠸ࡿṇつศᕸࡢศᩓࡢ
㛵ᩘࠊq, v, u ࡣ⎔ቃࢹ࣮ࢱࡢࢧࣥࣉࣜࣥࢢᆅⅬ㸦ࡢ㛵ᩘ㸧ࠊC ࡣࡑࢀࡽ࡜⏕ṧ⋡ࡢ୰ࡢࣃ࣓࣮ࣛࢱࡢ㛵ᩘ࡛࠶ࡿࠋ

-6-εύʔεੑͷԾఆͷ΋ͱͰͷଟ߲෼෍ͷద߹౓‫ݕ‬ఆ౷‫ྔܭ‬બ୒ͷࢦ਑ͱͦͷԠ༻๺ւಓ‫ڭ‬ҭେɾࡳຈɹछࢢ৴༟

๺ւಓ‫ڭ‬ҭେɾ۴࿏ɹؔ୩༞ཬ

1. φ-μΠόʔδΣϯεʹ‫ͮ͘ج‬ଟ߲෼෍ͷద߹౓‫ݕ‬ఆ౷‫ྔܭ‬ɽ
φ-μΠόʔδΣϯε ͱ‫ݺ‬͹ΕΔɼ෼෍ؒͷඇྨࣅ౓ͷई౓͕ɼCsiszár [3] ͱ Ali andSilvey [1] ʹΑͬͯఏҊ͞Εͨɻ̎ͭͷ཭ࢄ෼෍ɼp = (p1 , . . . , pk ) ͱ q = (q1 , . . . , qk )
ͷؒͷ φ-μΠόʔδΣϯε ଌ౓͸ɼDφ (p, q) = kj=1 qj φ (pj /qj ) ʹΑͬͯఆٛ͞ΕΔɻ͜͜Ͱɼφ(t) ͸ t > 0 ʹରͯ͠ఆٛ͞ΕΔ͍͔ͭ͘ͷੑ࣭Λຬ࣮ͨ͢ತؔ਺Ͱ͋Δɻ
X = (X1 , . . . , Xk ) Λଟ߲෼෍ Multk (n, π) ʹै͏֬཰ม਺ϕΫτϧͱ͢Δɻ͜͜Ͱɼ kj=1 Xj = n, kj=1 πj =1, 0 < πj < 1, (j = 1, . . . , k), Ͱ͋Γɼπ = (π1 , . . . , πk )
͸ະ஌ͷ֬཰ϕΫτϧͰ͋Δɻ͋Δ‫ݻ‬ఆ͞Εͨ֬཰ϕΫτϧ p = (p1 , . . . , pk ) ʹରͯ͠ɼ‫ؼ‬ແԾઆ H0 : π = p Λ‫ݕ‬ఆ͢ΔͨΊʹɼφ-μΠόʔδΣϯε౷‫ ྔܭ‬Kφ ͕
Zografos et al. [6] ʹΑͬͯಋೖ͞ΕͨɻKφ ͸ Kφ = 2nDφ (X/n, p) ʹΑͬͯఆٛ͞
ΕΔɻ͜͜Ͱɼφ(1) = φ (1) = 0 ͱ φ (1) = 1 Ͱ͋Δ (Pardo et al. [4])ɻ
φ ͱͯ͠ɼತؔ਺ φa (t) = {a(a + 1)}−1 {ta+1 − t + a(1 − t)} (a = 0, −1); = t ln t +1 − t (a = 0); = − ln t − 1 + t (a = −1) ΛબͿͱɼφ-μΠόʔδΣϯεଌ౓͸ɼύϫʔμΠόʔδΣϯεଌ౓ʹͳΔ (Read and Cressie [6])ɻΏ͑ʹɼφ-μΠόʔδΣϯε౷‫ ྔܭ‬Kφ ͷ଒͸ɼCressie and Read [2] ʹΑͬͯఏҊ͞ΕͨύϫʔμΠόʔδΣϯε౷‫ ྔܭ‬Ra ͷ଒Λ‫ؚ‬ΜͰ͍ΔɻύϫʔμΠόʔδΣϯε౷‫ ྔܭ‬Ra ͷ଒͕ɼPearson ͷ
X 2 ౷‫( ྔܭ‬a = 1), ର਺໬౓ൺ౷‫( ྔܭ‬a = 0) ౳Λ‫ؚ‬ΜͰ͍ΔɻZografos et al. [7] ͸ɼ‫ؼ‬ແԾઆ H0 ͷ΋ͱͰɼφ-divergence ౷‫ ྔܭ‬Kφ ͸͢΂ͯɼ઴ۙతʹࣗ༝౓ ν = k − 1ͷΧΠೋ৐෼෍ʹै͏͜ͱΛࣔͨ͠ɻ

2. Second-order correction termɽ‫ؼ‬ແԾઆ H0 ͷ΋ͱͰɼ౷‫ྔܭ‬ͷ‫·఺ݪ‬ΘΓͷ s ࣍ϞʔϝϯτΛ E(Kφs |H0 ) = E(Fνs )+mφ (s)+ o(n−1 ),n(s = 1, 2, . . .), ͱͯ͠ධՁ͢Δ͜ͱͷͰ͖Δ౷‫ ྔܭ‬Kφ Λߟ͑Δɻ
͜͜ͰɼFν ͸ɼࣗ༝౓ ν ͷΧΠೋ৐෼෍ʹै͏֬཰ม਺ͱ͢Δɻ͜ͷ࣌ɼmφ (s) ͸ɼsecond-order correction term ͱ‫ݺ‬͹ΕΔɻ΋͠ɼmφ (s) ͷઈର஋͕̌ʹ͚ۙΕ͹ɼͦͷ౷‫ྔܭ‬ͷ෼෍͕ࣗ༝౓ ν ͷΧΠೋ৐෼෍ʹ͍ۙͱߟ͑Δ͜ͱ͕Ͱ͖ΔɻΏ͑ʹɼ mφ (s) ͷઈର஋Λ‫͢ࢉܭ‬Δ͜ͱʹΑͬͯɼͲͷ౷‫͕ྔܭ‬ΧΠೋ৐෼෍ʹ͍͔ۙΛௐ΂Δ͜ͱ͕Ͱ͖ΔɻຊใࠂͰ͸ɼφ-μΠόʔδΣϯε౷‫ྔܭ‬ͷதͰɼ‫ؼ‬ແԾઆͷ΋ͱͰΧΠೋ৐෼෍ʹ࠷΋͍ۙ౷‫ྔܭ‬Λߟ͑Δɻ
3. Second-order correction term ʹؔ͢ΔఆཧKφ ͷ‫ؼ‬ແ෼෍ͷΧΠೋ৐෼෍΁ͷऩଋͷ଎͞ʹؔ͢ΔఆཧΛಋ͘ɻ‫ؼ‬ແԾઆ H0͕‫ݕ‬ఆ౷‫ ྔܭ‬Kφ ͷ mφ (s) ʹؔͯ͠ɼҎԼͷఆཧΛಘΔɻఆཧ 1. φ(t) Λ̐ճඍ෼ՄೳͰɼφ(4) (t) ͕ t = 1 Ͱ࿈ଓͰ͋Δͱ͢Δɻpi =O(k −1 ), (i = 1, . . . , k) ͱԾఆ͢Δͱɼmφ (s) = 0, (s = 1, 2, . . .) Λຬͨ͢ φ ʹରͯ͠ɼH0 ͷ΋ͱͰɼk Λແ‫ݶ‬େʹൃࢄͤ͞Δͱɼ4φ (1) + 3φ(4) (1) ͸̌ʹऩଋ͢Δɻ

-7-‫ؼ‬ແԾઆ H0 ͕ɼp = (1/k, . . . , 1/k) ͷ৔߹ɼsymmetric ͳ‫ؼ‬ແԾઆͱ‫ݺ‬͹ΕΔɻఆཧ 1 ͸ɼPardo [5, p.183] ͕ symmetric ͳ‫ؼ‬ແԾઆͷ΋ͱʹ͓͍ͯɼs = 1, 2, 3 ͷ৔߹ʹࣔͨ݁͠ՌΛɼsymmetric ͳ‫ؼ‬ແԾઆΛ‫ؚ‬ΊͨɼΑΓҰൠతͳ‫ؼ‬ແԾઆͱ s ≥ 4Λ‫ؚ‬Ίͨ͢΂ͯͷ s ʹ֦ுͨ͠΋ͷͰ͋Δɻφ ͱͯ͠ φa Λ༻͍Δ͜ͱʹΑͬͯύϫʔμΠόʔδΣϯε౷‫ʹྔܭ‬ద༻͢Δͱɼ࣍ͷ‫ܥ‬ΛಘΔɻ
‫ ܥ‬1. pi = O(k −1 ), (i = 1, . . . , k) ΛԾఆ͢ΔͱɼύϫʔμΠόʔδΣϯε౷‫ ྔܭ‬Raͷ଒ʹରͯ͠ɼa ʹ͍ͭͯͷ 2 ࣍ํఔࣜ mφa (s) = 0, (s = 1, 2, . . .) ͷ̎ͭͷղ͸ɼH0ͷ΋ͱͰɼk Λແ‫ݶ‬େʹൃࢄͤ͞Δͱɼa = 1 ͱ a = 2/3 ʹऩଋ͢Δɻ‫ ܥ‬2. pi = O(k −1 ), (i = 1, . . . , k) ͱԾఆ͢Δͱɼ౷‫ ྔܭ‬R1 ͱ R2/3 ͷ s ࣍Ϟʔϝϯτͷ mφ (s) ͸ɼmφ1 (s) = As k s + O(k s−1 ), mφ2/3 (s) = Bs k s + O(k s−1 ), (s =
1, 2, . . .) ͱධՁ͞ΕΔɻ͜͜ͰɼAs = (1/2)s(s − 1) (S/k 2 ) − 1 , Bs = (s/27)(s − 1)(S/k 2 ) − 3(2s − 3) (s = 1, 2, . . .) Ͱ͋Δɻͨͩ͠ɼS = ki=1 1/pi ͱ͢Δɻ‫ ܥ‬2 ΑΓɼ1 ≤ S/k 2 < 33/29 ͷͱ͖ɼ|A1 | < |B1 | ͱ |As | < |Bs |, (s ≥ 2) ͕੒Γཱͪɼଞํɼ39/29 ≤ S/k 2 ͷͱ͖ɼ|As | > |Bs |, (s ≥ 2) ͕੒Γཱͭ͜ͱ͕Θ͔Δɻ
͜ΕΒͷٞ࿦ΑΓɼ΋͠ɼΑΓ଎͘ΧΠೋ৐‫ݶۃ‬෼෍ʹऩଋ͢Δ‫ݕ‬ఆ౷‫ྔܭ‬Λબͼ͚ͨΕ͹ɼύϫʔμΠόʔδΣϯεʹ‫ͮ͘ج‬౷‫ ྔܭ‬Ra ͷதͰ͸ɼR2/3 ΍ R1 Λਪ঑͢Δɻ͞Βʹɼଟ߲෼෍ͷద߹౓‫ݕ‬ఆʹ͓͚Δ R1 ͱ R2/3 ͷؒͷൺֱʹؔͯ͠͸ɼ1 ≤ S/k 2 < 33/29 ͷ৔߹ʹ͸ R1 Λਪ঑͢Δɻ
Second-order correction term Λ༻͍ͨ͜ͷΑ͏ͳٞ࿦͸छʑͷ෼ׂදͷಠཱੑ‫ݕ‬ఆ΁ͷԠ༻͕ՄೳͰ͋Δɻ

4. εύʔεੑͷԾఆଟ߲෼෍ͰɼαϯϓϧαΠζ n ͕૿Ճ͢Δͱಉ࣌ʹΧςΰϦʔ਺ k ͕૿Ճ͢Δͱ͍͏৚݅΍ɼΧςΰϦʔ਺ k ͷ૿Ճʹ൐ͬͯαϯϓϧαΠζ n ΋૿Ճ͢Δ৚݅ɼ෼ׂ
දͰɼαϯϓϧαΠζ͕૿Ճ͢Δͱ͖ηϧ਺͕૿Ճ͢Δͱ͍͏৚݅͸ɼଟ߲෼෍ͷεύʔεੑͷԾఆ΍෼ׂදͷεύʔεੑͷԾఆͱ‫ݺ‬͹ΕΔɻεύʔεੑͷԾఆ͸ɼ࣮ࡍͷ໰୊Ͱ͸ɼΧςΰϦʔ਺ʹରͯ͠αϯϓϧαΠζ͕গͳ͍৔߹΍ɼηϧ਺ʹରͯ͠αϯϓϧαΠζ͕গͳ͍৔߹ͱͯ͠ղऍ͞ΕΔɻ
ຊใࠂ͸͜ͷΑ͏ͳಛ௃ͷ͋Δσʔλʹର͢Δ౷‫ྔܭ‬બ୒ͷࢦ਑ͱͳΔɻࢀߟจ‫ݙ‬[1] Ali, S. M. and Silvey, D., J. Roy. Statist. Soc. B 28 (1966) 131–142. [2] Cressie,
N. and Read, T. R. C., J. Roy. Statist. Soc. B 46 (1984) 440–464. [3] Csiszár,
I., Studia Sci. Math. Hungar. 2 (1967) 299–318. [4] Pardo, L., Pardo, M. C. andZografos, K., J. Japan Statist. Soc. 29 (1999) 213–228. [5] Pardo, L., Statisticalinference based on divergence measures, Chapman & Hall/CRC, 2006. [6] Read, T.
R. C. and Cressie, N. A. C., Goodness-of-ﬁt statistics for discrete multivariate data,
Springer-Verlag, New York, 1988. [7] Zografos, K., Ferentions, K. and Papaioannou,
T., Commun. Statist.-Theory Meth. 19 (1990) 1785–1802.

-8-

Simultaneous estimation of multiplicative Poisson means in two-way contingency
tablesYuan-Tsung Chang (Mejiro University), Shinozaki Nobuo (Keio University)

Abstract
Shrinkage estimation of Poisson means is considered when observations are given in the form of a two-way contingency table. Assuming a multiplicative Poisson model, estimators which shrink to the speciﬁed values or an
order statistic in one dimension and in two dimensions are considered and are shown to dominate the maximum
likelihood estimator (MLE) under normalized squared error loss.

1

Introduction

We consider two-way multiplicative model where xij , i = 1, . . . , I, j = 1, . . . , J, are independent random Poisson
random variables with means
λij = λαi βj , i = 1, . . . , I, j = 1, . . . , J,
I
J
where αi ≥ 0 and βj ≥ 0 satisfy i=1 αi = 1 and j=1 βj = 1, respectively. We denote the one-dimensionalfrequencies and the total frequency by
xi+ =

J


xij , i = 1, . . . , I, x+j =

j=1

I


xij , j = 1, . . . , J, x++ =

i=1

I 
J


xij .

i=1 j=1

As discussed in Hara and Takemura (2006) complete suﬃcient statistics are x1 = (x1+ , . . . , xI+ ) and x2 =(x+1 , . . . , x+J ). The MLE of λij is xi+ x+j
L
x++
λ̂M
=
ij
0

if

x++ = 0

if

x++ = 0.

They have given a class of improved estimators which shrink the MLE toward the origin under the normalized
squared error loss. The simple one isdxi+ x+jHT, i = 1, . . . , I, j = 1, . . . , J.
1−
δij =
x++
x++ + dThe following lemma is a special case of Lemma 2.1 of Hara and Takemura (2006) and is useful to evaluate therisk of the shrinkage estimators when normalized squared error loss is concerned.
Lemma 1.1. If g(x1 , x2 ) is a real-valued function satisfying E|g(x1 , x2 )| < ∞ and g(x1 , x2 ) = 0 when xi+ = 0
or x+j = 0, theng(x1 , x2 )(x++ + 1)
Eg(x1 + eIi , x2 + eJj ) ,=E
λij(xi+ + 1)(x+j + 1)
where eIi (eJj ) is I × 1 (J × 1) unit vector with i-th (j-th) component 1.
Next section we consider one-dimensional shrinkage to an order statistic or a speciﬁed point.

2

One-dimensional shrinkage to an order statistic or a speciﬁed pointFirst, we consider one-dimensional shrinkage to an order statistic.
Let x()+ be the -th smallest observation among x1+ , . . . , xI+ . We assume that I ≥  + 2 and consider thefollowing estimator which shrinks xi+ toward x()+ when xi+ ≥ x()+ :(xi+ − x()+ )+x+j(1)
xi+ − ϕ(W ), i = 1, . . . , I, j = 1, . . . , J,
δij =
x++
W +d
where W =I

i=1 (xi+ − x()+ )+, a+ = max(0, a) and d is a positive constant. Then we have the following.

-9-

Theorem 2.1. Suppose that ϕ(W ) is a non-decreasing function satisfying 0 ≤ ϕ(W ) ≤ 2(I −  − 1) and that d ≥I(1)Lsup ϕ(W )/2. Then δij , i = 1, . . . , I improves upon the MLE λM
ij , i = 1, . . . , I under the loss function
i=1 (λ̂ij −λij )2 /λij for any j = 1, . . . , J.
L
to a speciﬁed non-negative values, bi .
Next we consider the estimators shrink λ̂Mij
Let bi ≥ 0, i = 1, . . . , I be given numbers and we propose the following shrinkage estimator which shrinks xi+
to bi when xi+ ≥ bi :x+j(xi+ − bi )+
(2)
δij =
, i = 1, . . . , I, j = 1, . . . , J,
xi+ − ϕ(N, W )x++
W + d(N )where W =I

i=1 (xi+ − bi )+and N = #{i|xi+ ≥ bi }. Then we have the following.

Theorem 2.2. Suppose that ϕ(N, W ) is a non-decreasing function of W and satisﬁes 0 ≤ ϕ(N, W ) ≤ 2(N − 1)+(2)
for any 0 ≤ N ≤ I. Suppose that d(N ) ≥ supW ϕ(N, W )/2. Then δij , i = 1, . . . , I improves upon the MLE

I
L
2
λ̂M
ij , i = 1, . . . , I under the loss function
i=1 (λ̂ij − λij ) /λij for any j = 1, . . . , J.
It may be noticed that the shrinkage is made only when N ≥ 2.
Remark Theorems 2.1 and 2.2 can be generalized directly to the case of Poisson multiplicative model for a
multi-way contingency tables.
Next section we also consider two-dimensional shrinkage to order statistics or to a speciﬁed point.

3

Two-dimensional shrinkage to order statistics.

Let x()+ and x+(m) be the -th and m-th smallest observation among x1+ , . . . , xI+ and x+1 , . . . , x+J , respectively.
We assume that I ≥  + 2 and J ≥ m + 2 and consider the estimator which shrinks xi+ toward x()+ whenxi+ ≥ x()+ in the ﬁrst dimension and shrinks x+j toward x+(m) when x+j ≥ x+(m) in the second dimensionLsimultaneously. To improve upon the MLE λ̂Mij , we propose the following estimator :
(3)
δij =

1
x++



(xi+ − x()+ )+xi+ − ϕ1 (W1 )W1 + d1(x+j − x+(m) )+x+j − ϕ2 (W2 ),W2 + d2

i = 1, . . . , I, j = 1, . . . , J,

I
J
where W1 = i=1 (xi+ − x()+ )+ and W2 = j=1 (x+j − x+(m) )+ and d1 and d2 are positive constants. Then we
have the following.
Theorem 3.1. Suppose that ϕ1 (W1 ) and ϕ2 (W2 ) are non-decreasing functions satisfying 0 ≤ ϕ1 (W1 ) ≤ I −−1
and 0 ≤ ϕ2 (W2 ) ≤ J − m − 1, respectively. If d1 ≥ (I −  − 1)/(I − ) sup ϕ1 (W1 ) and d2 ≥ (J − m −(3)
L1)/(J − m) sup ϕ2 (W2 ). Then δij , i = 1, . . . , I, j = 1, . . . , J improves upon the MLE λ̂Munder the loss function
ij
I J
2
(
λ̂
−
λ
)
/λ
.
ij
ij
i=1
j=1 ijNext, we consider two-dimensional shrinkage to a speciﬁed point.
Let bi ≥ 0, i = 1, . . . , I and cj ≥ 0, j = 1, . . . , J be given numbers. Assuming that I, J ≥ 2, we shrink xi+ to bi
Lwhen xi+ ≥ bi and x+j to cj when x+j ≥ cj . To improve upon the MLE λ̂Mij , we propose the following estimator
(4)
δij =

1



x++

(xi+ − bi )+xi+ − ϕ1 (N1 , W1 )W1 + d1 (N1 )
I

+

J


(x+j − cj )+x+j − ϕ2 (N2 , W2 ),W2 + d2 (N2 ), W2 =

j=1 (x+j − cj )i = 1, . . . , I, j = 1, . . . , J,+

, N1 = #{i|xi+ ≥ bi , i = 1, . . . , I} and N2 = #{j|x+j ≥IJci , j = 1, . . . , J}. Although it may be natural to put the condition i=1 bi = j=1 cj , we do not need it in thefollowing.
where W1 =i=1 (xi+ − bi )Theorem 3.2. Suppose that ϕi (Ni , Wi ) is a non-decreasing function of Wi and satisﬁes 0 ≤ ϕi (Ni , Wi ) ≤ (Ni −1)+
(4)
for any Ni ≥ 0, and that di (Ni ) ≥ (Ni − 1)+ /Ni supWi ϕi (Ni , Wi ), for any Ni ≥ 0, i = 1, 2. Then δij improves


I
J
L
under the loss function i=1 j=1 (λ̂ij − λij )2 /λij .
upon the MLE λ̂Mij
It may be noticed that the shrinkage in the i-th dimension is made only when Ni ≥ 2.
For the rest of the content, please refer to ”Simultaneous estimation of Poisson means in two-way contingency
tables under normalized squared error loss“, JJSDS (2022) vol.5, issue 2 p577-628.

- 10 -

正方分割表における点対称性からの隔たりを測る
尺度の推定について
東京理科大学大学院 理工学研究科 中島 亮
東京理科大学大学院 理工学研究科 星野 瞭太
東京理科大学 理工学部 藤澤 健吾
東京理科大学 理工学部 田畑 耕治

1. APS モデルからの隔たりを測る尺度の推定について
行変数 X と列変数 Y が順序のある同じ分類からなる R × R 正方分割表において, (i, j)
セル確率を pij (i = 1, . . . , R; j = 1, . . . , R) とする. 点対称 (PS) モデルは次のように定義される (Wall and Lienert, 1976):

pij = pi∗ j ∗(1 ≤ i, j ≤ R),

ただし, i∗ = R + 1 − i, j ∗ = R + 1 − j である.
アナザー点対称 (APS) モデルは次のように定義される (Kurakami et al., 2017):

pij = pi∗ j ∗(i + j = R + 1).

モデルが与えられたデータに適合しない場合, モデルからの隔たりを測る尺度に関心がある. Iki and Tomizawa (2019) は pij + pi∗ j ∗ > 0 (i = 1, . . . , R; j = 1, . . . , R) を仮定
し, APS モデルからの隔たりを測る尺度を次のように提案した:
1 
2pij
ΦAP S =
pij log,
Δ log 2
pij + pi∗ j ∗

i+j=R+1

ただし, Δ =i+j=R+1 pij である.
 
観測度数 nij がサンプル数 n (n = i j nij ) の多項分布に従うと仮定し, p を R2 × 1
多項確率ベクトルとする. すなわち,

p = (p11 , p12 , . . . , p1R , p21 , p22 , . . . , p2R , . . . , pR1 , pR2 , . . . , pRR )t ,ここで t は転置を表す. p̂ij を標本比率 (p̂ij = nij /n) とし, pij を p̂ij で置き換えたベクト
ルを p̂ とする. ΦAP S の推定量 Φ̂AP S の漸近バイアスを次のように与えた:
 2
∂ ΦAP S
1t
E(Φ̂AP S − ΦAP S ) =tr(D(p) − pp ) ,2n∂p∂pt
ここで D(p) は p の i 番目の要素を i 番目の対角要素とする対角行列, tr は行列のトレースを表す. このとき, 次の推定量を提案した:
1∂ 2 Φ̂AP S
Φ̃AP S = Φ̂AP S − tr(D(p̂) − p̂p̂t ) ,2n∂ p̂∂ p̂t
ここで, [∂ 2 Φ̂AP S /∂ p̂∂ p̂t ] は [∂ 2 ΦAP S /∂p∂pt ] の pij を p̂ij で置き換えたものである.

- 11 -

2. RGS, CPS モデルからの隔たりを測る尺度の推定について
逆グローバル対称 (RGS) モデルは次のように定義される (Kurakami et al., 2017):



ΔU = Δ L ,ただし, ΔU =i+j<R+1 pij , ΔL =i+j>R+1 pij である.
条件付き点対称 (CPS) モデルは次のように定義される (Tomizawa, 1986):

pij = τ pi∗ j ∗(i + j < R + 1).

Iki and Tomizawa (2019) は RGS モデルと CPS モデルからの隔たりを測る尺度をそれ
ぞれ次のように提案した:


1
ΔcU
ΔcL
c
c
+ ΔL log
ΔU log,
ΦRGS =
log 2
1/2
1/2


Δp∗ij
Δp∗i∗ j ∗1 ΦCP S =+ pi∗ j ∗ log
pij log,
Δ log 2 i+j<R+1
ΔU
ΔL
ただし, ΔcU = ΔU /Δ, ΔcL = ΔL /Δ, p∗ij = pij /(pij + pi∗ j ∗ ) である.
APS モデルの場合と同様にして, 尺度の推定量をそれぞれ次のように提案した:∂ 2 Φ̂RGS
1
Φ̃RGS = Φ̂RGS − tr(D(p̂) − p̂p̂t ) ,2n∂ p̂∂ p̂t
1∂ 2 Φ̂CP S(D(p̂) − p̂p̂t ) .
Φ̃CP S = Φ̂CP S − trt2n∂ p̂∂ p̂Φ̂AP S , Φ̂RGS , Φ̂CP S の漸近バイアスを n 倍したものを ΛAP S , ΛRGS , ΛCP S とする. このとき, 次の定理及び系を得た:
定理 1. ΛAP S は ΛRGS と ΛCP S の和に等しい.
系 1. Φ̃AP S は Φ̃RGS と Φ̃CP S の和に等しい.
定理 2. Φ̂AP S , Φ̂RGS , Φ̂CP S の漸近バイアスの比率は R(R − 1)/2 : 1 : (R + 1)(R − 2)/2
である.
参考文献

[1] Iki, K. and Tomizawa, S. (2019). Measure of departure from point symmetry anddecomposition of measure for square contingency tables. Journal of Statistical Theory
and Applications, 19, 526-533.
[2] Kurakami, H., Negishi, K., and Tomizawa, S. (2017). On decomposition of pointsymmetry for square contingency tables with ordered categories. Journal of Statistics:
Advances in Theory and Applications, 17, 33-42.
[3] Tomizawa, S. (1986). Four kinds of symmetry models and their decompositions in asquare contingency table with ordered categories. Biometrical Journal, 28, 387-393.
[4] Wall, K. D. and Lienert, G. A. (1976). A test for point-symmetry in J-dimensional
contingency-cubes. Biometrical Journal, 18, 259-264.

- 12 -

正方分割表におけるモーメントに基づく
周辺同等性からの隔たりを測る尺度
足立 匠 1
1

安藤 宗司 2

田畑 耕治 2

東京理科大学大学院 理工学研究科 情報科学専攻
2

東京理科大学 理工学部 情報科学科

順序カテゴリ r × r 正方分割表において，行変数と列変数をそれぞれ X と Y とし，(i, j) セル確率を

pij = Pr(X = i, Y = j) とする (i = 1, . . . , r; j = 1, . . . , r)．正方分割表解析では，観測度数が主対角セルに集中する傾向があるため，行変数と列変数の独立性は成り立たないことが多い．そのため，独立性に代わり，
行変数と列変数の対称性に関して関心がある．行変数と列変数の対称性を表す代表的なモデルとして，周辺同
等 (MH) モデル (Stuart [1]) がある．

MH モデルは次のように表される．
ただし，pi· =r

t=1 pit ,p·i =pi· = p·i(i = 1, . . . , r),r

s=1 psi である．MH モデルは様々な表現があることが知られている．行変

数 X と列変数 Y の周辺累積確率を用いると，MH モデルは次のようにも表される．

FiX = FiY
ただし，

FiX =

i


(i = 1, . . . , r − 1),

ps· = Pr(X ≤ i),

FiY =

s=1

i


p·t = Pr(Y ≤ i).

t=1

周辺累積確率 FiX と FiY の差を考えると，MH モデルは次のようにも表される．

G1(i) = G2(i)(i = 1, . . . , r − 1),

(1)

ただし，

G1(i) =ir

pst = Pr(X ≤ i, Y ≥ i + 1),

G2(i) =

s=1 t=i+1ri

pst = Pr(X ≥ i + 1, Y ≤ i).

s=i+1 t=1

与えられたデータに対して MH モデルの当てはまりが悪い場合，(i) MH モデルよりも制約の弱いモデルを適用すること，(ii) MH モデルの当てはまりが悪い原因を分析すること，(iii) MH モデルからの隔たりの程度を測ることに関心がある．MH モデルの表現として式 (1) に注目して，これら (i) から (iii) に関する先行研究を述べる．
式 (1) に基づく周辺非同等モデルとして，Tahata and Tomizawa [2] は，m 次パラメータ周辺同等 (MH(m))
モデルを提案した．既知の m (m = 1, . . . , r − 1) に対して，MH(m) モデルは次のように定義される．
m−1

G1(i) =k

ψki G2(i)(i = 1, . . . , r − 1).

k=0

ψ0 = ψ1 = · · · = ψm−1 = 1 のとき，MH(m) は MH モデルに一致する．MH(1) モデルは拡張周辺同等(EMH) モデル (Tomizawa [3])，MH(2) モデルは一般化周辺同等モデル (Tomizawa [4]) に一致する．
式 (1) に基づく MH モデルよりも制約の弱いモデルとして，Tahata and Tomizawa [2] は，k 次モーメント
周辺一致 (k-MME) モデルを提案した．既知の正の整数 k に対して，k-MME モデルは次のように定義される．

E[X k ] = E[Y k ],- 13 -ただし，

E[X k ] =rsk ps· ,E[Y k ] =

r


s=1

tk p·t .

t=1

さらに，Tahata and Tomizawa [2] は，k-MME モデルの別表現を次のように与えた．
r−1


(i + 1) − ikkG1(i) =

i=1

r−1



(i + 1)k − ik G2(i) .

i=1

MH モデルの当てはまりが悪い原因を分析するために，Tahata and Tomizawa [2] は，「MH モデルが成り立つことと，全ての k = 1, . . . , r − 1 に対して k-MME モデルが成り立つことは必要十分である」という分解定理を与えた．この分解定理から，どのモーメントが一致していないことにより MH モデルの当てはまりが悪くなったのかを特定することが可能になる．

Tomizawa, Miyamoto and Ashihara [5] は，G1(i) + G2(i) > 0 (i = 1, . . . , r − 1) を仮定し，式 (1) に基づ
く MH モデルからの隔たりの程度を測る尺度を次のように提案した．

r−1 G∗1(i)
G∗2(i)1  ∗∗Ψ=+ G2(i) log
G1(i) log,
log 2 i=1
Q∗i
Q∗i
ただし，

G∗1(i) =

G1(i)
,
Δ

G∗2(i) =

G2(i)
,
Δ

Q∗i =

G∗1(i) + G∗2(i)
2

,

Δ=

r−1


(G1(i) + G2(i) ).

i=1

尺度 Ψ は次の性質がある．

• 0 ≤ Ψ ≤ 1，
• Ψ = 0 ⇔ G1(i) = G2(i) (i = 1, . . . , r − 1) ⇔ MH モデル，
• Ψ = 1 ⇔「G1(i) = 0 かつ G2(i) > 0」または「G1(i) > 0 かつ G2(i) = 0」(i = 1, . . . , r − 1)．
分解定理は，MH モデルの当てはまりが悪い原因を特定するには有用であるが，MH モデルからの隔たりの程度を測ることは難しい．一方，尺度 Ψ は MH モデルからの隔たりの程度を測ることには有用であるが，MHモデルの当てはまりが悪い原因を特定することは難しい．本講演では，MH モデルの当てはまりが悪い原因を特定することが可能な MH モデルからの隔たりの程度を測る尺度を提案した．提案尺度により，与えられたデータに対して MH モデルの当てはまりが悪い場合，そのデータの特徴をより詳細に分析することが可能になると期待される．

参考文献[1] A. Stuart. A test for homogeneity of the marginal distributions in a two-way classiﬁcation. Biometrika,
42(3/4):412–416, 1955.
[2] K. Tahata and S. Tomizawa. Generalized marginal homogeneity model and its relation to marginal
equimoments for square contingency tables with ordered categories. Advances in Data Analysis andClassiﬁcation, 2(3):295–311, 2008.
[3] S. Tomizawa. Diagonals-parameter symmetry model for cumulative probabilities in square contingency
tables with ordered categories. Biometrics, 49(3):883–887, 1993.
[4] S. Tomizawa. A generalization of the marginal homogeneity model for square contingency tables withordered categories. Journal of Educational and Behavioral Statistics, 20(4):349–360, 1995.
[5] S. Tomizawa, N. Miyamoto, and N. Ashihara. Measure of departure from marginal homogeneity forsquare contingency tables having ordered categories. Behaviormetrika, 30(2):173–193, 2003.

- 14 -

国民生活基礎調査による日本の所得分布と所得格差
西埜晴久 ∗

17 Dec. 2022†

はじめに

1

国民生活基礎調査は，厚生労働省が行っている統計である．昭和 61(1986) 年から始まり，3
年ごとに大規模調査が行われている．近年では，2019 年，2016 年，2013 年，2010 年が大規模調査が行われた年である．大規模調査の年の間の各年に簡易調査が行われている．なお，
2020 年調査の簡易調査は調査が中止されている．
大規模調査の年では，世帯票，健康票，所得票などがある．そのうちの所得票の個票デー
タを用いて所得分布の分析を行う．直近の大規模調査は 2019 年であり，その際の集計客体
数は世帯票・健康票で 217,179 世帯，所得票で 22,288 世帯である．
以下に公表されている所得票の数値を掲げる．

Table 1: 所得票の数値
年
調査期間
集計客体数
1998 年（H10） H9 年 1 月〜12 月
30506
2010 年（H22） H21 年 1 月〜12 月
26115
H30 年 1 月〜12 月
22288
2019（R1）

平均所得
657.7
549.6
552.3

中央値
536
438
437

なお，平均所得，中央値は総世帯で単位は (万円) である．

所得分布の推定2

1998 年の個票データの総所得から平均値，中央値，ジニ係数を直接計算すると平均値 657.96,
中央値 536, ジニ係数 0.4042 となる．これらの数値は公表されている値とほぼ同じである．
また，個票データを用いて，一般化ベータ分布，Singh-Maddala および第 2 種の一般化ベー
タ分布の所得分布に用いられる確率分布を最尤法で推定し，パラメータの推定値を得ることが出来た．
拡大乗数について：2000 年代以降の調査では，拡大乗数を用いて調整することで，平均所得などを計算している．拡大乗数とは，ある地域の国勢調査で得られる世帯数と所得票の調査世帯数との比である．実際には県別でかつ県内でも指定都市とそれ以外に分けて拡大乗数を算出している．一般的には，都市部の世帯の標本サイズが郡部の世帯の標本サイズに比べ小さいので，拡大乗数が大きくなる傾向がある．つまり，拡大乗数とは，各地域に対してウェイトを与えていることになるので，そのウェイト（拡大乗数）を考慮して対数尤度を求めた方が良い．
Table 2 には個票データ（個別）と 10 分位のデータで推定したジニ係数を掲げた．
∗†

広島大学大学院人間社会科学研究科
科研費シンポジウム「統計科学の開拓」金沢大学サテライト・プラザ- 15 -

Table 2: ジニ係数

直接

GG
SM
GB2

H10
0.40420個別
10 分位
0.4096 0.4033
0.4011 0.3835
0.3991 0.3824

H22
0.40326
個別
10 分位
0.4060 0.3912
0.4005 0.3970
0.4016 0.3908

R01
0.40910個別
10 分位
0.4120 0.4037
0.4065 0.3987
0.4076 0.3954

Table 3: ξ の推定結果
H10
H22
R01
ξ の推定値
上位 1%
0.2124 0.4056 0.1899
上位 0.5%0.0418 0.1312 0.1952
集計客体数
30506
26115
22288
21
20
5000 万円以上 72

所得分布の上裾の推定

3

Moriguchi and Saez (2008)1 は，税務統計より所得の上位 1%が占める割合と，上位 0.1%の占める割合を調べた．第二次世界大戦前の不平等度は高かったが，戦後の不平等度は小さ
く，米国が 1980 年代以降上位所得者の割合が上昇したのに比べ日本はそれほど上昇していないことを示した．そこでこの場合には分布の上裾がどうなっているかに関心がある．よって分布の上裾の部分を．一般化パレート分布 GP D(ξ, σ, u) を当てはめて形状パラメ―タ ξの推定を行った．形状パラメ―タ ξ の推定結果は Table 3 に掲げてある．

今後の課題

4

GB2 の推定は不安定なので，初期値の選び方などもう少し安定して推定する方法を考えた
い．また，今回取り上げた以外の分布を推定する．
所得票のデータと世帯票のデータを接続することで各世帯の世帯数を把握することが出来る．世帯数が分かると所得を世帯人数の平方根で割った等価所得を求めることが出来る．
年につれて世帯人員数（世帯の構成員数）が減少する傾向があるので，等価所得を用いると世帯人数の減少の影響を除いて所得の変化を調べることができる．
ただし，等価可処分所得については，所得票に公表されている年もある．
「等価可処分所得」とは，下記によって算出される．
等価可処分所得 ＝ （当初所得 ＋ 社会保障給付 − 税金 − 社会保険料 − 掛金など） ÷ √世帯人員数

1

Moriguchi, C. and Saez, E. (2008) “The Evolution of Income Concentration in Japan, 1886-2005:
Evidence from Income Tax Statistics” The Review of Economics and Statistics, 90, 713–734.

- 16 -

ૢ࡞ม਺๏ͷཧղ΁ɿ‫ྔܭ‬ੜ෺ͱ‫ࡁܦྔܭ‬ͷᬍ᫯
(Toward understanding the Instrumental Variables Method in Biometrics
and Econometrics)
2022 ೥ 12 ݄ 17 ೔
ࠃ༑௚ਓ1‫ݴݤ‬༿ (Key Words): ౷‫ܭ‬తҼՌਪ࿦, ྟচࢼ‫ݧ‬, ੓ࡦධՁ, Noncompliance, ATE, LATE, ૢ࡞ม਺๏, ߏ଄ํఔࣜ, ࣝผੑ, TSLS, LIML,
Mendelian Randomization, Ҩ఻ࢠӸֶཁ໿ : ҼՌؔ܎ (causality) ͸౷‫ܭ‬ՊֶΛ‫ؚ‬ΊॾՊֶʹͱͬͯ͸‫ج‬ຊత͔ͭॏཁͳ෼ੳର৅Ͱ͋Δɻ‫ྔܭ‬ੜ෺ͱ‫ࡁܦྔܭ‬ͷ෼໺Ͱ͸͜ͷؒɺ౷‫ܭ‬తҼՌਪ࿦ (statistical causal inference) ͕੝ΜʹԠ༻͞Ε͍ͯΔɻຊใࠂͰ͸·ͣ Rubin (1974) ʹ࢝·Δ൓࣮Ծ૝ (counter-factual) Ϟσϧͱ
Angrist, Imbens and Rubin (1996, ུͯ͠ AIR) ʹΑΔૢ࡞ม਺๏ (instrumental variables method) ͷԠ༻ͷҙຯΛઆ໌ͨ͠ɻ࣍ʹ఻౷తͳ‫ܦྔܭ‬ࡁֶ (econometrics) ʹ͓͚Δಉ࣌ํఔࣜͱߏ଄ํఔࣜ (structural equation) Λ؆୯ͳྫΛ༻͍ͯઆ໌͠ɺߏ଄ํఔࣜΛ༻͍ͨ౷‫ܭ‬తҼՌؔ܎ͷղऍΛड़΂ɺͦͷ౷‫ܭ‬తਪఆ๏Λٞ࿦ͨ͠ɻߏ଄ํఔࣜͷਪఆͰ͸ OLS
๏ (࠷খೋ৐๏) ͸ҰகੑΛ࣋ͨͳ͍ͷͰɺૢ࡞ม਺๏ (IV ๏) ͱͯ͠ͷ
Wald ๏ɺLIML(੍‫ݶ‬৘ใ࠷໬๏, ෼ࢄൺ࠷খ๏)ɺTSLS(2 ஈ֊࠷খೋ৐๏)ɺGMM(ҰൠԽੵ཰๏) ͳͲͷ௕ॴͱ୹ॴΛઆ໌͠ɺ͞Βʹ‫ྔܭ‬ੜ෺ͱ‫͚͓ʹͲͳࡁܦྔܭ‬Δ౷‫ܭ‬తҼՌ෼ੳͷ͞ΒͳΔ՝୊Λల๬ͨ͠ɻ
·ͨ࠷‫ʹޙ‬Ҩ఻ࢠӸֶʹ͓͚Δ MR(Mendelian Randomization, ϝϯσϧɾϥϯμϜԽղੳ) ʹ͓͚Δૢ࡞ม਺๏ͷར༻ʹؔ͢Δ࠷ۙͷల։ʹ͍ͭͯ΋ٞ࿦ͨ͠ɻ
ͳ͓ɺ࣌ؒͷؔ܎Ͱे෼ͳઆ໌͕Ͱ͖ͣʹলུͨ͠ͱ͜Ζ΋গͳ͘ͳ͔ͬ
ͨɻ‫ڀݚ‬ใࠂͷৄࡉͳ಺༰Λ஌Γ͍ͨ৔߹ʹ͸ະఆߘͰ͸͋Δ͕ɺʮૢ࡞ม਺๏ͷཧղ΁ɿ‫ྔܭ‬ੜ෺ͱ‫ࡁܦྔܭ‬ͷᬍ᫯ʯ
http://www.kunitomo-lab.sakura.ne.jp/2022-11-4DP.pdf
Λࢀর͞Ε͍ͨɻ

1ಛ೚‫ڭ‬त, ౷‫਺ܭ‬ཧ‫ ˟ ॴڀݚ‬190-8562 ౦‫ཱ౎ژ‬઒ࢢ྘ொ 10-3

- 17 -

多変量分布間の回帰モデル
東京大学経済学研究科 岡野遼
東京大学総合文化研究科 今泉允聡

1. はじめに説明変数と結果変数が共に分布の形で与えられるような回帰を考える. このような回帰は分布間回帰（distribution

regression）と呼ばれ, 複雑データ解析の一分野として近年関心がもたれている. 例えば, ある年のある国における年
齢別死亡率は, 横軸を年齢, 縦軸を死亡率とする密度関数として表示することで, 一次元分布に値をとるデータとみなすことができる. [1], [2] では世界 37 各国の 2013 年の年齢別死亡率分布を 1983 年のそれに回帰することで, これらの国々の年齢別死亡率分布が時間を経てどのように変化したかを分析している.

W を Rd 上の分布で有限な 2 次モーメントを持つもの全体とし, W に 2-Wasserstein 距離 dW を与えてできる距
離空間 (W, dW ) を Wassersetin 空間と呼ぶ. F を W × W 上の分布とし, (ν1 , ν2 ) ∼ F とする. 分布間回帰は ν1 を説明変数, ν2 を結果変数とするような回帰問題として定式化される.
分布間回帰のモデルは近年いくつか提案されているが, それらのほとんどは分布の次元を d = 1 の場合に限定している. 例えば, [1] は Wasserstein 空間の tangent bundle, [2] は最適輸送写像を用いた分布間回帰モデルをそれぞれ提案しているが, これらのモデルは d = 1 のときに最適輸送問題が closed form な解を持つことを利用しており, d ≥ 2の一般の分布に対してそのまま拡張することは困難である. 本研究では, 多変量分布であっても分布のクラスを適切に制限すれば, 最適輸送問題が closed form な解な解を持つことに注目し, 新たな分布間の回帰モデルを提案する.

2. 背景：ガウス分布間の最適輸送問題
以下では簡単のため, 分布のクラスをガウス分布族を制限した場合を考えるが, より一般にある種の楕円分布に制限しても同様に議論ができる. G を Rd 上の Gauss 分布全体とする. 非退化な共分散行列を持つ二つのガウス分2
布 μ1 = N (m1 , Σ1 ), μ2 = N (m2 , Σ2 ) ∈ G が与えられた時, それらの間の最適輸送写像 tμμ1 及び Wasserstein 距離dW (μ1 , μ2 ) はそれぞれ以下のように陽に表されることが知られている：
−1/2

1/2

1/2

−1/2

tμμ21 (x) = m2 + Σ1 [Σ1 Σ2 Σ1 ]1/2 Σ1 (x − m1 ),

1/2
1/2
dW (μ1 , μ2 ) = m1 − m2 2 + tr[Σ1 + Σ2 − 2(Σ1 Σ2 Σ1 )1/2 ].
また, Σ∗ を非退化な共分散行列とすると, μ∗ = N (m∗ , Σ∗ ) ∈ G における空間 (G, dW ) の接空間は, 内積空間

Tμ∗ = (Rd × Sym(d), Gμ∗ ) によって与えられることが知られている ([3]). ここで, Sym(d) はサイズ d × d の対称行列全体で, z = (a, V ), w = (b, W ) ∈ Rd × Sym(d) に対して, それらの間の内積を Gμ∗ (z, w) = a b + tr(V Σ∗ W )によって定めている. さらに, μ∗ における exponential map expμ∗ : Tμ∗ → G は expμ∗ (a, V ) = N (a + m∗ , (V +I)Σ∗ (V + I)) によって与えられ, μ∗ における logarithmic map logμ∗ : G → Tμ∗ は, logμ∗ N (m, Σ) = (m −−1/2m∗ , Σ∗

1/2

1/2

−1/2

[Σ∗ ΣΣ∗ ]1/2 Σ∗− I) によって与えられる. 一般に, 多変量分布間の最適輸送問題は解を陽に表すことができないが, 分布のクラスをガウス分布族 G に制限することで, このように最適輸送写像, Wasserstein 距離及び接空間の点との対応を陽に与えることが可能となる.

3. 提案するモデル
F を G × G 上の同時分布とし, (ν1 , ν2 ) ∼ F , νj = N (mj , Σj ), j = 1, 2 とする. また, 各 j = 1, 2 に対し, νj のFréchet 平均を νj⊕ = N (mj⊕ , Σj⊕ ) とする. 以下では, ν1 を説明変数, ν2 を結果変数とするような, G から G への回- 18 -

帰を考える. [1] と同様に空間 (G, dW ) の幾何を用いて, ガウス分布 ν1 及び ν2 をそれぞれ非線形制約のない行列に変換し, 分布値データ間の回帰を行列データ間の回帰に帰着させるというアプローチを取る.
まず, 各 j = 1, 2 に対し, Fréchet 平均 νj⊕ において (G, dW ) の接空間をはり, 分布 νj = Nj (mj , Σj ) を logνj⊕ νj ∈Sd という行列に変換する. ここで, Sd は Sd = {(a, B) ∈ Rd×(d+1) : a ∈ Rd , B ∈ Sym(d)} という非線形制約の無い行列の集合である. そして, 変換後の行列 X = logν1⊕ ν1 と Y = logν2⊕ ν2 の間に, 線形回帰モデル

Y = X, B + E,

E[E|X] = 0を仮定する. ここで, 4 次のテンソル B ∈ Rd×(d+1)×d×(d+1) が回帰パラメータであり, X, B ∈ Rd×(d+1) はcontracted tensor product と呼ばれる行列とテンソルの間の積である. さらに, 分布の次元 d が大きい時は, 回帰係
数 B の要素数は膨大になるため, B は低ランクを持つことを仮定する. 最後に, 空間 (G, dW ) とその接空間の間の対応が一対一になるために, 確率 1 で X, B ∈ logν2⊕ G となることも仮定する.

4. 推定量の構成とその性質について
実際には確率分布が直接観測されることは稀で, 我々は各分布からの離散観測のみが得られる場合が多い. そこで,
以下ではまずガウス分布のペア (ν1i , ν2i ) ∼ F , i = 1, ..., n が潜在的に生成され, その後潜在的な各分布から d 次元ベ
クトル Wjir ∼ νji , r = 1, ..., N が離散観測されるという設定を考える. 回帰パラメータ B が識別されるようにパラメータ空間 Θ をとり, 観測値 Wjir ∼ νji , r = 1, ..., N, i = 1, ..., n, j = 1, 2 に基づいて B を推定する.
推定量の構成法は以下の通りである. まず, 離散観測 Wjir , r = 1, ..., N を用いて, 潜在的な分布 νji の推定量ν̂ji = Nj (m̂ji , Σ̂ji ) を構成し, 推定された分布 ν̂ji , i = 1, ..., N を用いて, 経験 Fréchet 平均 ν̂j⊕ を計算する. 次に, 経験 Fréchet 平均 ν̂j⊕ で接空間をはり, 推定された分布 ν̂ji を X̂i = logν̂1⊕ ν̂1i ∈ Sd , Ŷi = logν̂2⊕ ν̂2i ∈ Sd という行列に変換する. 最後に, 規準関数を M̂n,N (B) = n−1ni=1 Ŷi −X̂i , B 2F とおき, 最小二乗法 arg minB∈Θ M̂n,N (B)

によってパラメータを推定する. ここで,  · F はフロべニウスノルムを表す.
こうして得られる推定量について, 一致性, 収束レート及び漸近正規性に関する以下の結果が示せる.

Theorem 1 真の Fréchet 平均 ν1⊕ , ν2⊕ が既知であると仮定し, その状況で得られる推定量を B̂n,N とする. また, パ
ラメータの真値を B0 とする. この時, いくつかの正則条件のもと, B̂n,N − B0 F = Op (n−1/2 + N −1/4 ) が成立.
さらに, N は n の数列で, N (n) = nq (q > 2) あれば, n → ∞ のとき√

n(vec(B̂n,N ) − vec(B0 )) はある正規分布

N (0, V ) に分布収束する.
発表ではこの他, Calgary weather data を用いた実データの解析例なども紹介した.

参考文献[1] Yaqing Chen, Zhenhua Lin, and Hans-Georg Müller. Wasserstein regression. Journal of the American Statistical Association, pages 1-14, 2021.
[2] Laya Ghodrati and Victor M Panaretos. Distribution-on-distribution regression via optimal transport maps.
Biometrika (to appear), available at arXiv preprint arXiv:2104.09418.
[3] Asuka Takatsu. Wasserstein geometry of gaussian measures. Osaka Journal of Mathematics, 48(4):1005 ‒1026, 2011.

- 19 -

GMANOVA モデルでの仮定緩和のさらなる可能性(報告書)
中京大学 教養教育研究院 永井 勇
n 個の各個体に対して, 全ての個体で測定時点を揃った状態で p 回測定して得られる経時測定データと呼ばれるデータの分析を本講演では考えた．このようなデータの分析の主な目的は,データに隠れている経時変動と呼ばれる時間的な変動を上手く捉えることである. そのための分析には, Pothoﬀ and Roy (1964) で提案された次の一般化多変量分散分析 (GMANOVA)
モデルがよく使われる;Y = 1n μ X  + AΞX  + E.

(1)

ここで 1n は全ての成分が 1 の n 次元ベクトル, Y は各行が各個体で測定して得られた経時測定データからなる n × p 行列, A は各行が各個体の特徴を表す測定時点に無関係な k 個の説明変数のデータからなる n × k 説明変数行列とし, A 1n = 0k (0r は全ての成分が 0 の

r 次元ベクトル) を満たしているとし, X は後述のように, q は解析者が用いる関数によって決まるような, 各行が測定時点の関数からなる p × q 行列であり, これらは既知である. また,

rA = rank(A), rX = rank(X) とする. さらに, μ は q 次元の未知ベクトル, Ξ は k × q の未知行列であり, E は E[E] = 0n 0p , Cov[vec(E)] = Σ ⊗ In の n × p の誤差行列とし, Σ は未知の p × p 正定値行列とする. ここで, A が満たしている A 1n = 0k は, 各説明変数 (各列) ごとにそれぞれ中心化されていることを表している.
このモデルで, p 回の測定時点を t1 , . . . , tp (t1 < t2 < · · · < tp ) とし, 例えば X の i 行目を(t0i , t1i , . . . , tq−1
) とすると, μ や Ξ は t1 , . . . , tp の (q − 1) 次多項式の係数からなるベクトルとi行列をそれぞれ表しており, これらを推定することが測定時点 t1 , . . . , tp の (q − 1) 次多項式を用いて経時変動 E[Y ] を推定することに対応していることを講演した. また, 過剰適合の問題もあるが, より柔軟な関数を X に用いて経時変動を推定することも可能であり, その際は,用いる関数の重み付き和で経時変動を推定することとなり, μ や Ξ は重みの部分に対応していることも講演で触れた. これらのことから分かるように, GMANOVA モデルであるモデル(1) における未知のベクトル μ や行列 Ξ を推定することで経時変動が推定できる.
したがって, この GAMOVA モデルにおいて, 経時測定データからなる行列 Y や各個体の説明変数を各行に持つ説明変数行列 A などから, 未知の μ, Ξ を推定することが, 経時測定
データの分析における目的である経時変動の推定において重要であることから, 本講演ではこれらを rA や rX の様々な状況における推定について着目した. これらの推定は, 次のリスク関数を最小にすることでよく行われている;
R(μ, Ξ|Σ) = tr (Y − E[Y ])}Σ−1 (Y − E[Y ]) ,ここで E[Y ] = 1n μ X  + AΞX  である.

- 20 -

このリスク関数 R(μ, Ξ|Σ) を最小にするような μ と Ξ を求めることで得られるそれぞれの推定量 μ̂, Ξ̂ は, rA = k かつ rX = q であれば, μ̂ = (X  Σ−1 X)−1 X  Σ−1 Y  1n /n とΞ̂ = (A A)−1 A Y Σ−1 (X  Σ−1 X)−1 の形で陽に得られることなどを講演では触れた. しかしながら, n < k などの状況で起きる rA = k かつ rX = q を満たさない場合はこれらの推定量が陽に得られないため, リスク関数の計算を工夫するなどして最小化する必要があることから, 本講演で様々な状況での推定量について講演した.
そこで, rA = k かつ rX < q の仮定の下での推定量を永井 (2022a, 2022b) で提案したことを講演で報告した. また, 多変量線形回帰モデル (モデル (1) で X = Ip としたモデル) において永井 (2021, 2022c) で提案した推定量をモデル (1) へ拡張することで, rA < k かつ

rX = q の仮定の下での推定量を, 本講演ではサブとして提案した. さらに，本講演のメインとして，rA < k かつ rX < q の仮定の下での推定量を提案した. その際に用いた共通するア
イデアにおいて，非常に強い仮定を置くと Koll and von Rosen (2005; Def. 4.1.3) で提案された Extended GMANOVA モデルと対応していることにも触れた．そこで提案した推定量は，一つだけ陽に得られ，他は順番に求まることを述べた．また，メインで提案した推定量の問題点にも触れ，そこを改良するアイデアとして永井 (2019) の手法が用いれると考えられることにも触れた．

GMANOVA モデル (1) において，様々な仮定とそれぞれで提案してきた推定量が以下のようにまとめられることを講演した.

表; rA や rX の仮定と推定手法について
仮定

rA = k & r X = q
rA = k & r X < q
rA < k & r X = q
rA < k & r X < qμ の推定量や Ξ の推定量

それぞれの不偏推定量(X  Σ−1 X)−1 X  Σ−1 Y  1n /n, (A A)−1 A Y Σ−1 X(X  Σ−1 X)−1
永井 (2022a)
永井 (2022b)
本講演でサブとして提案
本講演のメインとして提案
アイデア段階として触れた

引用文献:[1] Kollo, T. & von Rosen, D. (2005). Advanced Multivariate Statistics with Matrices,
Springer.
[2] Pothoﬀ, R. F. & Roy, S. N. (1964) A generalized multivariate analysis of variance model
useful especially for growth curve problems. Biometrika, 51, 313–326.
[3] 永井 勇 (2019) バランス型経時測定データにおける Extended GMANOVA モデルの解釈と新たな推定法, 多様な分野における統計科学に関する諸問題.
[4] 永井 勇 (2021) 高次元小標本における多変量線形回帰モデルでの推定法, 2021 年度統計関連学会連合大会.
[5] 永井 勇 (2022a) GMANOVA モデルにおける新たな推定方法とその解釈, . 多様な分野における統計科学の理論とその応用.
[6] 永井 勇 (2022b) GMANOVA モデルにおける新たな推定方法と解釈, 大規模複雑データの理論と方法論〜新たな発展と関連分野への応用〜.
[7] 永井 勇 (2022c) 説明変数がランク落ちしている状況での多変量線形回帰における不偏推
定量, 2022 年度統計関連学会連合大会.

- 21 -

多変量正規母集団における条件付き独立性検定について
松内 直輝（神戸大学大学院理学研究科）
首藤 信通（神戸大学大学院理学研究科）
本報告においては, 多変量正規母集団における条件付き独立性検定問題について


議論した. 特に, 多変量正規性を持つ p 次元確率変数ベクトル X = (X 1 , X 2 , X 3 )

に対し, 母集団から得られた標本ベクトルを基にして X 1 が与えられた下で, X 2と X 3 の条件付き独立性を帰無仮説とする仮説検定を行うための尤度比検定を構成した. ただし, X i は X の pi 次元分割ベクトル, p = p1 + p2 + p3 である.
本報告では, まず X 1 が与えられた下で, X 2 と X 3 が条件付き独立であることと同値となる条件を求めた. p 変量正規母集団 Np (μ, Σ) の下では, 平均べクトル μ,
分散共分散行列 Σ について, X の分割と対応する分割を⎞⎛⎞⎛
Σ11 Σ12 Σ13
μ1
⎟
⎟
⎜
⎜
⎟
⎜
⎟
⎜
μ = ⎜ μ2 ⎟ , Σ = ⎜ Σ21 Σ22 Σ23 ⎟
⎠
⎠
⎝
⎝
μ3
Σ31 Σ32 Σ33とする. 本報告で扱う仮説検定問題の帰無仮説と

Σ23·1 ≡ Σ23 − Σ21 Σ−1
11 Σ13 = O23
が同値であることを示すことができた. ここに, μi は μ の pi 次元分割ベクトル, Σij
は Σ の pi × pj 分割行列である.
また, 以下のようなパラメータの変換
⎛
⎞
⎛
⎞
Σ11
Σ−1
Σ
Ψ
Ψ
Ψ
12
13
11 12
⎜
⎟
⎜ 11
⎟
Σ−1
(12) Σ(12)3 ⎟⎜ −1 −1⎜⎟⎜ Σ21 Σ11
⎟ → ⎜Ψ21 Ψ22 Ψ23 ⎟ ,
Σ22·1
⎝
⎠
⎝
⎠
−1
Σ3(12) Σ(12)Σ33·(12)Ψ31 Ψ32 Ψ33
⎞
⎛
⎛
⎞
μ1
⎟
⎜
η
⎟
⎜
⎟
⎜ 1 ⎟
⎜
μ2 − Ψ21 μ1⎟⎜⎜⎛ ⎞ ⎟ → ⎜ η2 ⎟
⎟
⎜
⎟
⎜
⎝
⎠
μ
⎟
⎜
1
⎝
⎠
⎠
⎝ μ3 − Ψ3(12)η3μ2- 22 -

を考えると X 1 が与えられた下で, X 2 と X 3 の条件付き独立性に関する仮説検定
問題はH0 : Ψ23 = O23

vs. H1 : Ψ23 = O23と書き換えることができる.
具体的に, Π : Np (μ, Σ) から得られた N 個の p 次元標本ベクトル y 1 , . . . , y N が得られたとしてこの仮説検定問題における尤度比検定統計量を構成すると尤度比
検定統計量は
−2 log λ = N log(det(S33·1 )) − log(det(S33·(12) )−1
S23·1 ))
= N log(det(S33·1 )) − log(det(S33·1 − S32·1 S22·1となることが示された. ここに, Sij は

1 
y ,
ȳ =
N j=1 j
N



1 
S=
y j − ȳ y j − ȳ ,
N j=1
N

−1
としたとき, S の pi × pj 分割行列である. また, Sij·1 = Sij − Si1 S11
S1j ,

⎞

⎛
S(12) = ⎝S11 S12
S21 S22⎞⎛⎠ , S(12)3 = ⎝

S13
S23

⎠ , S33·(12) = S33 − S3(12) S −1 S(12)3
(12)

である.
最後に数値実験を行い, 本報告で議論している条件付き独立性の下で擬似的に発生させたデータセットから尤度比検定統計量を求め, シミュレーションで与えられる上側 100α% 番目の値と, 自由度 p2 p3 のカイ二乗分布の上側 100α%点の比較を行った. 数値実験の結果から p や p1 の値よりも p2 p3 の値が強く第 1 種過誤の精度に影響を与えることを確認した.
今後の課題としてより少ない標本でも検定の精度が保てるように分布の修正を行うことや, より緩い仮定の下での仮説検定の構成を行うことなどが考えられる.

- 23 -

ଟॏੵ෼ϒϥ΢ϯӡಈͷ౷‫ֶܭ‬తԠ༻ʹ͍ͭͯాத উਓʢҰ‫ڮ‬େֶ໊༪‫ڭ‬तʣ

[0,1] ্Ͱఆٛ͞Εͨϒϥ΢ϯӡಈ {W (t)} ͓ΑͼϑϥΫγϣφϧɾϒϥ΢ϯӡಈ {BH (t)}ͷ 2 ࣍൚ؔ਺ʹؔ࿈͢Δ౷‫ͯؔ͠ʹྔܭ‬ɼҎԼͷ (1)ɼ(2)ɼ(3) ͷใࠂΛͨ͠ɽ(1) g ॏੵ෼աఔ
Fg (t) = t  t100···

 tg−10

W (tg ) dtg dtg−1 · · · dt1

ʹରͯ͠ɼ

Vg =
ͷ෼෍ͷಛੑؔ਺͸ɼ 1
0(g = 1, 2, ...)

Fg2 (t) dtE(eiθVg ) = Dg (2iθ)

−1/2Ͱද͞ΕΔ͜ͱΛࣔͨ͠ɽ͜͜ͰɼDg (λ) ͸ɼؔ਺ Cov(Fg (s), Fg (t)) ͷϑϨουϗϧϜ
ߦྻࣜ (FD) Ͱ͋ΓɼFD Λಋग़͢Δํ๏Λઆ໌ͨ͠ɽg = 1, 2, 3 ͷ৔߹ͷ FD ͸࣍ͷΑ͏ʹͳΔɽD1 (λ) =D2 (λ) =ω =1(1 + cos λ1/4 cosh λ1/4 ),2
1
2(1 + cos λ1/6 + cos λ1/6 ω + cos λ1/6 ω 2 ) + cos λ1/6 cos λ1/6 ω cos λ1/6 ω 2 ,
9 √
1 + 3i,21
3 cos a cos b cos c cos d + 2 cos a cos b + cos b cos c + cos c cos d16+ cos d cos a + cos a cos c + cos b cos d + 3√
2
+
sin a sin b (1 + cos c cos d) + sin b sin c (1 + cos d cos a)16+ sin c sin d (1 + cos a cos b) − sin d sin a (1 + cos b cos c) ,1+ia = λ1/8 , b = aω, c = aω 2 , d = aω 3 , ω = √ .
2D3 (λ) =·ͨɼg ͕େ͖͘ͳΔ৔߹ʹ͸ɼ࠷খ‫ݻ‬༗஋ʢDg (λ) = 0 ͷ࠷খࠜʣ͚ͩΛ࢖ͬͨ෼෍Ͱ΋Α͍ۙࣅΛ༩͑Δ͜ͱΛࣔͨ͠ɽ(2) ଟॏ୯Ґࠜ౷‫ྔܭ‬F 2 (1)/2
Rg =  1 g
F 2 (t) dt(g = 1, 2, ...)

0ͷ෼෍ʹ͍ͭͯ΋ؔ࿈͢Δಛੑؔ਺Λ‫ٻ‬Ίɼີ౓ؔ਺Λ਺஋ੵ෼ʹΑΓ‫ٻ‬Ίͨɽ෼෍ؔ
਺͸ɼ
  1

P (Rg < x) = P x

01
Fg2 (t) dt − Fg2 (1) > 0
2- 24 -

= P

 1  1
0

1 1
+
=
2 π

0



Kg (s, t; x) dW (s) dW (t) > 0

 ∞
0






−1/2
1
Im Dg (2iθ; x)dθθʹΑΓ‫͢ࢉܭ‬Δ͜ͱ͕Ͱ͖Δɽ͜͜ͰɼDg (λ; x) ͸ɼ
 1

g
g
1
1
Kg (s, t; x) =
x
(u
−
s)(u
−
t)
du
−
(1
−
s)(1
−
t)
(g!)2
2
max(s,t)ͷ FD Ͱ͋ΔɽFD ͷಋग़ʹ͸ɼੵ෼ํఔࣜ

f (t) = λ

 1
0

Kg (s, t; x)f (s) ds͕ɼඍ෼ํఔࣜ

f (2g+2) (t) + (−1)g λxf (t) = 0
͓Αͼ 2g + 2 ‫ݸ‬ͷ‫ڥ‬ք৚݅

f (1) = f  (1) = · · · = f (g−1) (1) = 0,

f (g) (1) =(−1)g−1 λ  1(1 − s)g f (s) ds,
2g!
0f (g+1) (0) = f (g+2) (0) = · · · = f (2g+1) (0) = 0.
ͱಉ஋ͱͳΔ͜ͱΛ༻͍ͨɽ͜͜ͰɼE(Rg ) = g + 1 ͱ͍͏‫ڵ‬ຯਂ͍݁Ռ͕ಘΒΕͨɽ(3) ௨ৗͷ୯Ґࠜ౷‫ͱྔܭ‬ଟॏ୯Ґࠜ౷‫ྔܭ‬ͷؒʹհࡏ͢ΔϑϥΫγϣφϧ୯Ґࠜ౷‫ྔܭ‬
B 2 (1)/2
QH =  1 H
2BH(t) dt0ͷ෼෍ʹ͍ͭͯߟ࡯ͨ͠ɽ͜͜ͰɼBH (t) ͸ϑϥΫγϣφϧɾϒϥ΢ϯӡಈ (fBm) Ͱ͋ΓɼfBm ͸ฏ‫ ۉ‬0 ͷਖ਼‫ن‬աఔͰɼ‫ڞ‬෼ࢄؔ਺͕

1
Cov(BH (s), BH (t)) = [s2H + t2H − |s − t|2H ]2Ͱ༩͑ΒΕΔɽSH =SH ͷۙࣅͱͯ͠ɼ
TH =

 1
0

 1
0

2
BH
(t) dt ͓Αͼ QH ͷ෼෍ಋग़͸ະղܾ໰୊Ͱ͋Δɽ͜͜Ͱ͸ɼ2CH(t) dt,CH (t) =2(1 − H) t2H−1ΛఏҊͨ͠ɽTH ʹؔ࿈ͨ͠ FD ͸ɼDH (λ) = Γ(1 − ν)J−ν (η)

 t
0u1/2−H dW (u)

 
−ν

η
2



ͱͳΔɽ͜͜ͰɼJ−ν (η) ͸ୈ 1 छϕοηϧؔ਺ɼ·ͨɼη = 2(1 − H)λ/(H + 1/2),(2H − 1/2)/(H + 1/2) Ͱ͋Δɽ·ͨɼQH ͷۙࣅͱͯ͠ɼ

C 2 (1)/2
RH =  1 H
2CH(t) dt0ΛఏҊͯ͠ɼͦͷ෼෍ಛੑʹ͍ͭͯ΋ߟ࡯ͨ͠ɽ

- 25 -

ν=

On Asymptotic Distribution in Martingale Convergence ofSupercritical Branching Processes with Poissonian oﬀsprings
Bat-Erdene A.* , Kawasaki S.* , Li J.** , Altantsetseg E.****

Faculty of Science and Engineering, Iwate University** School of Computer Science, Chongqing University*** School of the Engineering and Applied Sciences, National University of Mongolia

Abstract
A branching process is a mathematical model of Erdos-Renyi random graphs. For a normalized
process, a martingale convergence theorem holds. However, its asymptotic distribution has notbeen known so far. We propose a characterization of the distribution via analysis on a functional
equation for the Laplace transform of the distribution. It turns out that a numerical analysis mostly
coincides well with the theoretical characterization. Applications of the Branching process include
a biological population, nuclear chain reactions, and the spread of computer software viruses incommon. Mathematical models of these applications play a central role in ﬁguring out the mainprocess and predicting future extensions. Our funding works in common cases and provides an
explanation for the characterization of the processes.

Problem. Let Wn = Zn /λn , n = 0, 1, 2, · · · , where Zn is the supercritical branching process withthe mean λ which is greater than 1. Then, {Wn } is known [2] to form a martingale and converges toa random variable W∞ on R+ = [0, ∞) a.s. as n → ∞.
However, not much is known about the distribution of W∞ so far in previous studies. The paper
aims to characterize the distribution by analyzing a functional equation that holds for the Laplacetransform of the distribution.

Analysis and result. W∞ is known to have a density function on (0, ∞) [2, Corollary 12.1], which

we denote by w : (0, ∞) → R+ , while W∞ has a point mass at the origin [2, Theorem 6.2]. Thus, we
may write W∞ ∼ ηδ(x) + w(x) on x ∈ R+ and δ(x) being the Kronecker delta function.
Let ϕ : R+ → R+ be the Laplace transform of W∞ ,e−ux w(x) dx.
(1)
ϕ(u) =(0,∞)

Mclaurin’s expansion was used to approximate the function ϕ(u) described in below.
(2)

ϕ(u) =

 ϕ(n) (0)
n∈N0

n!

un

Then, upon calculating ϕ(n) (0), n ∈ N recursively. To the evaluate ϕ(n) (0), an approximation of Kn
as we deﬁned in below is necessary(3)

ϕ(n) (0) = −Kn (−ρ)n−1 ,n∈Nfor appropriate functions Kn = Kn (λ).
We write 1 + λ1 + · · · + ( λ1 )n−2  ρn−1 below. We note that ρn −→ ρ. Through the recursive- 26 -

relationn − 1 −(n−1−l)
λKn−l (λ)Kl (λ)
l

n−1


1

Kn =

ρn−1

l=1
n−1

=

(4)

(n − 1)!  −(n−1−l) Kn−l (λ)
Kl (λ)λ·· (n − l),ρn−1(n − l)!
l!
l=1

In the result stage, we propose an approximation function Kn (λ) as follows:
Theorem 1. For arbitrarily ﬁxed λ ∈ [1, ∞),−(an+b) · 1 + O λ−n(5)Kn  cn × (n − 1)! λ(λ + 1)

,

as

l→∞

for some cn > 0.
A derivative ϕ (u), instead of ϕ(u) has considered. That is,ϕ (u) =(6) ϕ(n+1) (0)
n∈N0

n!

un = −

 Kn+1 (λ)
n∈N0

n!

(−ρu)n .

Applying the expression (5) to (6), we have that the main term of ϕ (u) is given by
const. ×

 


−a n
− ρ λ(λ + 1)
u = const. ×n∈N0

u + ρ−11
a .
λ(λ + 1)



Now we recall that ϕ (u) is the Laplace transform of −x w(x). Thus, by the inverse Laplace transformof ϕ (u), we found that w(x) has the corresponding main component of the probability density function
given by 
aconst. × x−1 exp − λ(λ + 1) x ,x ∈ (0, ∞).

Figure 1: a. The function Kn (λ), its relation and values of an equation pn = an + b. b. Asymptotic
linearity, pn = an + b at some ﬁxed λReferences[1] Hofstad, R. (2016). Random Graphs and Complex Networks, Cambridge Univ Press.
[2] Athreya, K. B., Ney, P. E. (2004). Branching Processes, Dover publ.
[3] Harris, T. E. (1963). The Theory of Branching Processes, Dover publ.

- 27 -

Forest Construction of Gaussian and Binary Variables based on WBIC
Ashraful Islam, Joe Suzuki
Graduate school of Engineering Science, Osaka University, Japan

Mutual information (MI) is a metric that determines the association between two random variables bymeasuring the amount of information that one variable holds about the other. It quantiﬁes the dependency
between them. A higher mutual information value indicates a stronger relationship between the randomvariables. It is related to the concept of entropy and is used in various ﬁelds such as information theory,
statistics, and machine learning. Estimated joint probabilities from observed samples in each variable
category combination are used to calculate MI between discrete variables. But the conventional method
is ineﬃcient for estimating the MI of a mixture of discrete and continuous random variables because theconditional probabilities for discrete variables given continuous variables cannot be determined. In thispaper, we examine a new MI method that can accommodate a mixture of discrete and continuous randomvariables. The estimation of free energy for new MI is not easy, but we can use WBIC to calculate thecorresponding free energies. After that, the Chow-Liu algorithm can be modiﬁed to take into accountthe new MI in order to build a forest rather than spanning trees.
This section will outline the paper’s most signiﬁcant ﬁndings and conclusions. Chow and Liu (1968)
considered estimating mutual information between two discrete variables. When only one of the variables
is Gaussian and the rest are discrete, Edwards, De Abreu, and Labouriau (2010) recalculated mutual
information using the ANOVA model. If X is a Gaussian random variable, and Y and Z are discretevariables, getting maximum likelihoods of the graphical models such as X − Z and Y − X are simple,
but obtaining maximum likelihoods for other types of graphical models as Y − X − Z by maximizing
observational mutual information estimates is diﬃcult(Suzuki, 2017). The state-of-the-art avoids thisproblem by restricting the class to forests that separate discrete and continuous nodes.
On the other hand, the Bayesian method makes the assumption that there is inﬁnite number ofhistograms, but the number of histograms that are required for accuracy is unclear to us. Each cluster
also has a large sample size variance. Adjusting parameters is diﬃcult when many clusters have nosample. The proposed approach succeeds where the other methods fail in capturing the relationshipbetween multiple gene expressions and SNPs. Some previous works that relate to our research are as
follows:
Suzuki (1993) considers MDL method and obtains total lengths for all the variables, and estimatesthe mutual information as:(α − 1)(β − 1)
log2 n(1)2n
For large n, Jn ≤ 0 indicates the independence of X and Y (Suzuki, 2012).
Edwards et al. (2010) assumed about a situation where some random variables are continuous andsome are Gaussian variables. Assume, X (1) , X (2) , . . . , X (N ) are Gaussian. And, mutual information:
Jn = In −1(2)
In (i, j) = − log2 1 − ρ̂2 ,2
To use the mutual information estimator, Edwards et al. (2010) imposed the limitation that noGaussian variables could be intermediate between the discrete ones.
Suzuki (2017) took into account the Dirichlet distribution with the proportionality x θ(x)a(x)−1 ,and derived the formula for mutual information as follows:
Jn =Qn (xn , y n )1
log2 n n n nn
Q (x ) Q (y )(3)Suzuki (2017) has found that −log2 Qn (xn , y n ) is nothing but the free energy of the mixture of xn andy . In the ﬁeld of data science, free energy refers to the amount of uncertainty or surprise present in asystem or dataset. In statistical mechanics and information theory, free energy is employed to quantify
n- 28 -

the disorder or randomness of a system. But obtaining free energy is not easy. In this circumstance,
WBIC (Watanabe, 2021) is a way to ﬁgure out free energy.
The WBIC algorithm estimates the free energy of a system using Bayesian statistics, taking into account the uncertainty in themodel’s parameters. Again, we use the Kernel Hilbert–Schmidt independence criterion (HSIC) to verify the Mutual Information we
get from the new mutual information method. In the context ofstatistical inference, the HSIC was developed to assess how muchthe embedding of the Hilbert space depends on the underlying distribution. To obtain the free energy from the mixture of discreteand Gaussian, we use stan. Stan speciﬁes statistical models probabilistically. Stan uses Markov Chain Monte Carlo methods like theNo-U-Turn sampler for continuous variable Bayesian inference.
Let, for large n, X and Y are discrete and Gaussian randomvariables, respectively, and their joint pdf,
8

11

10

12

6

4

5

1

2

3

7

9

f (x, y) =

2

i=1

px (1 − p)1−x . √

1
1
exp{
(y − μi )2 }
2σi
2πσi

Figure1:
Forest’P.sojae.survey’ datafor

And we propose the new mutual information:
Qn (X n , Y n )
1
e−Fxy
1
1
log n n=log
= (Fx +Fy −Fxy )n
Q (X ).Qn (Y n )ne−F x .e−F yn(4)
We calculate the MI from equation 4 using simulated data fora mixture of discrete and Gaussian variables.
When we use n=500, we get the following results for the mutual information: -0.008828827, HSIC: 0.000759246, and p-value:
0.2397602. These ﬁndings both conﬁrm and imply that the variables X and Y are free to move around without inﬂuence from one
another.
We can see that the suggested mutual information estimator
Figure 2: Forest using simulated
works
for all types of variables: discrete, continuous, or a mixdata
ture of discrete and continuous variables. After calculating themutual information from the new method, we applied it to the
Chow-Liu algorithm to ﬁnd the forest. Figure 1 shows the forest
for P.sojae.survey data with 12 variables; we can see that the ﬁrst variable is independent of the others.
That is why it is separated from the other variables. Figure 2 shows a forest with 75 simulated observations with a mixture of discrete and continuous variables with 100 samples each. Therefore, dependingon the previous results and the ﬁgures, the current method can construct a forest from a mixture of
discrete and continuous nodes, depending on the independence of the variables. In future work, we will
apply the proposed mutual information method to multiple gene expressions and SNPs (single-nucleotide
polymorphism) data.
Jn =References
Chow, C., & Liu, C. (1968). Approximating discrete probability distributions with dependence trees.
IEEE transactions on Information Theory, 14 (3), 462–467.
Edwards, D., De Abreu, G. C., & Labouriau, R. (2010). Selecting high-dimensional mixed graphicalmodels using minimal aic or bic forests. BMC bioinformatics, 11 (1), 1–13.
Suzuki, J. (1993). A construction of bayesian networks from databases based on an mdl principle. In
Uncertainty in artiﬁcial intelligence (pp. 266–273).
Suzuki, J. (2012). The bayesian chow-liu algorithm. In The sixth european workshop on probabilistic
graphical models (pp. 315–322).
Suzuki, J. (2017). A novel chow–liu algorithm and its application to gene diﬀerential analysis. International Journal of Approximate Reasoning, 80 , 1–18.
Watanabe, S. (2021). Waic and wbic for mixture models. Behaviormetrika, 48 (1), 5–21.

- 29 -

多次元空間における埋め込み 1 次元曲線の同時信頼領域

はじめに:

千葉大・融合理工学府

山添 滉弥

千葉大・理学研究院

内藤 貫太

道路交通網上の交通事故や，脳神経細胞網における樹状突起スピンの分布といった，枝分かれした

1 次元曲線上で観測されるデータをネットワークデータを呼ぶ．ネットワークデータに対して，与えられたネッ
トワーク上における密度推定の手法についてはすでに研究されている (Liu and Ruppert, 2021; McSwiggan

et al., 2017)．しかしながら，ネットワークを構成する 1 次元曲線の同時信頼領域をデータから構築する方法については，いまだ研究がないようである．本研究ではこの問題を多次元空間に埋め込まれた 1 次元曲線の同時信頼領域を構築する問題として捉え，ノンパラメトリック回帰のアプローチから取り組んだ．
問題:

I = [a, b] ⊂ R を閉区間とし，ϕi : I → R (i = 1, . . . , d) を滑らかな関数とする．このとき，Rd に埋め込まれた 1 次元曲線を⎤
ϕ1 (t)⎥⎢
ϕ(t) = ⎣ ... ⎦ ∈ Rd⎡

ϕd (t)
で表す．確率ベクトル Y ∈ Rd と共変量 T ∈ I の組 (T, Y ) からの n 個の観測値 (T1 , Y1 ), . . . , (Tn , Yn ) に基づ

き 1 次元曲線 ϕ の同時信頼領域を構築したい．つまり，Mϕ = ϕ(t) ∈ Rd | t ∈ I ⊂ R としたとき，十分小さい α > 0 に対して，

P (D ⊃ Mϕ ) ≥ 1 − αとなるような有界閉領域 D ⊂ Rd の構築を目指す．
設定: 1 次元曲線の推定の関連研究である Hastie and Stuetzle (1989) を参照しつつ，モデルの仮定を行う．Rdに埋め込まれた 1 次元曲線 ϕ に対して，任意の点 t における勾配ベクトル ϕ (t) を ϕ (t) = [ϕ1 (t) · · · ϕd (t)]Tで表す．また，n1 (t), . . . , nd−1 (t) を {ϕ (t)/ ϕ (t) , n1 (t), . . . , nd−1 (t)} が Rd の正規直交基底となるようにとる．さらに，T を I 上の密度 fT に従う確率変数とし，V = [V1 · · · Vd−1 ]T を原点中心半径 r の (d − 1) 次
元球 Brd−1 上の密度 fV に従う確率ベクトルとする．
i.i.d.

(T1 , V1 ), . . . , (Tn , Vn ) ∼ fT × fV と E[V1 ] = 0 を仮定する．確率ベクトル Y の d 次元標本 Y1 , . . . , Yn として，

Yi = ϕ(Ti ) + N (Ti )Viを考える．ただし，Yi = [Yi1 · · · Yid ]T , Vi = [Vi1 · · · Vid−1 ]T , N (t) = [n1 (t) · · · nd−1 (t)] である．
同時信頼領域の構築:

共変量 T = [T1 · · · Tn ]T と Y の第 j-成分のデータ Yj = [Y1j · · · Ynj ]T との組(T , Yi ) (i = 1, . . . , d) に基づき，t における ϕi (t) の局所線形推定量 ϕi (t)(i = 1, . . . , d) が得られる．これらを並べることにより ϕ(t) の推定量 ϕ(t) が構築される．さらに ϕ(t) を用いてその勾配ベクトルを求める：⎤
ϕ1 (t)⎥⎢
ϕ(t) = ⎣ ... ⎦ ,⎡ϕ (t) =ϕd (t)- 30 -dϕ(t).
dt局 所 線 形 推 定 量 に つ い て は Wand and Jones (1995); Fan and Gijbels (1996) な ど を 参 照 さ れ た い ．

n1 (t), . . . , nd−1 (t) を {ϕ (t)/ ϕ (t) , n1 (t), . . . , nd−1 (t)} が Rd の正規直交基底となるようにとり，ϕ (a+)d−1+ N (a)v 0 ≤ R ≤ r, v ∈ B√r2 −R2 ,Da (r) = ϕ(a) − R ϕ (a+)

DJ (r) = ϕ(t) + N (t)v t ∈ J, v ∈ Brd−1 ,ϕ (b−)d−1√
Db (r) = ϕ(b) + R + N (b)v 0 ≤ R ≤ r, v ∈ B r2 −R2 ,ϕ (b−)

N (t) = n1 (t) · · · nd−1 (t)とする．3 つの排反な領域を用いて同時信頼領域 D(r) ⊂ Rd をD(r) = Da (r) ∪ DJ (r) ∪ Db (r)として定義する (図 1 参照)．本講演では ϕ(t) と D(r) の理論的性質および，信頼係数 1 − α から D(r) の r を決定する方法について解説を与えた．また，共変量 T がデータとして得られていない場合において，Y1 , . . . , Yn
から T1 , . . . , Tn を作成する方法についても報告した．さらに，提案領域が同時信頼領域として機能していることを確認するシミュレーション実験の結果や，実データへの適用例についても報告した．
1.6

1.4

y

1.2

1.0

0.8

0.6
0.6

図1

0.8

1.0

x

1.2

1.4

1.6


 から作られた同時信頼領域 D(r)推定曲線 ϕ(t)の例

参考文献
Fan, J. and Gijbels, I. (1996). Local Polynomial Modelling and Its Applications: Monographs on Statistics
and Applied Probability 66. Taylor & Francis.
Hastie, T. and Stuetzle, W. (1989). Principal curves. Journal of the American Statistical Association, 84,
502–516.
Liu, Y. and Ruppert, D. (2021). Density estimation on a network. Computational Statistics & Data
Analysis, 156, 107128.
McSwiggan, G., Baddeley, A., and Nair, G. (2017). Kernel density estimation on a linear network. Scandinavian Journal of Statistics, 44, 324–345.
Wand, M. P. and Jones, M. C. (1995). Kernel Smoothing. Chapman & Hall.

- 31 -

多標本問題における 2-step 単調欠測データの下での平均ベクトルの同等性検定の検出力について
神戸大学・理・院
神戸大学・理・院

勝又
首藤

真
信通

本報告では, k 個の母集団 Π(g) : Np (μ(g) , Σ) (g = 1, . . . , k) から, 2-step 単調欠測
データが得られた下で帰無仮説 H0 : μ(1) = · · · = μ(k) に対する尤度比検定統計量を導出し, その第 1 種の過誤や検出力について数値的評価を与えた.
(g)
具体的には
, 第 g 母集団
Π(g) (g = 1, . . . , k) から N1 個の p 次元標本ベクトル


(g,1)

xj



(g,1)
= x(g,1)
, x2j
1j





(g)

(g)

(g,2)

(j = 1, . . . , N1 ), N2 個の p1 次元標本ベクトル x1j (j =(g)

1, . . . , N2 ) が観測されたと仮定し, 尤度比検定統計量が

k



1  1  (g,1)
(g,2)
(g,1)
(g,2)
−1
−2 log Λ =N log 1 +
+
x
+
x
x
x
Δ̂
1·
1·
1·
1·
11
N g=1 N (g)
1  (·,1)
(·,2)
(·,1)
(·,2)
−1
− 2 x1· + x1·
Δ̂11 x1· + x1·
N





k
(g,1)
(·,1) (·,1)

 (·,1) −1 x̄(g,1)− x̄1·−x̄(g) x̄1·
1·
1·
N1
W+ N1 log 1 +
(g,1)
(·,1)
(g,1)
(·,1)
−
x̄
−
x̄
x̄
x̄
2·
2·
2·
2·
g=1


k

 
−1 


(g)
(g,1)
(·,1)
(·,1)
(g,1)
(·,1)
N1 x̄1· − x̄1·W11
x̄1· − x̄1·
− N1 log 1 +g=1(g)Ni(g) (g,i) (g,i)(g,i)(g)−1 Ni(g,i)
となることを示した. ここで, x1· = j=1
x1j , x̄1· = Ni
j=1 x1j , N1(g)

(·,1)
(g,1)
(·,1)
(·,1)
x̄1· = N1−1 kg=1 x1· , Δ̂11 = N −1 (W11 + W (·,2) ), W11 = kg=1 j=1
 N2(g) (g,2)
(g,1)
(g,1)
(g,1)
(g,1)
(g,2)
(g,2)
(g,2)
(x1j − x̄1· )(x1j − x̄1· )
(x1j − x̄1· )(x1j − x̄1· ) , W (·,2) = kg=1 j=1(g) (g)
(g,1)
(g,2)
(g,1)
(g,2)
(g)
(g)
+ kg=1 N1 N2 /N (g) (x̄1· − x̄1· )(x̄1· − x̄1· ) , N (g) = N1 + N2 ,



(g)
(g)
N = kg=1 N (g) , N1 = kg=1 N1 , N2 = kg=1 N2 である.
また, k = 3, p = 6, p1 = 2, 3, 4, 5 の場合における上記の尤度比検定統計量について, 自由度 p(k − 1) = 12 のカイ二乗分布の上側 5% 点を棄却限界値として用いた
(1)
(2)
(3)
場合の仮説検定方式と, 同様の設定の下で N2 = N2 = N2 = 0 における尤度比
検定統計量について同じ棄却限界値を用いた場合の仮説検定方式の第 1 種の過誤を比較した. 表 1 はその結果の一部であり, 前者の第 1 種の過誤は ASL(miss), 後者の第 1 種の過誤は ASL(comp) で表す.

- 32 -

表 1: k = 3, p = 6, p1 = 2 における第 1 種の過誤の比較
(1)
(2)
(2)
(3)
(3)
N2
N1
N2
N1
N2
ASL(miss) ASL(comp)
10
10
10
10
10
10
0.085274
0.089827
20
10
20
10
20
10
0.062790
0.063582
30
10
30
10
30
10
0.057406
0.057700
40
10
40
10
40
10
0.054916
0.055446(1)
N1

数値実験の結果, 完全データのみで尤度比検定を構成するよりも欠測データを利用して尤度比検定を構成する方が, 第 1 種の過誤がより正確に制御されることを確認した. また, サンプルサイズ N が小さいとき, p1 が大きくなるほど, 第 1 種の過誤がより正確に制御されることを確認した.
次に, Σ = Ip の場合において, 欠測パターンに対応するマハラノビス二乗距離を
(1)
(2)
(1)
(2)
それぞれ δi2 = (μi − μi ) (μi − μi ) (i = 1, 2) とする. これらのマハラノビス
二乗距離を変化させたときに k = 3, p = 6, p1 = 2, 3, 4, 5 における検出力が p1 の影
響を受けるかについて確認した. 表 2 はその結果の一部である.
(1)

(2)

(3)

(1)

(2)

(3)

表 2: δ12 = (0.0)2 , N1 = N1 = N1 = N2 = N2 = N2 = 10 における検出力
の比較
δ12
δ22
p1 = 5p1 = 4
p1 = 3p1 = 222(0.0) (0.1) 0.072707 0.079148 0.083119 0.087148(0.0)2 (0.5)2 0.122120 0.128863 0.134764 0.138954(0.0)2 (1.0)2 0.311521 0.320074 0.328623 0.336122(0.0)2 (1.5)2 0.621688 0.632784 0.642896 0.652537(0.0)2 (2.0)2 0.875459 0.883849 0.891225 0.898671(0.0)2 (2.5)2 0.977481 0.980229 0.982739 0.985267(0.0)2 (3.0)2 0.997588 0.998116 0.998584 0.998942

数値実験の結果, サンプルサイズ N が小さいとき, p1 が小さくなるほど検出力が高くなる傾向を確認することができた. また, 欠測データでも同様にサンプルサイ
ズ N が増加するにつれて, 検出力が高くなる傾向を確認した.
今後の課題として, サンプルサイズ N が小さい状況下でも, ほぼ正確な仮説検定となるような分布の修正を行うことなどが挙げられる. また, 本報告で提案する仮説検定では平均ベクトルの差がどの群間にあるかを推測することができないので,これを行うには欠測データの下での平均ベクトルに関する多重比較法の構成について考える必要がある.

- 33 -

Subexponentiality of densities of inﬁnitely divisible distributions
on the whole real line
Muneya Matsui

Abstract
We show the equivalence of three properties for an inﬁnitely divisible distribution:
the subexponentiality of the density, the subexponentiality of the density of its
Lévy measure and the tail equivalence between the density and its Lévy measuredensity, under monotonic-type assumptions on the Lévy measure density. The key
assumption is that tail of the Lévy measure density is asymptotic to a non-increasing
function or is almost decreasing. Our conditions are natural and cover a rather wideclass of inﬁnitely divisible distributions. Several signiﬁcant properties for analyzing
the subexponentiality of densities have been derived such as closure properties of[ convolution, convolution roots and asymptotic equivalence ] and the factorization
property. Moreover, we illustrate that the results are applicable for developing
the statistical inference of subexponential inﬁnitely divisible distributions which are
absolutely continuous.

Introduction
Let f, g be probability density functions on R and denote by f ∗ g the convolution of f and g: ∞f ∗ g(x) =f (x − y)g(y)dy,−∞∗nand denote by f the nth convolutions with itself. Throughout the paper, for functions α, β : R →R+ , α(x) ∼ β(x) means that limx→∞ α(x)/β(x) → 1. We study the following characteristics fordensities.

Deﬁnition 0.1. (i) f is (right-side) long-tailed, denoted by f ∈ L, if there exists x0 > 0 such thatf (x) > 0, x ≥ x0 and for any ﬁxed y > 0 f (x + y) ∼ f (x).
(ii) f is (right-side) subexponential on R, denoted by S, if f ∈ L and f ∗2 (x) ∼ 2f (x).
(iii) f with dist. F is weakly (right-side) subexponential on R, denoted by S+ , if f ∈ L and thefunction f+ (x) = 1R+ (x)f (x)/F (0), x ∈ R is subexponential, i.e. f+ ∈ S. Here F (x) = 1 − F (x).
Deﬁnition 0.2. (i) We say that a density f : R → R+ is asymptotic to a non-increasing function(a.n.i. for short) if f is locally bounded and positive on [x0 , ∞) for some x0 > 0, and(0.1)

sup f (t) ∼ f (x)
t≥x

and

inf

x0 ≤t≤xf (t) ∼ f (x).

(ii) We say that a density f : R → R+ is almost decreasing (al.d. for short) if there exists x0 > 0and K > 0 such thatf (x + y) ≤ Kf (x),for all x > x0 , y > 0.

Notice that the al.d. property includes the a.n.i. property, and the latter is satisﬁed by theregularly varying functions with negative indices.
We will investigate properties of the above sort, particularly on inﬁnitely divisible distributions
μ on R. The characteristic function (ch.f.) of μ is ∞1μ(z) = exp(eizy − 1 − izy1{|y|≤1} )ν(dy) + iaz − b2 z 2 ,
(0.2)
2
−∞
∞
where a ∈ R, b ≥ 0 and ν is the Lévy measure satisfying ν({0}) = 0 and −∞ (1 ∧ x2 )ν(dx) < ∞.
Throughout this paper, we always assume that the Lévy measure ν of μ has a density, and wedenote by ID(R) the class of all inﬁnitely divisible distributions on R.

- 34 -

Main contents
We separate the cases depending on weather ν(R) < ∞ or ν(R) = ∞. Note that we use notation galso for the (non-proper) density of a Lévy measure. The followings are the results for the absolutely
continuous case (ν(R) = ∞).
Theorem 0.3. Let μ ∈ ID(R) with ν(dx) = g(x)dx such that ν(R) = ∞. Let f0 (x) be a density ofμ0 ∈ ID(R) with a = b = 0 and ν(dx) = 1{|x|≤1} g(x)dx. Suppose that g1 (x) = 1{x>1} g(x)/ν((1, ∞))is bounded, and there exists γ > 0 such that
lim eγx f0 (x) = 0.

(0.3)

x→∞

For a density f of μ we consider the following properties.
(i)f ∈ S+(ii)

g1 ∈ S +(iii)

g1 ∈ Land&f is al.d.
lim f (x)/g1 (x) = ν((1, ∞)).

x→∞(a) If g is a.n.i., then we can choose f such that (i), (ii) and (iii) are equivalent.
(b) If g is al.d., then we can choose f such that (ii) ⇔ (iii) implies (i).
We could remove several conditions in Theorem 0.3 by assuming the absolute integrability of thespectrally positive part μ+ (z).
0.4. Let μ ∈ ID(R) with and ν(dx) = g(x)dx such that g1 (x) is bounded. Suppose that
Theorem∞∞|μ(z)|dz< ∞, which implies −∞ |μ(z)|dz < ∞, so that μ has a bounded continuous density+−∞
f . Then the following relations hold between the properties (i), (ii) and (iii) of Theorem 0.3.
(a) If g is a.n.i., then we can choose f such that (i), (ii) and (iii) are equivalent.
(b) If g is al.d., then we can choose f such that (ii) ⇔ (iii) implies (i).
We apply our results to the consistency proof of the maximum likelihood estimation (MLE forshort) for μ ∈ ID(R) which is absolutely continuous. For simplicity we put a = b = 0 in μ(z) of(0.2) and assume that the spectrally positive part μ+ (z) is absolutely integrable.
Let f (x; θ) be the density of μ with θ a parameter vector and g(x; θ) be a density of the corresponding Lévy measure ν. Let (X1 , . . . , Xn ) be a random sample from f (x; θ0 ) with θ0 ∈ Θ whereΘ is a compact parameter space. Deﬁne the likelihood functionnMn (θ) = n−1
log f (Xi ; θ).
i=1

MLE θ̂n maximizes the function θ → Mn (θ). We say that a function α(x; θ) is identiﬁable ifa.e.
α(· ; θ) = α(· ; θ  ) every θ = θ ∈ Θ, i.e. α(x; θ) = α(x; θ ) does not hold. For convenience, weonly consider the symmetric or positive-half case, but we can easily generalize the result in the
non-symmetric two-sided case. We use the function g1 deﬁned in Theorem 0.4.
Proposition 0.5. Let μ ∈ ID(R) given by (0.2) with a = b = 0 such that μ+ (z) is absolutely integrable. Let g(x; θ) be a symmetric or positive-half density of ν. Suppose (i) : g(x; θ) is identiﬁable,θ → g(x; θ) is continuous in θ for every x, and (supθ∈Θ | log g1 (x; θ)|)g1 (x; θ0 )dx < ∞ with Θ acompact set such that θ0 ∈ Θ. Suppose (ii) : g1 (x; θ) is bounded and a.n.i., and g1 ∈ S. Then MLEpθ̂n satisﬁes θ̂n → θ0 .
References[1] Matsui, M. (2022) Subexponentialiy of densities of inﬁnitely divisible distributions (submitted).

- 35 -

ʮΨϯϚ෼෍ͷ‫ܗ‬ঢ়ύϥϝʔλͷϕΠζਪఆʯ౦‫ژ‬େֶ‫ֶࡁܦ‬෦ɹɹɹೖߐɹ‫܆‬༧ߘͷ௨ΓɺΨϯϚ෼෍ͷ‫ܗ‬ঢ়ύϥϝʔλͷϕΠζਪଌʹؔ͢Δ‫ߨ͍ͯͭʹڀݚ‬ԋͨ͠ɻߨԋ‫ʹޙ‬͸ɺҎԼͷ࣭͓ٙΑͼٞ࿦͕͋ͬͨɻ• ఏҊख๏ͷΞΠσΟΞͷ֬ೝ͕͋ͬͨɻPolya-inverse gamma ෼෍ͱ͍͏ෳࡶͳ෼෍Λ༻͍Δઌߦ‫ڀݚ‬ͱ͸ҟͳΓɺఏҊख๏Ͱ͸ϕʔλ෼෍ʹΑΔ؆қͳࠞ߹ද‫ݱ‬Λ‫͍͚ͯͭݟ‬Δ͜ͱɺͦͷਖ਼౰ԽΛ MH ๏ʹΑͬͯߦ͍ͬͯΔ͜ͱɺͦͯ͠ MH ๏ͷड༰֬཰͕ߴ͍͜ͱΛ਺ֶతʹ֬ೝͰ͖Δ͜ͱΛ੔ཧͨ͠ɻ

• ओఆཧͷূ໌ʹ͍ͭͯิ଍ͨ͠ɻΨϯϚؔ਺ͷ‫਺ٯ‬ͷྦྷ৐ʹରͯ͠ɺ෼ࢠɾ෼฼ʹҾ਺ΛͣΒͨ͠ΨϯϚؔ਺Λ͔͚Δ͜ͱͰϕʔλؔ਺Λ࡞Γͩ͢͜ͱɺ࢒ͬͨΨϯϚؔ਺ͷੵʹ Gauss’ multiplication

formula Λద༻͢Δ͜ͱͰ୯ҰͷΨϯϚؔ਺ʹॻ͖௚ͤΔ͜ͱΛࢦఠͨ͠ɻ• ఏҊख๏͕ଞͷύϥϝτϦοΫͳ෼෍ʹର͢Δ MCMC ๏ʹԠ༻ՄೳͰ͋Δ͜ͱʹؔ࿈ͯ͠ɺҰൠԽ͞Εͨଟ߲෼෍΁ͷద༻Մೳੑ͕ࣔࠦ͞Εͨɻ• ࣄલ෼෍ͱͯ͠཭ࢄ෼෍Λ༻͍ΔΞϓϩʔνͱͷൺֱΛ‫ٻ‬ΊΒΕɺಛʹෆ࣮֬ੑͷධՁʹ͓͍ͯ཭ࢄ෼
෍ʹΑΔΞϓϩʔν͸ਖ਼֬ੑΛ͓ܽͦ͘Ε͕͋Δ఺Λࢦఠͨ͠ɻ• ఏҊख๏ͷଟ߲σΟϦΫϨ෼෍΁ͷԠ༻ʹؔ࿈ͯ͠ɺߴ࣍‫ݩ‬ɾେ‫ن‬໛σʔλͷ৔߹ͷ‫ࢉܭ‬Մೳੑʹ͍࣭ͭͯ໰͕͋ͬͨɻཧ࿦্͸જࡏม਺ͷαϯϓϦϯάͷฒྻԽ͕ՄೳͰ͸͋Δ͕ɺ࣮ࡍʹ‫Ͳ͕ؒ࣌ࢉܭ‬Ε͚͔͔ͩΔ͔ʹ͍ͭͯ͸‫ڀݚ‬՝୊ͱͨ͠ɻ

• ΨϯϚ฼ूஂ෼෍ͱ͍͏୯७ͳઃఆͰ͋ΔͷͰɺ‫ܗ‬ঢ়ύϥϝʔλͷਪଌʹ͓͍ͯ΋ɺϕΠζਪఆྔ͓Αͼ MSE ͸ຊ౰ʹղੳతʹ‫͍ͳ͖Ͱࢉܭ‬ͷ͔ͱ͍͏࣭໰ʹରͯ͠ɺ‫ط‬஌ͷಛघؔ਺Ͱͷ‫ه‬ड़͸೉͍͠ͱฦ౴ͨ͠ɻ• ઌߦ‫ڀݚ‬ͷ Polya-inverse gamma ෼෍͕ແ‫ݶ‬෼ղՄೳͰ͋Δ͜ͱ͕ࢦఠ͞Εɺؔ࿈ͯ͠ Polya-gamma෼෍ͳͲͷઌߦ‫ڀݚ‬Λ঺հͨ͠ɻ

- 36 -

カーネル型推定量を利用した推測について
中央大学理工学部

前園宜彦

1. カーネル型推定量
X1 , X2 , ..., Xn を互いに独立で同じ分布にしたがう確率変数とし，分布関数を FX ，密度関数を fX とする．この分布関数 FX の推定量としては経験分布関数
Fn (x) =n
1
I(Xi ≤ x),n i=1

x ∈ R,がある．ここで I(A) は定義関数である．経験分布関数は階段関数となり連続ではない．
連続になるような推定量でよく利用されるのがカーネル型推定量である．
カーネル型分布関数推定量（Nadaraya (1964)）はFX (x) =

n
1
K
n i=1



x − Xi,hx ∈ R,

∞で与えられる．ただし k(·) はカーネル関数とし， −∞k(v)dv = 1 の条件を満たす．バンv
ド幅は n → ∞ のとき h → 0 かつ nh → ∞ とする．このとき K(v) = −∞k(u)du である．適当な条件の下で n → ∞ のとき

Bias[FX (x)] =V ar[FX (x)] =



∞
h2 
fX (x)
v 2 k(v)dv + o(h2 ),
2
−∞

h
1
2h
FX (x)[1 − FX (x)] −r1 fX (x) + o
n
n
n



∞
となる．ただし r1 = −∞vk(v)K(v)dv で，多くの場合非負の値をとる．
2. 確率点推定
p-確率点は Q(p) = F −1 (p) で定義される．経験分布関数に基づく p-確率点 Q(p) の推定
量はQ(p) = Fn−1 (p) = inf{x; Fn (x) ≥ p}で与えられる．漸近平均二乗誤差は

AM SE(Q(p)) ={Q (p)}2 p(1 − p)p(1 − p)=nf 2 (Q(p))
nとなる．Falk(1984) や Maesono & Penev(2011) では，カーネル関数 k(·) を利用した滑らかな推定量
x−p
1 1 −1
Qp,h =Fn (x)Kdxh 0hの漸近分布について研究している．漸近平均二乗誤差は


AM SE(Qp,h ) ={Q (p)}2 p(1 − p) 2h  2
1 1+
Q (p) − +nn2 2で与えられる．

- 37 -

 1−1

K 2 (x)dx

3. 順位検定の連続化への応用
X1 , X2 , . . . , Xn を互いに独立で同じ母集団分布 F (x − θ) にしたがう無作為標本とする．
ここで対応する確率密度関数は f (−x) = f (x) を満たす原点対称な分布とする．θ は未知
母数で，帰無仮説 H0 : θ = 0 に対して対立仮説 H1 : θ > 0 の検定問題を考える．この問題に対する順位検定統計量の有意確率については Lehmann & D’abrera(2006) でも指摘されているように，検定統計量の取り得る値が細かいほど有意確率が小さくなる傾向がある．
ψ(x) = 1 (x ≥ 0), = 0 (x < 0) とおくとき，符号検定 S とウィルコクソンの符号付き順位検定 W はS=nψ(Xi ),

W =

i=1



ψ(Xi + Xj )1≤i≤j≤nである．これら S, W の連続化として
n




Xi
S = n − nFn (0) = n −
K −
,
h
i=1

! = n(n + 1) −W2

1≤i≤j≤nK −X i + Xj2hが考えられる．この統計量の分布は帰無仮説の下でも母集団分布に依存することになるが，! は S, W を滑らかにしたものに漸近平均と漸近分散は F (·) に依存しない．この S, W
なっている．
4. 境界バイアス縮小す推定量に基づくノンパラメトリック検定
F を仮定された分布関数とするとき，次の検定問題を考える．
帰無仮説 H0 : FX = F

対立仮設 H1 : FX = F

このような一般的な設定の下でよく利用される検定として

KSn = sup |Fn (x) − F (x)|,CvMn = n

x∈R

 ∞

−∞

[Fn (x) − F (x)]2 dF (x)を検定統計量とするコルモゴロフ・スミルノフとクラーメル・フォンミーゼス検定統計量が知られている．
Rizky & Maesono(2022) は全単射関数 g を用いたバイアスを縮小するカーネル型推定量n1
FX (x) =
W
n i=1



g −1 (x) − g −1 (Xi )hを用いて" = sup |FX (x) − F (x)|,
KS
x∈R

" =n
CvM

 ∞
−∞



, x∈Ω

[FX (x) − F (x)]2 dF (x)を提案している．これに対して" →p 0,|KSn − KS|" | →p 0|CvMn − CvM" と CvM" を利用した有意確率は，通常のコが成り立つ．このことより，検定統計量 KS
ルモゴロフ・スミルノフ検定統計量およびクラーメル・フォンミーゼス検定統計量に対する数表を利用して検定できる．シミュレーションの結果では，検出力が高いことが示された．

- 38 -

On a generalization of Clayton-Oakes model by R. L. Prentice
Hideatsu Tsukahara∗
January 3, 2023
When we have two failure times T1 and T2 and our primary interest is in the association
between them, we need families of bivariate survivor functions S(t1 , t2 ) for statistical modeling.
Postulating some desirable properties on certain conditional hazard functions, Clayton [1] derived such a family which is now called Clayton-Oakes model. In terms of the associated copula,
it is written asC(u1 , u2 ) = (u−θ + v −θ − 1)−1/θ ∨ 0,

θ ∈ [−1, ∞) \ {0}.

This family is an example of Archimedean copulas, which are closely related to the frailty modelin survival analysis as shown in Oakes [3]. The cross ratio function of T1 and T2 is deﬁned byθ∗ (t1 , t2 ) :=S(t1 , t2 )∂12 S(t1 , t2 )∂1 S(t1 , t2 )∂2 S(t1 , t2 )Oakes [3] shows that the cross ratio function is a function of survivor function only if and only
if the associated copula is Archimedean. The Clayton-Oakes model amounts to assuming thatthe cross ratio function is constant.
Prentice [4] proposed a generalization of the Clayton-Oakes model; he wrote “A trivariate
generalization with unrestricted marginal and pairwise marginal survivor functions isS(t1 , t2 , t3 ) = S(t1 , t2 , 0)−θ + S(t1 , 0, t3 )−θ + S(0, t2 , t3 )−θ− S(t1 , 0, 0)−θ − S(0, t2 , 0)−θ − S(0, 0, t3 )−θ + 1

−1/θ

∨ 0,

where −1 ≤ θ < ∞.” However, this is a functional equation in S(t1 , t2 , t3 ), and it is not clearthat there exists a 3-dimensional survivor function S(t1 , t2 , t3 ) satisfying this equation. So the
problem should be correctly posed in the following way.
Problem Given three bivariate sf’s S12 (t1 , t2 ), S13 (t1 , t3 ), S23 (t2 , t3 ) satisfying the compatiblility conditions
S12 (t1 , 0) = S13 (t1 , 0) =: S1 (t1 ),S12 (0, t2 ) = S23 (t2 , 0) =: S2 (t2 ),S13 (0, t3 ) = S23 (0, t3 ) =: S3 (t3 ),is the function G deﬁned byG(t1 , t2 , t3 ) = S12 (t1 , t2 )−θ + S13 (t1 , t3 )−θ + S23 (t2 , t3 )−θ− S1 (t1 )−θ − S2 (t2 )−θ − S3 (t3 )−θ + 1
∗

−1/θ

∨0

Faculty of Economics, Seijo University, 6–1–20 Seijo, Setagaya-ku, Tokyo, 157-8511, Japan, E-mail:tsukahar@seijo.ac.jp- 39 -

a proper survivor function?
The problem can be stated in terms of copulas as well. Let C12 (u1 , u2 ), C13 (u1 , u3 ), C23 (u2 , u3 )be the copula associated with S12 (t1 , t2 ), S13 (t1 , t3 ), S23 (t2 , t3 ); namely
S12 (t1 , t2 ) = C12 (S1 (t1 ), S2 (t2 )),S13 (t1 , t3 ) = C13 (S1 (t1 ), S3 (t3 )),S23 (t2 , t3 ) = C23 (S2 (t2 ), S3 (t3 )).
The problem is then reduced to whether the functionG(u1 , u2 , u3 ) = C12 (u1 , u2 )−θ + C13 (u1 , u3 )−θ + C23 (u2 , u3 )−θ
−θ
−θ
− u−θ
1 − u2 − u3 + 1

−1/θ

∨0

is a (proper) copula for any given copulas C12 , C13 and C23 .
One can easily observe that the answer is “no”; if the answer were yes, then the Fréchet class
F(F12 , F13 , F23 ) (the set of trivariate distributions with three given bivariate margins) would
always be non-empty (see Joe [2, Section 3.4]).
In this expository note, we begin with the review of the original Clayton-Oakes model andthe key concept of cross ratio, give an explicit counterexample to the above problem, and discuss
under which conditions Prentice’s generalization gives a valid survivor function.

References[1] D. G. Clayton. A model for association in bivariate life tables and its application in epidemiological studies of familial tendency in chronic disease incidence. Biometrika, Vol. 65,
pp. 141–151, 1978.
[2] H. Joe. Multivariate Models and Dependence Concepts. Chapman and Hall, London, 1997.
[3] D. Oakes. Bivariate survival models induced by frailties. Journal of the American Statistical
Association, Vol. 84, pp. 487–493, 1989.
[4] R. L. Prentice. Higher dimensional clayton–oakes models for multivariate failure time data.
Biometrika, Vol. 103, pp. 231–236, 2016.

- 40 -

Kendall の順位相関係数を固定した下での最小情報コピュラ
助田 一晟 ∗清智也 †1 最小情報コピュラ
一般に d 次元コピュラとは全ての周辺分布が [0, 1] 上の一様分布である d 次元の分布関数を指し, 任意の d 次元分布は Sklar の定理によりコピュラと周辺分布を用いて表現することができることから, コピュラは 2 変数間の従属性を完全に記述しているとされる. 所与の制約条件下で, 独立コピュラ ([0, 1]2 上の一様密度を持つコピュラ) に

Kullback-Leibler 情報量の意味で最も近いコピュラを最小情報コピュラと呼ぶ. このコピュラは以下の最適化問題の解として与えられるものである（例えば Bedford et al. [1]）.

 1 1minimize
 1
s.t.
0

0

0

p(x, y) log p(x, y)dxdy

p(x, y)dy = 1,

 1 100 10p(x, y)dx = 1

hk (x, y)p(x, y)dxdy = μk

(1.1)
(1.2)
(1.3)

最小情報コピュラの離散近似版として最小情報チェス盤コピュラが考えられている [3]. 最小情報コピュラと同様に, エントロピー最大化原理に則った枠組みで定式化を行うことで, 最小情報チェス盤コピュラは以下の最適化問題の解と捉えることができる.

minimizen nπij log πij

(1.4)

1
1 
,
πij =
n i=1
n

(1.5)

i=1 j=1

s.t.

n


n

πij =

j=1
n 
n


hk,ij πij = μk

(1.6)

i=1 j=1

ただし, h1 , . . . , hK は所与の関数, hk,ij は領域 Dij の中心での hk (x, y) の値, θ = (θ1 , . . . , θK ) は未知パラメータである. この最適化問題は有限次元凸最適化問題となっており, 最適解が一意に存在する.

2 ２次の最小情報コピュラへの拡張
従来の最小情報コピュラは所与の制約として 1 次の母数（ある関数の期待値の形で表されるもの）のみを扱っていた. たとえば, Spearman の順位相関係数は 1 次の母数であり, Spearman の順位相関係数を固定した状況は従来の最小情報コピュラの枠組みで扱うことができる. 一方, Spearman の順位相関係数と比較されることが多い Kendall の
順位相関係数は 2 次の母数であるため, Kendall の順位相関係数を固定した状況は従来の最小情報コピュラの枠組みで扱うことができない. そこで本研究では所与の制約式を 2 次の母数である Kendall の順位相関係数としたときの最適解（最小情報チェス盤コピュラ）を考察する.
∗ 東京大学大学院情報理工学系研究科数理情報学専攻† 東京大学大学院情報理工学系研究科数理情報学専攻- 41 -

従来の最小情報チェス盤コピュラと同様に, Kendall の順位相関係数を τ (∈ [−1, 1]) に固定した下での最小情報
チェス盤コピュラは以下のように最適化問題の解として定式化される. ただし最後の制約はチェス盤コピュラ上のKendall の順位相関係数 [2] を固定している.
(M P ) minimize

n 
n


πij log πij(2.1)

i=1 j=1

s.t.

n


πij =

j=1

n
1 
1
,
πij =
n i=1
n

⎛

1
⎜ 2
T
1 − Tr(ΞΠΞΠ ) = τ, Π = (πij ), Ξ = ⎜⎝. . .
2

(2.2)
0
1

...
...
2

⎞
0
0⎟
⎟
⎠
1

(2.3)

まず, 各領域の確率質量（分割表モデルにおける各マス目の値）を移動させる操作を以下のように導入する. チェ
ス盤コピュラ全体の空間を一様コピュラ U = n12 I を原点とするベクトル空間 (の部分空間) とみなし, (n − 1)2 本から成る斜交基底 (Tij )i=1,...,n−1, j=1,...,n−1 を考える. ただし, Tij は i 行 j 列と i + 1 行 j + 1 列が +1, i 行 j + 1 列と i + 1 行 j 列が −1, それ以外が 0 となっている行列である. 任意のチェス盤コピュラ Π は以下のように一意な係数

aij を用いて表記できる : Π = U +

i,j aij Tij , (i = 1, . . . , n − 1, j = 1, . . . , n − 1)

この操作に対する最適化問題 (M P ) の停留条件を考えることで, 以下の主張を得た.

Lemma 1. 任意のチェス盤コピュラ Π と任意の成分移動操作 Tij と Ti j  に対し,
Π := Π + Tij − r Ti j  , r :=と定める. このとき,

πi,j + πi+1,j + πi+1,j + πi+1,j+1π+ πi +1,j  + πi +1,j  + πi +1,j  +1i ,j T

1 − Tr(ΞΠΞΠT ) = 1 − Tr(ΞΠ ΞΠ ) + O( 2 )

(2.4)

(2.5)

Lemma 2. チェス盤コピュラ Π = (πij ) に対し, 任意の成分移動操作 Tij を微小に行なった時の情報量の変分はlog

πi,j πi+1,j+1
πi+1,j πi,j+1

Theorem 1 (2 次の最小情報チェス盤コピュラの不変量). 2 次の最小情報チェス盤コピュラでは任意の (i, j) の組
み (i, j = 1, . . . , n − 1) に対し, 以下の値が一定.

πi,j πi+1,j+1
1log
πi,j + πi+1,j + πi+1,j + πi+1,j+1
πi+1,j πi,j+1
また, 発表ではその他 Total Positivity や幾何構造等の諸性質や, 提案した最小情報コピュラモデルを実際の株価データにモーメント法を基にフィッティングさせた結果を紹介した.

参考文献[1] Bedford, T. and Wilson, K. J., On the construction of minimum information bivariate copula families, Ann.
Inst. Stat. Math., 66, 703–723, 2014.
[2] Durrleman, V., Nikeghbali, A. and Roncalli, T., Copulas Approximation and New Families, Available at
SSRN, 2000.
[3] Piantadosi, J., Howlett, P., and Boland, J., Copulas with maximum entropy, Optim. Lett., 6, 99–125, 2012.

- 42 -

スピアマンランク行列による主成分分析のロバストネス
千葉大・融合理工学府

渡邊 宏大

千葉大・理学研究院

内藤 貫太

はじめに: ロバストな主成分分析については，これまで多くの手法が提案されている．本発表では，Marden(1999) の Rank から派生したスピアマンランク行列に基づく主成分分析のロバスト性について報告した．Han
and Liu (2018) で議論されている Kendall’s tau との比較についても考察した．
設定と定義:y ∈ Rd の Spatial sign をS(y) =

y
y2
0

,

y = 0,y=0で定義する．X, Y , Z, X1 , . . . , Xn を有限な分散共分散行列 Σ を持つ d 次元確率分布 F にしたがう互いに独立な確率ベクトルとする．Spatial sign に基づくロバストな主成分分析の手法が，Marden (1999)，Han andLiu (2018) で議論されている．
母集団での Kendall’s tau matrix K ∈ Rd×d は
K = EX,Y S(X − Y )S(X − Y )Tで定義され，その推定量は

)=
K


2
S(Xi − Xj )S(Xi − Xj )Tn(n − 1) j<iで定義される．一方，母集団でのスピアマンランク行列 ΣR ∈ Rd×d は

*
+


T
ΣR = EX,Y ,Z S(X − Y )S(X − Z)T = EX EY [S(X − Y )] EZ [S(X − Z)]で定義され，その推定量は
n

1
ΣR =
n(n − 1)(n − 2) i=1

S(Xi − Xj )S(Xi − Xk )Tj=i,k=i,j=kで定義される．
問題: 重要な事実として，F が楕円分布のとき，Σ，K および ΣR は同じ直交行列で対角化されることが知られている (Marden (1999)，Han and Liu (2018) を参照)．Marden (1999) においても，ΣR の推定量 ΣR によるロバストな主成分分析が議論されているが，その ΣR は

ΣR =


1)
2
K + 1−nn- 43 -

ΣR) は K に，ΣR は ΣR に確率収束することが示される．このことから，ΣRと分解される．n → ∞ のとき，K) の比較は，漸近的には ΣR と K) の比較となる．我々の問題は特に，ΣR と K) はどちらがロバストな主成とK分分析を提供するのか？ということになる．
影響関数: ロバストネスの指標として影響関数を考える (影響関数については，Hampel et al. (1986) を参照)．

F を楕円分布とする．比較の対象として，K と ΣR の最大固有値に関する影響関数および固有ベクトルに関する影響関数を与える．

K は直交行列 Γ = [γ1 · · · γd ] と対角行列 A = diag(a1 , . . . , ad )，a1 > · · · > ad > 0 によって K = ΓAΓTとスペクトル分解されるとする．上で述べた重要な事実から，ΣR も同じ直交行列によって ΓΛΓT ，Λ =diag(λ1 , . . . , λd ) とスペクトル分解できる．ただし，対応する固有値については，λ1 > · · · > λd > 0 となるかは定かではない．したがって，γ1 が属している固有値を λM ，M ∈ {1, . . . , d}，ΣR の最大固有値と対応する固有ベクトルをそれぞれ λmax ，γL ，L ∈ {1, . . . , d} とする．このとき，K と ΣR の最大固有値に関する影響関
数 IFK (x, a1 , F )，IFΣR (x, λmax , F ) はそれぞれ

IFK (x, a1 , F ) = γ1T A(x)γ1 − 2a1 ,(1)IFΣR (x, λmax , F ) = γLT B(x)γL − 3λmax(2)で与えられる．ここで，
A(x) = 2EX S(x − X)S(x − X)T ,
B(x) = EX,Y S(X − Y )S(X − x)T + S(X − x)S(X − Y )T + S(x − X)S(x − Y )Tである．また，K と ΣR の固有ベクトル γ1 に関する影響関数 IFK (x, γ1 , F )，IFΣR (x, γ1 , F ) はそれぞれd1γ T A(x)γ1 · γk ,a1 − ak k
k=2
1
γ T B(x)γ1 · γk
IFΣR (x, γ1 , F ) =λM − λk k
IFK (x, γ1 , F ) =

(3)
(4)

k=M

で与えられる．
ロバストネスの比較と適用例： 本講演においては，いくつかの 2 次元楕円分布のもとで，固有値，固有ベクトルの影響関数 (1)，(2)，(3)，(4) の曲面を描画することおよび固有空間の類似度の比較を通して，ΣR による主成分分析のロバスト性を報告した．外れ値を含む実データへの適用結果についても報告した．

参考文献
Hampel, F. R., Ronchetti, E. M., Rousseeuw, P. J., and Stahel, W. A. (1986). Robust Statistics : The
Approach Based on Inﬂuence Functions. Wiley.
Han, F. and Liu, H. (2018). Eca: High-dimensional elliptical component analysis in non-gaussian distributions. Journal of the American Statistical Association, 113:252–268.
Marden, J. I. (1999). Some robust estimates of principal components. Statistics and Probability Letters,
43:349–359.

- 44 -