科研費シンポジウム「多様な分野における統計科学に関する諸問題」（Problems related to statistical science in various ﬁelds）
日時：2019 年 9 月 14 日（土）〜9 月 16 日（月）(Date: Saturday, 14, – Monday, 16, September, 2019)
場所：コープシティ花園 4F ガレッソホール a(Place: Coop City Hanazono 4F Garesso Hall a)（TEL: 025-248-7511）
科学研究費・基盤研究（A) （課題番号：15H01678）「大規模複雑データの理論と方法論の総合的研究」
研究代表者：青嶋 誠（筑波大学）
Makoto Aoshima (University of Tsukuba)
開催責任者：蛭川 潤一（新潟大学）
Junichi Hirukawa (Niigata Universiry)Program
Saturday, 14, September
Reception 13:00-13:25
Opening 13:25-13:30 Junichi HIRUKAWA (Niigata University)
Afternoon Session I (in Japanese) 13:30-15:00
Chair Kou Fujimori (Waseda University)
1. 13:30-14:15 黒木 裕鷹 (KUROKI, Yutaka)
東京理科大学大学院工学研究科 (Graduate School of Engineering, Tokyo University of Science)「レーティング手法を用いたネットワーク中心性の構築とその統計的性質」(Constructing Network Centrality Measure based on Rating Methods and its Statistical Properties)
2. 14:15-15:00 KOTSUBO, Takuto
Graduate School of Engineering, Tokyo University of Science
Hidden Markov models for cylindrical data and its application for animal movement analysis
Coﬀee Break 15:00-15:15
Afternoon Session II (in Japanese) 15:15-16:45
Chair Takayuki Shiohama (Tokyo University of Science)
3. 15:15-16:00 Nakayama Yugo
Graduate School of Pure and Applied Sciences, University of Tsukuba
Support vector machine and optimal choice of parameters for high-dimensional imbalanced data
4. 16:00-16:45 新村秀一 (Shinmura Shuichi)
成蹊大学名誉教授 (Emeritus Professor of Seikei Univ.)「高次元遺伝子データ解析の呪いからの解放」(Release of Curse of High-dimensional Data Analysis)

Sunday, 15, September
Morning Session (in English) 9:15-11:30
Chair Muneya Matsui (Nanzan University)
5. 9:15-10:00 Yuichi Goto, Marc Hallin, Masanobu Taniguchi
Waseda University, Université libre de Bruxelles, Waseda University
Kolmogorov-Smirnov Tests for Laplace Spectral Density Kernels
6. 10:00-10:45 Kou Fujimori, Sota Sakamoto and Yasutaka Shimizu
Waseda University
Generalized maximum composite likelihood estimator for determinantal point processes
7. 10:45-11:30 Junichi Hirukawa and Kou Fujimori
Niigata University and Waseda University
Weak convergence of the partial sum of I(d) process to a fractional Brownian motion in ﬁnite
interval representation
Lunch 11:30-13:15
Afternoon Session I (in English) 13:15-15:30
Chair Nobuaki, HOSHINO (Kanazawa University)
8. 13:15-14:00 Toshihiro Abe
Nanzan University
A closed form EM algorithm for a multivariate skew-normal model
9. 14:00-14:45 Muneya Matsui
Nanzan University
Asymptotics of maximum likelihood estimation for stable law with continuous parameterization
10. 14:45-15:30 張 元宗 (Chang Yuan Tsung), 篠崎 信雄 (Nobuo Shinozaki), William, E.
Strawderman
目白大学・社会学部・社会情報学科（Mejiro University, Department of Social Information, Facultyof Studies on Contemporary Society),
慶應大学・理工学部 (Faculty of Science and Technology, Keio University),

Department of Statistics and Biostatistics, Rutgers University
Pitman Closeness Domination in Predictive Density Estimation for Two Ordered Normal Means
Under α-Divergence LossCoﬀee Break 15:30-15:45

Afternoon Session II (Guest speakers session) 15:45-16:45
Chair Junichi Hirukawa (Niigata University)
11. 15:45-16:45 Konstantinos Fokianos
Lancaster University
Auto-Distance Covariance Function for Time Series Analysis

Monday, 16, September
Morning Session (in Japanese) 9:15-11:30
Chair Toshihiro Abe (Nanzan University)
12. 9:15-10:00 入江 薫 (Kaoru Irie)
東京大学経済学部 (Faculty of Economics, The Univeristy of Tokyo)「縮小事前分布と状態空間モデル」(Shrinkage priors and state space models)
13. 10:00-10:45 永井 勇 (Isamu Nagai)
中京大学 国際教養学部（School of International Liberal Studies, Chukyo University）「バランス型経時測定データにおける Extended GMANOVA モデルの解釈と新たな推定法」(Interpretation for the extended GMANOVA model in the balance type longitudinal data, andnew estimation method)
14. 10:45-11:30 星野 伸明 (Nobuaki, HOSHINO)
金沢大学経済学類 (School of Economics, Kanazawa University)「分散可変な一般化多項分布について」(Generalized multinomial distributions with scalable variance)
Closing 11:30-11:35 Junichi HIRUKAWA (Niigata University)レーティング手法を用いたネットワーク中心性の構築とその統計的性質黒木裕鷹
東京理科大学大学院工学研究科経営工学専攻
塩濱敬之
東京理科大学工学部情報工学科

概要
ソーシャル・ネットワーキング・サービス (SNS) の普及と定着により, 人間関係におけるコメントや意見のやり取りを記録したネットワークデータが大量に蓄積されている. ネットワー
ク分析において, 中心的な役割を果たす個人 (ノード) を特定することが重要な課題であり, これまで様々な中心性指標が提案されてきた. 代表的な中心性指標のほとんどはネットワークの構造に由来するが, 正規方程式に基づいたレーティングアルゴリズムである Massey のレーティン
グは, エッジの重みを考慮した中心性として解釈することができる. 本研究では, Massey のレー
ティングを中心性指標として再解釈し, その統計的性質を検討する.

1 はじめに近年, 価値や金銭のやりとり, 交友関係に着目した社会的ネットワークの分析が, グラフ理論や統計的手法を用いて盛んに研究されている. 特に, 情報技術の発展により現実世界に存在する巨大で複雑なネットワークの性質を明らかにする複雑ネットワーク分析が注目されている. 複雑ネット
ワークの分析を通じて, ネットワークを構成するコミュニティの理解や位置付け, 将来時点のネット
ワーク特性の予測が可能となることにより, ソーシャルネットワーキングサービス (SNS) やブロッ
クチェーン技術を利用した情報サービスの新たな価値創造が期待されている. 社会ネットワークを構成する個人や組織をノード (頂点), ノード間の関係性をエッジ (辺) といい, ノードとエッジで構成されるグラフはネットワークの社会的相互作用を記述する. 相互関係のあるデータにおいて重要な課題は, ネットワーク内で重要な役割を果たすノードを把握することにある. このような重要性の指標は中心性 (centrality) と呼ばれ, 主に社会科学の分野で様々な中心性指標が提案されてきた.
たとえば, Bonacich (1987) や Katz (1953) を参考にすること.
代表的な中心性指標には, 次数中心性, 媒介中心性, 近接中心性, 固有ベクトル中心性などが知られている. これらの中心指標はネットワークの構造から直接求めることができるがノードの特性やエッジの重みは考慮されない. たとえば, SNS における「賛成」「反対」の絶対数や, 投げ銭システ
ムに代表されるようなトークンの付与回数などを考慮した分析において, このような中心性指標を用いてノードの重要性を評価することはできない. 関係性のセマンティクスを考慮した, 非対称な有向グラフを用いたネットワーク指標を構成するには, レーティングシステムやランキングアルゴ
リズムを利用した中心性の指標の構築が必要となる. このような研究例として, スポーツやゲームの試合結果を利用するレーティングシステムや Page Rank に代表されるランキングシステムがあり, 検索エンジンの設計や推薦システム, 影響力のある拡散者の特定など, 広範囲に応用されている.

1

スポーツやゲームのデータ解析では, 「強さ」の指標であるレーティングの開発が盛んに行われており, 算出されたレートは結果の予測やトーナメントの出場プレーヤーの選出など, 様々な場面で利用される. また, Page Rank のアルゴリズム (Brin and Page, 1998) は, Google 検索エンジンにおいて Web ページをランク付けするために考案された. スポーツやゲームにおけるレーティング
は, 一対比較モデルを応用して勝敗確率を予測する Bradley-Terry モデル (Bradley and Terry, 1952),
チェスに端を発する Elo レーティング (Elo, 1978), ラプラスの継起則を応用した Colley のレーティ
ング (Colley, 2002) などが存在する. 全米大学体育協会 (NCAA) の主催するバスケットボールトー
ナメントの出場チームの選出には, アルゴリズムは非公開ではあるものの, シーズン中のデータに基づいた指標である Rating Percentage Index (RPI) が利用されている. 動的なレーティングの更新が可能な Elo レーティングはもともとチェスの実力を表すために考案されたものだが, 世界ラグビー協会は修正した手法を採用している. また, 正規方程式に基づいたレーティング手法として, Masseyのレーティング手法 (Massey, 1997) が知られている. Massey の手法は, 代表的なネットワーク中心性である Katz の中心性 (Katz, 1953) と関連する.
本研究では, Massey のレーティング手法をネットワーク中心性として再解釈し, その統計的性質を検証する. 本稿の構成は以下の通りである. 第 2 章では, 代表的な中心性指標をいくつか紹介し,各指標の特徴を述べる. 第 3 章では Massey のレーティング手法について説明し, 中心性指標としての解釈, 特に Katz の中心性との関連を考察する. また, レーティング手法の試合結果に対する感度に関する先行研究も取り上げる. 第 4 章ではデータ解析の例として, 2018 年テニス ATP ワールドツ
アーの試合結果を解析し, その統計的性質を明らかにする. 最後に, 今後の課題や展望について第 5章で述べる.

参考文献
Bonacich, P. (1987). Power and centrality: A family of measures. American journal of sociology,
92(5):1170–1182.
Bradley, R. A. and Terry, M. E. (1952). Rank analysis of incomplete block designs: I. the method of
paired comparisons. Biometrika, 39(3/4):324–345.
Brin, S. and Page, L. (1998). The anatomy of a large-scale hypertextual web search engine. Computer
networks and ISDN systems, 30(1-7):107–117.
Colley, W. N. (2002). Colley s bias free college football ranking method: The colley matrix explained.
Princeton University, Princeton.
Elo, A. E. (1978). The rating of chessplayers, past and present. Arco Pub.
Katz, L. (1953). A new status index derived from sociometric analysis. Psychometrika, 18(1):39–43.
Massey, K. (1997). Statistical models applied to the rating of sports teams. Bluefield College.

Hidden Markov models for cylindrical data andits application for animal movement analysis
Takuto Kotsubo
Graduate School of Engineering, Tokyo University of Science
Takayuki Shiohama
Department of Information and Computer Technology

Abstract
Analyzing animal tracking is a growing field in ecology and various models have been proposedin literature. State space model is often used to fit the animal tracking data, however the limitationsof the linear and Gaussian assumption are sometimes reported. An alternative approach for modeling animal movement data is Hidden Markov Models (HMM) with cylindrical distribution, where thecylindrical data consists of circular and linear variables. For linear part of the cylindrical distribution,
we adopt the generalized Palate-type distribution for heavy tailed observations in linear part, which
includes the Weibull-von Mises distributions on cylinder. Estimation for the model parameters aredone by using Expectation and Maximization (EM) algorithms, and we propose a modified M-stepprocedure and a Maximum at posteriori (MAP) estimation for model parameters to avoid local maximum solutions in EM algorithms.
Keywords: Animal movement, cylindrical data, EM-algorithms, maximum likelihood estimation.

1

Introduction

A cylindrical data refers to a bivariate data which consists of circular and linear variables. Circular
data usually deal with direction and has periodic properties with frequency 2π. Linear variables arise
in cylindrical data are usually defined on a positive real line. Cylindrical data arise often in various
fields in natural sciences, for example, wind directions and speeds in meteorology (Breckling (1989);
Ailliot et al. (2015)), wave directions and heights in oceanology (Ris et al. (1999)), ozone concentrations
and wind directions in environmental science (Camalier et al. (2007); Yi and Prybutok (1996)) and theturning angels and step lengths of the animal movement in ecology (Jonsen (2016); Adam et al. (2019)).
These data consist of time series nature, however time series modeling of cylindrical data are not fully
investigated in literature. The reason behind this is that there exists a few distributions on cylinder.
The well known cylindrical distribution is called the Johnson and Wehrly distribution (Johnson andWehrly (1978)). The drawbacks of the Johnson and Wehrly distribution are reported in Abe and Ley(2017), where they pointed out that the mode of the linear variables always defined at 0, which can not
be applied for many actual data analysis. Recently, Abe and Ley (2017) proposed cylindrical distribution
with combining sine-skewed von Mises and Weibull distribution to implement the skewness structuresin directional variables and to have the mode at some point on the domain of the probability distributionfunction. Imoto et al. (2019) extends the Abe-Ley distributions to have more heavy tail distributions on1

real part. There exists another class of the cylindrical distributions which includes the distributions ofMardia and Sutton (1978) and Kato and Shimizu (2008). The normalizing constants of these distributions
are somewhat complex and that is not expressed in an analytic form, hence in this study we focus on
applying the distributions of Abe and Ley (2017) and Imoto et al. (2019).
The rest of paper is organized as follows. Section 2 introduces our HMM models. Several useful distributional properties of cylindrical distributions are explained. Section 3 provides the maximum likelihood
estimation procedures which utilizing EM algorithms. Some technical improvements are also discussed
to avoid the local maximum solution in optimization together with to avoid unboundedness in likelihood
functions. In Section 4, some Monte Carlo simulations are performed to assess the performance of our
proposed procedures. Section 5 provides real data analysis of animal movement trajectories. Section 6
concludes our paper.

References
Abe, T. and C. Ley (2017). A tractable, parsimonious and flexible model for cylindrical data, with applications.
Econometrics and Statistics 4, 91–104.
Adam, T., C. A. Griffiths, V. Leos-Barajas, E. N. Meese, C. G. Lowe, P. G. Blackwell, D. Righton, and R. Langrock(2019). Joint modelling of multi-scale animal movement data using hierarchical hidden markov models. toappear in Methods in Ecology and Evolution.
Ailliot, P., J. Bessac, V. Monbet, and F. Pene (2015). Non-homogeneous hidden markov-switching models forwind time series. Journal of Statistical Planning and Inference 160, 75–88.
Breckling, J. (1989). The analysis of directional time series: applications to wind speed and direction, Volume 61.
Springer Science & Business Media.
Camalier, L., W. Cox, and P. Dolwick (2007). The effects of meteorology on ozone in urban areas and their use inassessing ozone trends. Atmospheric Environment 41(33), 7127–7137.
Imoto, T., K. Shimizu, and T. Abe (2019). A cylindrical distribution with heavy-tailed linear part. to appear inJapanese Journal of Statistics and Data Science.
Johnson, R. A. and T. E. Wehrly (1978). Some angular-linear distributions and related regression models. Journalof the American Statistical Association 73(363), 602–606.
Jonsen, I. (2016). Joint estimation over multiple individuals improves behavioural state inference from animal
movement data. Scientific Reports 6, 20625.
Kato, S. and K. Shimizu (2008). Dependent models for observations which include angular ones. Journal of
Statistical Planning and Inference 138(11), 3538–3549.
Mardia, K. and T. Sutton (1978). A model for cylindrical variables with applications. Journal of the Royal
Statistical Society: Series B (Methodological) 40(2), 229–233.
Ris, R., L. Holthuijsen, and N. Booij (1999). A third-generation wave model for coastal regions: 2. verification.
Journal of Geophysical Research: Oceans 104(C4), 7667–7681.
Yi, J. and V. R. Prybutok (1996). A neural network model forecasting for prediction of daily maximum ozone
concentration in an industrialized urban area. Environmental pollution 92(3), 349–357.

Support vector machine and optimal choice ofparameters for high-dimensional imbalanced data
Yugo Nakayama
Graduate School of Pure and Applied Sciences, University of Tsukuba

1

Introduction

In this talk, we considered the classification for high-dimensional data. Suppose we have independent andd-variate two populations, πi , i = 1, 2, having an unknown mean vector µi and unknown covariance matrix
Σi (≥ O). Let ∆ = kµ1 −µ2 k2 , where k·k denotes the Euclidean norm. We have independent and identically
distributed (i.i.d.) observations, xi1 , . . . , xini , from each πi . We assume ni ≥ 2, i = 1, 2 and n1 ≤ n2 . Let
x0 be an observation vector of an individual belonging to one of the two populations. We assume x0 andxij s are independent. Let n = n1 + n2 .
In recent years, the margin-based classification methods such as support vector machine (SVM) arestudied for high-dimensional data. It is well known that many classifiers, including SVM, overfit the datain the high-dimensional setting. Nakayama et al. (2019) claimed that SVM gives inconsistent properties for
imbalanced data:e(1) = 1 + oP (1) and e(2) = oP (1) as d → ∞;ore(1) = oP (1) and e(2) = 1 + oP (1) as d → ∞,

(1)
(2)

where e(i) denote the error rate of misclassifying an individual from πi into the other class for i = 1, 2. The
cause of inconsistency is due to the imbalance of sample sizes and a huge bias of high-dimensional space.
In overcome this difficulty, Nakayama (2019) proposed linear robust SVM (RSVM) to imbalanced data. In
this paper, we investigate the asymptotic properties of RSVM more deeply in the framework of n/d → 0 as
m → ∞, where m = min{d, n1 , n2 }.

2Asymptotic properties of robust SVM

In this section, we compare SVM and RSVM for high-dimensional data. First, we introduce a formulation of
SVM. Let (x1 , . . . , xn ) = (x11 , . . . , x1n1 , x21 , . . . , x2n2 ). #S denotes the number of elements in a set S. The
discriminant function of SVM is given by
y(x) =

Xj∈Ŝ

α̂j tj xTj x +


X
1 X
tj −
α̂j 0 tj 0 xTj xj 0 ,
nŜ
0

(3)

j ∈Ŝ

j∈Ŝ

where Ŝ = {j|α̂j 6= 0, j = 1, . . . , n} and nS = #S. α̂ = (α̂1 , . . . , α̂n ) is obtained by solving the following
maximization
!
n
n
n
X
1 XX
T
α̂ = argmax
αj −αj αj 0 tj tj 0 xj xj 0
2 j=1
α
j=1
k=1

subject to0 ≤ αj ≤ C, j = 1, . . . , n, and

n
X
j=1

1

αj tj = 0,(4)

where C > 0 is a regulariation parameter. Note that (3) becomes a hard-margin type when C → ∞. Thereexist some xj s satisfying that tj y(xj ) = 1 (i.e., α̂j 6= 0). Such xj s are called the support vector. SVM
classifies x0 into π1 if y(x0 ) < 0 and into π2 otherwise.
The regularization parameter C influences the properties of SVM. Let ∆∗ = ∆ + tr(Σ1 )/n1 + tr(Σ2 )/n2 .
For C, we consider the following conditions:2
lim sup≤ 1; or
(5)
n
∆
m→∞
1 ∗C
2
> 1.
(6)
lim infm→∞ n1 ∆∗ C
Let δ = tr(Σ1 )/n1 − tr(Σ2 )/n2 . Then we have the following result.
Theorem 1. Assume some regularity conditions andlim inf |δ|/∆ < 1.

(7)e(i) = 1 + oP (1) as m → ∞ for i = 1, 2.

(8)

m→∞

Then, under (5), (3) holds (8).

Under (6), we have the following result.
Theorem 2. Assume some regularity conditions andn21 C∆> 0; andm→∞nn0 −ntr(Σi ) − tr(Σi0 ) ii
lim sup+< 1,
n1 ni ∆C
ni ∆
m→∞

(9)

lim inf

i 6= i0 .

(10)Then, under (6), (3) holds (8).
SVM gives (1) and (2) for high-dimensional imbalanced data. We improved SVM by using a robustintercept to imbalanced data. Let Si = {j|α̂ij 6= 0, j = 1, . . . , ni } for i = 1, 2. Then we define robust SVM:
yr (x0 ; w) =X

α̂j tj xTj x0 +j6=j 0 (j,j 0 ∈S1 )j∈SNote thatPn1j=1 α̂1j =Pn2Xα̂1j xT1j x1j 0−
2(nS1 − 1)X
j6=j 0 (j,j 0 ∈S2 )α̂2j xT2j x2j 0
2(nS2 − 1)(11)

j=1 α̂2j (= A,b say) from (4). We have the following theorem.

Theorem 3. Under some regularity conditions, (11) holds (8).
RSVM gives (8) for imbalanced data and any C satisfying (9) without assuming (7) and (10).
In this talk, we generalized RSVM for kernel functions and discussed the choice of tuning parameters. Finally, we checked the performance of RSVM and the validity of the tuning parameter by numerical simulations
and actual data analysis.

References[1] Nakayama, Y. (2019). Robust support vector machine for high-dimensional imbalanced data. Communications in Statistics - Simulation and Computation, in press (doi: 10.1080/03610918.2019.1586922).
[2] Nakayama, Y., Yata, K., and Aoshima, M. (2019). Bias-corrected support vector machine with gaussian
kernel in high-dimension, low-sample-size settings. Annals of Institute of Mathematical Statistics, inpress (doi: 10.1007/s10463-019-00727-1).

2

報告書：多様な分野における統計科学に関する諸問題
成蹊大学名誉教授
新村秀一「高次元 Microarray データを用いて癌遺伝子の特定と癌の亜種を見つける研究」が 1970 年から行われてきた。これらの研究で用いられたデータが公開されているので、統計に限らず工学系の機械学習（AI）、
パターン認識、Bio 工学の新テーマとして研究されてきた。本研究に関する医学以外の多くの研究が「3つの困難」を指摘していることが示す通り、いずれも成功していない。
しかし、仮に症例数が n=100 で遺伝子数が p=10,000 の発現量データとすれば、2 群判別が最も適した手法である。筆者は 2015 年 10 月 28 日から 12 月 25 日の僅か 54 日間で、簡単にこの問題を解決した。用いたデータは、1999 年から 2004 年の間に米国の 6 医学研究プロジェクトが Science や The New
England Journal of Medicine に論文を発表し、研究に用いた Microarray の公開データである。これらはAlon と Singh は癌と健常の 2 クラス、Golub らの 4 種は異なった 2 種の癌のクラスである。異なった癌種の判別であるが、ほぼ共通の結果が得られた。
結果は非常に単純である。6 種のデータは線形分離可能なデータ(Linearly Separable Data, LSD)である（判別分析の新事実 Fact3）。この重要な事実であり信号が、これまでの研究では指摘していない。唯一、青嶋、矢田、石井らが、これ等の 6 種を含む 10 種以上の Microarray データで 2 群が 2 つの異なった球
体上に布置していることを高次元 PCA 等の研究で指摘している。さらに、筆者が開発した Matryoshka
Feature SelectioNMethod (Method2)で簡単に線形分離可能な含まれる遺伝子数が n 個以下の遺伝子のｋ組の部分空間(SMall Matryoshka, SM)と最小誤分類数(Minimum Number of Misclassifications, MNM)が1 以上の雑音の遺伝子の部分空間に分割できた（判別分析の新事実 Fact4）。各 SM は統計分析が得意とする小標本であるが、ロジスティック回帰以外、線形分離可能な事実が示されなかった。そこで、MNM基準による改定 IP-OLDF（RIP）の判別スコア(RIP Discriminant Score, RipDS)を変数とし、n*k 次元の信号データ(k<=n)を作成した。これを一元配置の分散分析、ｔ検定、相関分析、クラスター分析、PCA
で癌の遺伝子診断と Malignancy Indexes を世界で初めて提案でき、今年 5 月に Springer から 2 冊目の ”High-dimensional Microarray Data Analysis- Cancer Gene Diagnosis and Malignancy Indexes byMicroarray”を出版した。
以 上の 研 究 が 簡 単 に で き た の は 、 大 学 卒 業 以 来 行 っ て き た 判 別 分 析 の 新 理 論 (New Theory ofDiscriminant Analysis After R. Fisher, Springer 2016)が 2015 年に完成し、新理論がその応用問題として
Microarray を用いた癌の遺伝子解析（判別分析の問題 5）を簡単に解決できたためである。本来であれば癌の遺伝子研究の専門家でない筆者が「癌の遺伝子診断」までを行うことは適していないことは十分理解している。しかし、癌は遺伝子の病気であり、高次元の Microarray 空間で 2 群が完全に分かれていて、
かつ MNM＝０である k 個の SM に分割できる。そして信号データを作ることで、統計手法で線形分離可能であり有効と考えられる結果が数多く発見できた。しかし、これ等のどの結果が医学的に役に立つか否かは医学専門家の検証が必要である。残念ながら Golub らの研究後に、「NIH が乳がん以外の癌に関してこの種研究は成果が出ないと判断し、医学研究が終わったようである」。今、いかに医学専門家の

検証につなげるかを模索している。しかし NIH の報告を知らずに研究を続けている工学研究者の在り方
は、一般的に問題であろう。また、データが LSD であるのに、そのデータを学習標本に用いる AI 研究
が LSD の事実を指摘しない点だけが、まだ説明できていない。
大学卒業以来の研究テーマである判別の新理論を確立し、その応用として「高次元 Microarray の癌の遺伝子解析と診断」に成功した。そこでこれまでの研究を見直した結果、LSD である高次元データは、
ケース数 n 個以下の遺伝子のｋ組の小標本に必ず分割できるという事実が統計にとって一番重要と考えた。すなわち、我々は「高次元データの呪いから数理計画法(MP)の LP と IP で解放される。そして、分割された SM を統計分析することでこれまで見えてこなかった新しい世界が広がる（QP で定式化した
SVM ではできない）」ことになる。以上の重要なテーマを以下の 4 回のシンポジュームで報告したい。
1． 九州大学：「高次元遺伝子解析の呪いからの解放 1 -統計が 1970 年からこの問題を解決できなかった理由-」
2． 新潟大学：「高次元遺伝子解析の呪いからの解放 2 -癌の遺伝子診断-」
3．東京工業大学：「高次元遺伝子解析の呪いからの解放 3 -機械学習などの工学研究の問題点-」
4．秋田大学：「高次元遺伝子解析の呪いからの解放 4 -高次元データの分割法の最新結果-」
本発表では、統計的判別関数がなぜ役に立たなかったかの理由を報告する。正規分布を仮定した分散
共分散行列に基づく判別関数は、得点で試験の合否判定を行うと数学Ⅰa やⅡｂで誤分類確率が 2 割を超えるものがあることを既に報告した。そして Microarray では 17％も誤判別する。すなわち LSD 判別には全く役に立たないばかりか（判別分析の Problem2）、判別分析の基本統計の誤分類数(NM)も問題が多い（判別分析の Problem1）。このため、SM という小標本の判別分析でも分散共分散行列に基づく判別関
数の Fisher の LDF，QDF，正則化判別関数（RDA）、LASSO は全く役に立たない。これに対して、最尤推定法で求めたロジスティック回帰は、Microarray では役に立たないが、全ての SM を NM＝０で判別した。LASSO 研究者は、癌の遺伝子解析への応用を考えているが未だ成功していない。彼らは回帰係
数や判別係数を 0 にした非ゼロのモデルを用いれば、変数選択の NP-Hard を克服でき問題が解決できると単純に考えている。MNM 基準による RIP では、筆者の学位論文に示す通り、LSD でない Fisher のア
イリスデータで自然に多くの判別係数が 0 になることを示している。重要なのは LSD であるかどうかであり、次に判別係数が 0 になるかどうかという階層構造の優先度が全く理解されていない。ハードマー
ジン最大化 SVM(H-SVM)と RIP だけが理論的に LSD を判別できる。筆者の研究では、筆者の開発した3 種の最適判別関数(OLDF)とソフトマージン SVM(Ｓ-SVM)も LSD であることを確認した。医学研究
でＳＶＭを用いた研究もあるが Fact3 を報告したものがないことを考察する。そして残念なことに、RIPと Revised LP-OLDF だけが、ｋ組の SM と雑音に簡単に分割できるが、2 次計画法(QP)を用いた SVMができない理由を数理計画法の入門的な知識で説明する。一方、RIP が分割できる理由を筆者の定義したIP と LP モデルの実行可能領域と、線形代数の初歩的知識である連立方程式の解という初歩的知識の組
み合わせで明快に説明する。すなわち 1970 年から医学研究者の努力にかかわらず、また多くの統計や工学研究者の研究でも成果が出なかった理由を考察する。そして何故 MP の IP と LP モデルだけが高次元
データの呪いから我々を解放できたかを説明する。これによって、今後高次元の他のデータであっても、容易に統計分析の研究対象になることを説明する。研究は質が高く、2 群が LSD であるという検証しや
すいデータを用いたことによって、はじめて役に立つ研究を退官後に完成できたことは医学データを研究対象としたことが幸運であった。

Kolmogorov-Smirnov Tests
for Laplace Spectral Density Kernels
Yuichi Goto
(Waseda University)
Marc Hallin(Université libre de Bruxelles)
Masanobu Taniguchi (Waseda University)

1. Introduction
A new type of spectral density, called the Laplace spectral density kernel, was proposedby [1]. Laplace spectral density kernels characterize the collection of all marginal
bivariate distributions of a given stationary stochastic process, without requiring any
moment assumptions. In this presentation, we proposed a Kolmogorov-Smirnov (KS)test for Laplace spectral density kernels.

2. Setting and Main Results
Let {Xt : t ∈ Z} be a strictly stationary process with continuous and strictly increasing marginal distribution function F (·). We assume {Xt } satisﬁes the following
assumption.
Assumption 1. There exists constants ρ ∈ (0, 1) and K ∈ R such that, for arbitrary
intervals A1 , . . . , Aℓ ⊂ R and arbitrary t1 , . . . , tℓ ∈ Z,Cum I{Zt1 ∈A1 } , . . . , I{Ztp ∈Aℓ } ≤ KP(Zt ∈ A1 ) . . . P(Zt ∈ Aℓ )ρmaxi,j∈{1,...,ℓ} |ti −tj | .
The Laplace cross-covariance kernel of {Xt } is deﬁned, for all k ∈ Z and (x1 , x2 ) ∈R , as
cx1 ,x2 (k) := cF (x1 ),F (x2 ) (k) :=Cov I{Xt−k ≤x1 } , I{Xt ≤x2 } = Cov I{Ut−k ≤F (x1 )} , I{Ut ≤F (x2 )} ,2

where Ut := F (Zt ). Associated with this cross-covariance kernel, the Laplace spectral
density kernel is deﬁned, for (x1 , x2 ) ∈ R2 and λ ∈ [−π, π], as
1 X
fx1 ,x2 (λ) := fF (x1 ),F (x2 ) (λ) :=cx ,x (k) exp{−ikλ}.
2π k=−∞ 1 2∞

In order to estimate the Laplace spectrum of {Xt }, the Laplace periodogram is deﬁned, for (x1 , x2 ) ∈ R2 and λ ∈ [−π, π], asF (x ),F (x2 )Inx1 ,x2 (λ) := In,U 1(λ) :=

1 x1
d (λ)dxn2 (−λ),
2πn n

where
dxn (λ) :=

n−1
X
s=0

I{Xs ≤x} exp{−iλs} =n−1
X

I{Us ≤F (x)} exp{−iλs}.

s=0This research supported by Grant-in-Aid for JSPS Research Fellow Grant Number JP201920060 (Y.
Goto), and the Research Institute for Science & Engineering of Waseda University and JSPS Grantin-Aid for Scientiﬁc Research (S) Grant Number JP18H05290 (M. Taniguchi).

Generalized maximum composite likelihood estimators fordeterminantal point processes
Kou Fujimori*, Sota Sakamoto and Yasutaka Shimizu
Waseda University

In this talk, we dealt with the estimation problem for determinantal point processes (DPPs).
DPPs are the classes of spatial point processes with repulsive properties for each pair of realized
points and studied intensively in terms of statistical physics to capture the behavior of fermions.
Let (Ω, F, P) be a probability space. For the kernel function K : Rd × Rd → R, the point
process X called determinantal point process with kernel K if the measure;


̸=
∑
E
1{x1 ∈A1 ,...,xp ∈Ap }  , p = 2, 3, . . .
(x1 ,...,xp )∈X p

where Aj , j = 1, . . . , p are bounded Borel sets on Rd and X p is p-direct product of the point
process X, has the following density function (p-th order joint intensity function)
ρ(p) (x1 , . . . , xp ) = det[K](x1 , . . . , xp ),with[K](x1 , . . . , xp ) := (K(xi , xj ))1≤i,j≤p ∈ Rp×p .
We consider the case when the kernel K satisﬁes the following form:√
K(x, y) = ρ(x, λ)ρ(y, λ)Cα (x − y),where the function ρ(·, λ) is the intensity function with an unknown parameter λ > 0 and Cα is
the pair correlation function which satisﬁes Cα (0) = 1 with an unknown parameter α ∈ Rq . We
denote this kernel K = Kθ , where θ = (λ, α).
Our goal was to construct estimators for unknown parameter θ. Since well-known point
processes such as the Gibbs point process or Cox process only have the exact form of joint
intensities of second order, we often use the second order estimation function such as composite
likelihood functions. On the other hand, joint intensities of DPPs are given by the determinantof positive deﬁnite kernels, which allows us to compute the joint intensities of general order.
Therefore, we can introduce the two step composite likelihood approach by using p-th order joint
intensity function ρ(p) based on observations from ﬁxed window Dn ⊂ Rd which is the generalizedversion of 2nd order estimation function. For the intensity parameter λ, we considered thefollowing normalized quasi-likelihood function.
∫∫
Hn1 (λ) =log ρ(x, λ)N (dx) −ρ(u, λ)du,Dn

Dnwhere Dn ⊂ Rd is the observation window centered 0 which satisﬁes the following condition:|Dn | ≍ nd ,µd−1 (∂Dn ) ≍ nd−1 ,n→∞with | · | and µd−1 (·) is the d and d − 1-dimensional Lebesgue measure, respectively. We deﬁned
the following normalized p-th order composite likelihood function for every integer p ≥ 2 toestimate the parameter α:
∫ {
}
(p)
(p)
Hn2 (λ, α) :=log[ρθ (x1 , . . . , xp )] − log[Kw,p (r : θ)]pDn×wr (x1 , . . . , xp )N (p) (dx1 · · · dxp ),

(1)

(p)

where r > 0 is a tuning parameter, ρθ and Kw,p are respectively the joint intensity of p-thorder of DPP(Cθ ) and the modiﬁed K-function of p-th order:(p)ρθ (x1 , . . . , xp ) = det[Kθ ](x1 , . . . , xp ),p∏=ρ(xi , λ) det[Cθ ](x1 , . . . , xp )
i=1

∫

and
Kw,p (r : θ) =

p
Dn

(p)

ρθ (x1 , . . . , xp )wr (x1 , . . . , xp )dx1 · · · xp ,where wr is a bounded weight function whose support is given by
Srp := {(x1 , . . . , xp ) : |x1 − xj | ≤ r, 1 ≤ j ≤ p}.
For example, a simple choice of the weight function wr is given by
wr (x1 , . . . , xp ) = 1Srp (x1 , . . . , xp ).
Note that N (p) is a counting measure of p-th order induced by the point process X, i.e.,


p
̸=
∏
∑
N (p) Aj  =
1{x1 ∈A1 ,...,xp ∈Ap } ,j=1(x1 ,...,xp )∈X p

where Aj , j = 1, . . . , p are bounded Borel sets on Rd . Using the estimating functions Hn1 (λ)
and Hn2 (λ, α), we deﬁned the following two-step estimator for θ = (λ, α).
(p)

(p)

Deﬁnition 1. The estimator θ̂n = (λ̂n , α̂n ) is called generalized maximum composite likelihood
estimator if
λ̂n := arg sup Hn1 (λ),

(2)

λ∈Θλ
(p)

α̂n(p) := arg sup Hn2 (λ̂n , α).

(3)

α∈Θα

In particular for stationary case, i.e., the case when the kernel function satisﬁes the following
condition;Kθ (x, y) = λCα (x − y), x, y ∈ Rd ,we proved the consistency and the moment convergence of the estimators as the volume of theobservation window |Dn | → ∞ under some suitable conditions. Moreover, we presented theinformation criterion based on the proposed estimating method in this talk.

2














                 

     

          !               ! #     $     % % & &'                                 "            "    !

        % &          "   !

  (     )    ! %,---&      !     

! %,---&        ¾ %1 ,&! )   +  +%,--.&  )    !  *                %)& % & /2   %,-34&

"   

 

%

 





%



&



  











     %    ! 5  / 





 %¡&   7  !

 

  

&










 




  


%




1

%
&

&   

 



¾ %1 &



%1 &





6

,

 

6

,



,
6



&

,         %6 6 &6  ! %,---&      ! 8

%



&

 0







  

  

 



   
%

 %&  % & %  &&

5$   7   %&            







 
 






'   $  

6 &



              ,  

  !

   
 $

 

6  %6

µ

 

 


,

 %&



 


,






 



 
 
 



%

&



  
















 


,





% & 

 



 % & 9










 % &

 


 
    % & )    ! !     ! ! %,-34&!       $           10$ :66;:<.!

$ '!  2   $ =!  2 $ ! %,---&!              
? $ 2! 

5 @@$ )

   A
5 $

7       4$ >.,;>4.!

(! %611,&!                 7$ 4.<;4-.!

! %61,<&! *     B           #0        16$ ,.<;,-6!

6 

A closed form EM algorithm for a multivariate
skew-normal model
Toshihiro Abe (Nanzan University)
Hironori Fujisawa (The Institute of Statistical Mathematics, RIKEN)
Takayuki Kawashima (Tokyo Institute of Technology)
The most famous skew distribution is a skew-normal distribution, which was proposed by
Azzalini (1985, 1986). Various extensions to multivariate skew distributions have been proposed.
Azzalini & Dalla Valle (1996) introduced a two-dimensional skew distribution whose marginal
distribution is a univariate skew-normal distribution through additional parameters. Azzalini& Capitanio (1999) proposed a multivariate skew-normal distribution and investigated many
properties of multivariate skew-normal distributions. Arellano-Valle & Genton (2005) presented
some expressions of skew distribution and discussed basic classes of skew distribution, including
skew-symmetric, skew-elliptical, and skew-spherical distributions. It is known that the maximum likelihood estimation for skew-normal distributions presents a problematic aspect in theneighborhood of base distribution.
In this talk, we consider the maximum likelihood estimation via an EM algorithm (Dempsteret al., 1977) for a multivariate skew-normal distribution. Some stochastic representations were
proposed to obtain a multivariate skew-normal distribution. A stochastic representation often
enables us to make an EM algorithm in a convenient manner (McLachlan & Krishnan, 2007). Lin
et al. (2007) constructed an EM algorithm for a mixture of univariate skew-normal distributions,but the EM algorithm obtained there demands a solution of complicated estimating equation
for a skew parameter. The same problem happens even in a usual univariate skew-normal
distribution, so that it is hard to extend the idea to a multivariate case. Lin (2009) considered an
EM algorithm for a mixture of multivariate skew-normal distributions, although the multivariate
skew-normal distribution treated there was based on Sahu et al. (2003) and a diﬀerent typefrom Azzalini & Capitanio (1999). The EM algorithm demands multi-dimensional numericalintegrations, which are hard to compute in a high-dimensional case. So far, there has no EM
algorithm in a closed form.
Here we review the stochastic representation which implies the multivariate skew-normal
distribution (Azzalini & Capitanio, 1999). Suppose
)
(
)
(
Ω δ
Y
∗
∗
∼ Np+1 (0, Ω ), Ω =
.
(1)
δT 1
Y0
Then, U = sgn(Y0 )Y has the following density function:f (u) = 2Φ(αT u)ϕ(u; 0, Ω),(2)where α = Ω−1 δ/(1 − δ T Ω−1 δ)1/2 , Φ(z) is the cumulative density function of the standard normal distribution and ϕ(z; 0, Ω) is the normal density function with mean vector 0 and covariancematrix Ω. The random variable U with the density (2) is said to have the multivariate skewnormal distribution (with mean zero). The distribution of U + µ is expressed as SNp (µ, Ω, α).
In this talk, we construct an EM algorithm in a closed form, using the stochastic representation (1) with a slight modiﬁcation. We change the covariance structure of (Y , Y0 ) from Ω∗to()Ωτ Ω1/2 δΣ=(3),τ δ T Ω1/2
τ2

and then we consider the random variable X = sgn(Y0 )Y . Here we give two devices. The ﬁrst
one is essential when we construct an EM algorithm. We add the parameter τ in (3). We canconsider the covariance matrix Σ on the whole space of positive deﬁnite matrix, although Ω∗ has
the restriction τ = 1 in (1). This device makes it easy to optimize the objective function in the
M-step of the EM algorithm. The second one is a slightly diﬀerent expression of parameter. We
replace δ in Ω∗ by Ω1/2 δ in Σ, which changes the positive deﬁnite condition of the covariancematrix from δ ⊤ Ω−1 δ < 1 to the simple condition δ ⊤ δ < 1. We will show that
√ X = sgn(Y0 )Y
has a multivariate skew-normal distribution SNp (0, Ω, α) with α = Ω−1/2 δ/ 1 − δ ⊤ δ, which isthe same as the original multivariate skew-normal distribution if Ω1/2 δ is placed back to δ. A
remarkable point is that the distribution of X does not depend on the parameter τ , since thejoint distribution of (Y , Y0 ) includes the parameter τ . In this sense, the parameter τ can becalled the overparameter. The overparameter τ is not necessary when we consider a multivariate
skew-normal distribution, but it enables us to construct an EM algorithm in a closed form. The
R package snem has been developed to obtain the maximum likelihood estimate via the EM
algorithm.

References
Arellano-Valle, R.B. & Genton, M.G. (2005). On fundamental skew distributions. Journal ofMultivariate Analysis 96 (1), 93–116.
Azzalini, A. (1985). A class of distributions which includes the normal ones. Scandinavian Journal of Statistics 12 (2), 171–178.
Azzalini, A. (1986). Further results on a class of distributions which includes the normal ones.
Statistica 46 (2), 199–208.
Azzalini, A. & Capitanio, A. (1999). Statististial applications of the multivariate skew normal
distribution. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 61(3), 579–602.
Azzalini, A. & Dalla Valle, A. (1996). The multivariate skew-normal distribution. Biometrika
83 (4), 715–726.
Dempster, A.P., Laird, N.M. & Rubin, D.B. (1977). Maximum likelihood from incomplete data
via the em algorithm. Journal of the Royal Statistical Society. Series B (Methodological) 39,
1–38.
Lin, T.I. (2009). Maximum likelihood estimation for multivariate skew normal mixture models
Journal of Multivariate Analysis 100 (2), 257–265.
Lin, T.I., Lee, J.C. & Yen, S.Y. (2007). Finite mixture modelling using the skew normal distribution. Statistica Sinica 17, 909–927.
McLachlan, G. & Krishnan, T. (2007). The EM Algorithm and Extensions. John Wiley & Sons.
Sahu, S.K., Dey, D.K. & Branco, M.D. (2003). A new class of multivariate skew distributions
with applications to Bayesian regression models. Canadian Journal of Statistics 31 (2), 129–
150.

Asymptotics of maximum likelihood estimation for stable
law with continuous parameterization
Nanzan University

Muneya Matsui

Abstract Asymptotics of maximum likelihood estimation for α-stable law are analytically investigated with a continuous parameterization. The consistency and asymptotic
normality are shown on the interior of the whole parameter space. Although these asymptotics have been provided with Zolotarev’s (B) parameterization, there are several gaps
between. Especially in the latter, the density, so that scores and their derivatives are
discontinuous at α = 1 for β ̸= 0 and usual asymptotics are impossible. This is considerable inconvenience for applications. By showing that these quantities are smooth in thecontinuous form, we fill gaps between and provide a convenient theory. We numericallyapproximate the Fisher information matrix around the Cauchy law (α, β) = (1, 0). The
results exhibit continuity at α = 1, β ̸= 0 and this secures the accuracy of our calculations.

Main contents We use the following definitions and notations throughout. Denote ch.f.
of stable law in continuous form by

)
(
}
{
 exp − |σt|α 1 + iβ signt tan πα (|σt|1−α − 1) + iµt
if α ̸= 1
2()
φ(t) =
 exp − |σt| − iσt (2β/π) log |σt| + iµt
if α = 1,
where µ ∈ R, σ ∈ R+ , α ∈ (0, 2] , β ∈ [−1, 1] with R+ = (0, ∞). We denote thisparameter space by ΘM and its interior by Θ◦M . A parameter vector is denoted byθ = (θ1 , θ2 , θ3 , θ4 )′ = (µ, σ, α, β)′ . As usual f ′ , f ′′ mean the first and the second derivatives
with respect to (w.r.t.) x and fθ = (fθ1 , fθ2 , fθ3 , fθ4 )′ denotes a vector of partial derivativesof f w.r.t. θ. The second order partial derivatives w.r.t. x and θ are denoted byfθ′ i =

∂2f
∂2f
=
,
∂x∂θi
∂θi ∂xfθi θj =

∂2f
∂2f
=
,
∂θi ∂θj∂θj ∂θi

i, j = 1, . . . , 4,i.e. all derivatives will be shown to be interchangeable in our case. Moreover, denote thelog-likelihood function of f and its scores respectively byℓ(x; θ) = log f (x; θ)andℓθi (x; θ) =∂ℓ(x; θ)fθ (x; θ)= i,∂θif (x; θ)1 ≤ i ≤ 4,

where for convenience we sometimes write ℓ(x) and ℓθi (x) for these quantities. The second
order derivatives of score functions w.r.t. θ and x, denoted by( ′)1∂ℓθi (x; θ)= 2fθi (x; θ)f (x; θ) − fθi (x; θ)f ′ (x; θ) ,∂x
f (x; θ)()∂ℓθi (x; θ)1ℓθi θj (x) == 2
fθi θj (x; θ)f (x; θ) − fθi (x; θ)fθj (x; θ) ,∂θjf (x; θ)ℓ′θi (x) =1 ≤ i, j ≤ 4,

are also investigated. Again orders of partial derivatives w.r.t. (x, θ) are all exchangeable.
Based on the following technical Lemma 0.1 and Proposition 0.2, we establish asymptotics of maximum likelihood estimation as Theorem 0.3.
Lemma 0.1. For every x ∈ R, f (x; θ) : θ ∈ Θ◦M is twice continuously diﬀerentiable w.r.t.
θ, and fθi , i = 1, . . . , 4 is continuously diﬀerentiable w.r.t. x. Moreover fθi , fθ′ i , fθi θj i, j =1, . . . , 4 are jointly continuous in (x, θ) on R × Θ◦M .

The tails of f and its derivatives for suﬃciently large |x| satisfy
f = O(|x|−(1+α) ),fµ = −f ′ = O(|x|−(2+α) ),fσ = O(|x|−(1+α) ),fα = O(|x|−(1+α) log |x|),fβ = O(|x|−(1+α) ),fµ′ = −fµµ = O(|x|−(3+α) ),fσ′ = −fµσ = O(|x|−(2+α) ),fα′ = −fµα = O(|x|−(2+α) log |x|),fβ′ = −fµβ = O(|x|−(2+α) ),fσσ = O(|x|−(1+α) ),fσα = O(|x|−(1+α) log |x|),fσβ = O(|x|−(1+α) ),fαα = O(|x|−(1+α) log2 |x|),fαβ = O(|x|−(1+α) log |x|),fββ = O(|x|−(1+α) ).

Furthermore for α = 1, β ∈ (−1, 1), we havefσσ = O(|x|−3 log |x|),fββ = O(|x|−3 log |x|).

Proposition 0.2. Let θ ∈ Θ◦M . For every x ∈ R,ℓθi (x), ℓ′θi (x)andℓθi θj (x),i, j = 1, . . . , 4,

are well-defined and continuous in θ, and they are jointly continuous in (x, θ) on R × Θ◦M .
Concerning tail behaviors, we have for suﬃciently large |x|, x ∈ R,ℓµ (x) = O(|x|−1 ),ℓσ (x) = O(1),ℓα (x) = O(log |x|),ℓβ (x) = O(1),ℓµµ (x) = −ℓ′µ (x) = O(|x|−2 ),ℓµσ (x) = −ℓ′σ (x) = O(|x|−1 ),ℓµα (x) = −ℓ′α (x) = O(|x|−1 log |x|),ℓµβ (x) = −ℓ′β (x) = O(|x|−1 ),and moreover,ℓσσ (x) = O(1),ℓσα (x) = O(log |x|), ℓσβ (x) = O(1),ℓαα (x) = O(log2 |x|), ℓαβ (x) = O(log |x|), ℓββ (x) = O(1).
Theorem 0.3. Let θ̂n be the maximum likelihood estimator based on i.i.d. n observations
from stable law (Pθ : θ ∈ ΘM ). Assume that the true parameter θ0 is in the interior
θ0 ∈ Θ◦M and prepare an arbitrary compact set C ⊂ Θ◦M such that θ0 ∈ C. Then MLEθ̂n restricted on C is consistent and has asymptotic normality. In particular we have an
expression
n
√
1 ∑
√
ℓθ0 (Xk ) + oPθ0 (1),n (θ̂n − θ0 ) = Iθ−1
0
n
k=1

where

√
d
n (θ̂n − θ0 ) → N (0, Iθ−1
) as n → ∞, and Iθ0 is the Fisher information matrix.
0

Other remaining results including numerical works are given in the reference. By
applying our present work, we are examining Quasi-maximum likelihood estimation forstable law.

References[1] Matsui, M. (2019) Asymptotics of maximum likelihood estimation for stable law with continuous parameterization. arXiv:1901.09303.

2

Pitman Closeness Domination in Predictive Density Estimation for Two Ordered Normal
Means Under α-Divergence Loss
Yuan-Tsung Chang (Mejiro University) Nobuo Shinozaki (Keio University )William, E. Strawderman (Rutgers University)
We consider Pitman closeness domination in predictive density estimation problems when the underlying
loss metric is α-divergence {D(α)}, a loss introduced by Csiszàr (1967) given byZp̂(ỹ|y)Dα {p̂(ỹ|y), p(ỹ|ψ)} = fαp(ỹ|ψ)dỹ,(1)p(ỹ|ψ)where, for −1 ≤ α ≤ 1
 4
 1−α2 (1 − z (1+α)/2 ),fα (z) =z log z,− log z,|α| < 1α=1
α = −1.

(2)

Here KL loss corresponds to α = −1. The case α = 1 is sometimes referred to as reverse KL loss.
We consider the normal distributions case and generalized Bayesian predictive densities case.
1) Case of normal distributions
If the true density function of Y is N (µ, σ 2 ) and the estimated predictive density of Y is N (µ̂, σ̂ 2 ), Changand Strawderman (2014) have derived the general form of Dα loss and have shown that it is a concave
monotone function of quadratic loss and is also a function of the variances (observed, predicand, andplug-in). The general form is given as following:a) for −1 < α < 1,(µ̂−µ)24
22
2
2 −A(σ 2 ,σ̂ 2 )
2Dα (N (ỹ|µ̂, σ̂ ), N (ỹ|µ, σ )) =,(3)
1 − d(σ , σ̂ )e1 − α2whereσ (α−1)/2 τd(σ , σ̂ ) = (α+1)/2 ,
σ̂
2

2

2

2

A(σ , σ̂ ) =



1−α
2σ 2



(1 − α)τ 2
1−
2σ 2


> 0,

1
=
τ2




1+α 1−α
+
.
2σ̂ 2
2σ 2

Further, d(σ 2 , σ̂ 2 ) < 1 and A(σ 2 , σ̂ 2 ) > 0.

.



(µ̂−µ)2σ2σ2.
c) for α = −1, D−1 (N (ỹ|µ̂, σ̂ 2 ), N (ỹ|µ, σ 2 )) = 12σ̂ 2 − log σ̂ 2 − 1 +σ̂ 2b) for α = +1, D+1 (N (ỹ|µ̂, σ̂ ), N (ỹ|µ, σ )) = 12
2

2





σ̂ 2
σ̂ 2
σ 2 − log σ 2 − 1

2+ (µ̂−µ)
σ2

(4)
(5)

The underlying distributions considered are normal, including the distribution of the observables,
the distribution of the variable whose density is to be predicted, and the estimated predictive density
which will be taken to be of the plug-in type. We demonstrate {D(α)} Pitman closeness domination ofcertain plug-in predictive densities over others for the entire class of metrics simultaneously when related
Pitman’s closeness domination holds in the problem of estimating the mean.
Examples of Pitman closeness domination presented relate to the problem of estimating the predictive
density of the variable with the larger mean. More precisely, let X1 ∼ N (µ1 , σ12 ) and X2 ∼ N (µ2 , σ22 ) be
two independent random normal variables, where µ1 ≤ µ2 . Under the above restriction we wish to predicta normal population with mean equal to the larger mean, µ2 , and variance equal to σ 2 , Ỹ ∼ N (µ2 , σ 2 ).
We consider different versions of this problem, depending on whether the σi2 , i = 1, 2 are known or are
unknown but satisfy the additional order restriction, σ12 ≤ σ22 . The case of two ordered normal means
with known covariance matrix is also considered.
We also consider {D(α)} Pitman domination of certain generalized Bayesian (best invariant) procedures suggested by Corcuera and Giummole (1999) in next section.
2) Case of generalized Bayesian predictive densities
In this section we discuss improving the generalized Bayesian predictive densities suggested by Corcuera
and Giummole (1999) under D(α) loss.
Based on the data
Xij ∼ N (µi , σi2 ), i = 1, 2, j = 1, · · · , ni ,we predict the density Ỹ ∼ N (µi , σi2 ), i = 1, 2. We denote its density function by p(ỹ; µi , σi ), where µiand σi2 are unknown.

When −1 ≤ α < 1, Corcuera and Giummole (1999) have established that the best invariant predictive
density of p(ỹ; µi , σi ) based solely on xi1 , · · · xin2 is"!2 #−(2ni −1−α)/2(1−α)
1−αy − x̄i,(6)p̂α (ỹ; x̄i , σ̃i ) ∝ 1 +2ni + 1 − ασ̃iwhere x̄i is the sample mean and σ̃i2 = ((ni − 1)/ni )s2i is the sample variance. Corcuera and Giummole(1999) have also shown that p̂α (ỹ; x̄i , σ̃i ) is the generalized Bayesian predictive density for the priordensity f (µi , σi ) ∝ 1/σi , 0 < σi < ∞. It is to be noted that p̂α (ỹ; x̄i , σ̃i ) is not a normal distribution,although the plug-in density N (x̄i , s2i ) is the generalized Bayes rule when α = 1.
We consider the following two cases separately where order restrictions on µi and/or σi2 are present,i) when µ1 ≤ µ2 case and ii) when µ1 ≤ µ2 and σ12 ≤ σ22 .
OSwith µ̂CSor µ̂OSWe consider to improve p̂α (ỹ; x̄i , σ̃i ) or p̂α (ỹ; µ̂OSi ,ii , σ̃i ) by replacing x̄i with µ̂irespectively, where)n1 s22n2 s21OSGDµ̂1 = min X̄1 , µ̂≡X̄1 +X̄2 , µ̂OS
= max{X̄2 , µ̂GD },2
n1 s22 + n2 s21
n1 s22 + n2 s21

µ̂CS
=
1

µ̂CS
=
2

 OS
,
 µ̂1 12
X̄1 + n1n+nX̄2 , min X̄1 , n1n+n2
2

if s21 > s22

 OS µ̂2 ,

if s21 ≤ s22

if s21 ≤ s22

1
2
X̄1 + n1n+nX̄2 , max X̄2 , n1n+n2
2

if s21 > s22 .

The next lemma is usefully for improving the generalized Bayesian predictive densities (6).
Lemma Let f (·) be the probability density function of X ∼ N (0, τ 2 ). Assume
R ∞ that g(t) ≥ 0 issymmetric about the origin and is a strictly decreasing function of |t| such that −∞ g(x)f (x)dx < ∞.
R∞
Then −∞ g(y − x)f (y − µ)dy is a strictly decreasing function of |x − µ|.
Let µ̂i denote an estimator of µi , i = 1, 2 in general. Now we show that for any 1 ≤ α < 1,
Dα (p̂α (ỹ; µ̂i , σ̂i ), p(ỹ; µi , σi )) is a strictly increasing function of |µ̂i − µi |.
From Lemma, we see that for |α| < 1,Z ∞Dα (p̂α (ỹ; µ̂i , σ̂i ), p(ỹ; µi , σi )) ∝ 1 −g(ỹ − µ̂i )f (ỹ − µi )dỹ∞is a strictly increasing function of |µ̂i − µi |, where2 −(2ni −1−α)(1+α)/4(1−α)

y−x
1−α
g(y − x) = 1 +2ni + 1 − α
σ̂i
and



(1 − α)(y − µ)2f (y − µ) ∝ exp −.
4σ 2For α = −1,
D−1 (p̂−1 (ỹ; µ̂i , σ̂i ), p(ỹ; µi , σi )) = −Eỹp̂−1 (ỹ; µ̂i , σ̂i )logp(ỹ; µi , σi )is a strictly increasing function of |µ̂i − µi | from Lemma.
i) µ1 ≤ µ2 case
Theorem 1. The predictive density estimate p̂α (ỹ; µ̂OSi , σ̂i ), i = 1, 2 is closer to the predictive densityp(ỹ; µi , σi ) than p̂α (ỹ; x̄i , σ̂i ), respectively, under the {D(α)} metric for all −1 ≤ α < 1 and for every
estimator σ̂i if and only if µ̂GD is Pitman closer to µ than X̄i for all σ12 and σ22 when µ1 = µ2 = µ.
ii) µ1 ≤ µ2 and σ12 ≤ σ22 case
Theorem 2. The predictive density estimate p̂α (ỹ; µ̂CS
2 , σ̂2 ) is closer to the predictive density p(ỹ; µ2 , σ2 )than p̂α (ỹ; µ̂OS
,
σ̂
)
under
the
{D(α)}
metricforall−1≤ α < 1 and for every estimator σ̂22 .
2
2
CS
Theorem 3. The predictive density estimate p̂α (ỹ; µ̂1 , σ̂1 ) is not Pitman closer to p(ỹ; µ1 , s1 ) thanp̂α (ỹ; µ̂OS
1 , σ̂1 ) when µ2 − µ1 is sufficiently large, under the {D(α)} metric for all −1 ≤ α < 1 and for any
estimator σ̂12 .

Auto-Distance Covariance Function for Time Series Analysis
Konstantinos Fokianos
k.fokianos@lancaster.ac.uk
Department of Mathematics & Statistics, Lancaster University
There has been a considerable recent interest in measuring dependence by employing the concept of distancecovariance function, a new measure of dependence for random variables, introduced by Székely et al. (2007). This
tool has been recently defined to the context of multivariate time series by Zhou (2012), but without exploring
the interrelationships between the various time series components. In this work, we extend the notion of distancecovariance to multivariate time series by defining its matrix version.
We denote by {Xt : t = 0, ±1, ±2, . . . } a d-dimensional time series process, with components Xt;r , r =1, . . . , d. Suppose we have available a sample of size n, that is {Xt , t = 1, . . . , n}. We define the pairwise autodistance covariance function as a function of the joint and marginal characteristic functions of the pair (Xt;r , Xt+j;m ),(r,m)

for r, m = 1, . . . , d. Denote by φj

(r,m)

φj

(u, v) the joint characteristic function of Xt;r and Xt+j;m ; that isj ∈ Z,(u, v) = E [exp (i(uXt;r + vXt+j;m ))] ,(r,m)

and the marginal characteristic functions of Xt;r and Xt+j;m as φ(r) (u) := φj(r,m)(u, 0) and φ(m) (v) := φj(0, v)

respectively, where (u, v) ∈ R , and i = −1. The pairwise auto-distance covariance function (ADCV) between Xt;r2

2and Xt+j;m , Vrm (j), is defined as the positive square root of
(r,m)

2
Vrm
(j) =

1
π2

φj

Z
R2

(u, v) − φ(r) (u)φ(m) (v)
2

2

dudv,

2

|u| |v|j ∈ Z.

The auto-distance covariance matrix, V (j), is then defined bydV (j) = [Vrm (j)]r,m=1 ,j ∈ Z.

The pairwise auto-distance correlation function (ADCF) between Xt;r and Xt+j;m , Rrm (j), is a coefficient thatlies in the interval [0, 1] and also measures dependence and is defined as the positive square root of2Rrm(j) = p

2
Vrm
(j)
p
,
2 (0) V 2 (0)Vrrmm

for Vrr (0)Vmm (0) 6= 0 and zero otherwise. The auto-distance correlation matrix of Xt , is then defined asdR(j) = [Rrm (j)]r,m=1 ,j ∈ Z.

When j 6= 0, Vrm (j) measures the dependence of Xt;r on Xt+j;m . In general, Vrm (j) 6= Vmr (j) for r 6= m,
since they measure different dependence structure between the series {Xt;r } and {Xt;m } for all r, m = 1, 2, . . . , d.
1

Thus, V (j) and R(j) are non-symmetric matrices, but V (−j) = V 0 (j) and R(−j) = R0 (j). More properties can befound in Fokianos and Pitsillou (2018). The empirical pairwise ADCV, Vbrm (j), for j ≥ 0, is the non-negative squareroot of
2
Vbrm
(j) =

n−j
X
1
Ar B m ,2(n − j) t,s=1 ts ts

where Ar = Ats and B m = Bts are Euclidean distance matrices given byArts = arts − ārt. − ār.s + ār.. ,
P











Pn−j rPn−j rn−j r
rr
2
s=1 ats /(n − j), ā.s =t=1 ats /(n − j), ā.. =t,s=1 ats /(n − j) .
m
Bts
is defined analogously in terms of bm
ts = |Xt+j;m − Xs+j;m |. Fokianos and Pitsillou (2018) have shown that2
for a d-dimensional strictly stationary and ergodic process {Xt } with E |Xt;r | < ∞, for r = 1, . . . , d, then for allwith arts = |Xt;r − Xs;r |, ārt. =j ∈ Z,Vb (j) → V (j),almost surely, as n → ∞. In addition, under pairwise independence it holds that2nVbrm(j) → Z :=X

λk Zk2 ,kin distribution, as n → ∞, where {Zk } is an i.i.d sequence of N (0, 1) random variables, and (λk ) is a sequence ofnonzero eigenvalues.
Based on these results, we will be considering the following two problems:• Testing the iid hypothesis by employing the generalized spectral density approach as developed by (Hong,1999)• Developing consistent testing procedures for the detection of multiple change-points in a given stream ofdata.

References
Fokianos, K. and M. Pitsillou (2018). Testing independence for multivariate time series via the auto-distance correlation matrix.
Biometrika 105, 337–352.
Hong, Y. (1999). Hypothesis testing in time series via the emprical characteristic function: A generalized spectral density approach. J. Am. Stat. Assoc. 94, 1201–1220.
Székely, G. J., M. L. Rizzo, and N. K. Bakirov (2007). Measuring and testing dependence by correlation of distances. Ann. Stat. 35,
2769–2794.
Zhou, Z. (2012). Measuring nonlinear dependence in time-series, a distance correlation approach. J. Time. Ser. Anal. 33, 438–457.

2

報告書
入江 薫
（東京大学経済学部・講師）
科研費シンポジウム「多様な分野における統計科学に関する諸問題」において、三日目にあたる九月十六日の午前のセッションにて、「縮小事前分布と状態空間モデル」という演題で研究成果を発表した。本発表は Irie, K., “Bayesian Dynamic Fused LASSO”, arXiv
preprint をもとにしており、当日もこのタイトルを表示した。
発表は予稿の内容に沿って行い、上記の公表済みの研究について発表した。講演時間が十分にあることを踏まえ、背景となる先行研究を詳細にレビューしたうえで、統計的計算手法の導出（の概略）も含めて当研究を解説した。縮小事前分布の研究は時系列データのベ
イジアンモデリングの分野においてもさかんに行われているものの、先行研究のほとんどが動的パラメータの経時変化に対する縮小のみを考慮しており、通常の意味での縮小効
果、すなわち動的回帰係数パラメータをゼロに近づけるという意味での縮小効果をモデルに含めていない、という点を指摘した。一方で、複数の罰則項を用いて異なる縮小効果を同時に実現する Fused LASSO の方法は上記の問題を解決するように思えるが、動的パラメータの時間発展を記述するものではないため、予測分析に使用することができず実用上は問題の多いアプローチであることを説明した。このような状況で本研究は、Fused
LASSO のモデリングを用いたマルコフ過程を動的パラメータの事前分布として使用することを提案する。著しい特徴として、そのような事前分布を用いた状態空間モデルは潜在パラメータを導入することで条件付き線形ガウス型の状態空間モデルとなり、事後分布の計算が簡易なギブス・サンプラーで実行可能となるため、計算の段階においては効率性・スケーラビリティ・実行の簡便さにおいて多くの利点を持つ。以上のようなモデルの特性はシミュレーションデータを用いた分析で例証され、本講演では多数の図表を示しながら分析結果を紹介した。
当日の質疑応答では、シンポジウムの主旨にのっとり、本研究の分野であるベイジアン・モデリングの最近の研究課題にまでさかのぼって、非専門家の方に向けて解説を行う機会が得られた。非ベイズ的な l1 罰則項が二重指数事前分布に対応することはよく知られているが、事前分布のモデリングを発展されることでよりスパースな状況を想定する分布を構成する研究がなされており、特に Carvalho, Polson and Scott (2009, 2010)以降さかんに研究・利用されている馬蹄(horseshoe)事前分布の解説を行った。聴衆から多くの質問があって、ベイズ統計学においてはスタンダードになりつつあるこの事前分布のクラスについて詳細な解説を行うことができ、本シンポジウムの主旨に多少なりとも貢献することができたと考える。
また、当日は縮小推定に詳しい/専門とする先生が多く出席しており、研究の核となるモデリングの点について有意義な指摘・提案をしていただいた。たとえば、対立する二種類の縮小効果を調整する重みパラメータは時間に関して一定の定数であるが、この設定が動的なスパース性に対しては制約的であり、重みパラメータの動的なモデリングが望ましいとの指摘があった。この点については従来から認識していたが、もっぱら技術的な問題から無視していたため、今回の指摘は問題の重要さを実感する意味で非常に有意義であった。
また、モデリングに関する様々な提案がなされ、階層化、grouped LASSO、二階差分に対する fused LASSO など、魅力的な今後の研究課題として取り組ませていただくことになった。
発表した研究は現在査読を受けている。今回得られたコメントと同様の指摘を査読者から受けることも十分考えられる。本シンポジウムでのやりとりを活かして論文を改善し、早くに出版することを目指したい。

バランス型経時測定データにおけるExtended GMANOVA モデルの解釈と新たな推定法
中京大学 国際教養学部 永井 勇
各個体に対して経時的に測定することで得られるデータは経時測定データと呼ばれ, 様々な分野で収集・分析がされている. この経時測定データの分析においては, データに潜む経時的な変動 (経時変動) を上手く捉えることが一つの目的である. 本講演では, 全ての個体で測定時点が揃っている, つまり, n 個全ての個体に対して時点 t1 , . . . , tp で測定した経時測定データを考えた. このような経時測
定データはバランス型経時測定データと呼ばれる.
本講演では, このバランス型経時測定データにおける経時変動を t1 , . . . , tp の多項式などで捉えることを考えた場合, どのようにモデル化されるかを報告し, このモデルが Potthoﬀ and Roy (1964) により提案された次の GMANOVA (一般化多変量分散分析) モデルと対応している事をまず報告した;

Y = 1n m′ X ′ + AΞX ′ + E,
ここで, Y は各行が各個体の経時測定データからなる n × p 既知行列, 1n は全ての成分が 1 の n 次元ベクトル, m は q 次元未知ベクトル, X は p × q 既知行列 (詳しくは後述), A は各行が各個体の性別などからなる n × k 既知 (説明変数) 行列, Ξ は p × q 未知行列, E は n × p 誤差行列である. モデル化の話より, (1n m′ X ′ + AΞX ′ ) が Y の経時変動を表しており, 例えば X の j 行目を (1, tj , . . . , tq−1)とjすることが, 経時変動を測定時点 t1 , . . . , tp の (q − 1) 次多項式で推定することに対応している.
ただし, このモデルには, 全ての説明変数に対して同一の次数の多項式などしか用いることができないという問題がある. 一方で, 1n m′ X ′ の部分を除いた形のモデルに対しては, Kollo and von Rosen(2005) の Def. 4.1.3 で拡張が考えられている. そこで本講演では, 上記の GMANOVA モデルに Kolloe
ら (2005) の拡張したモデルを用いた次のモデル (拡張 GMANOVA モデル) を提案した;r∑
Y = 1nµ′X0′ +AiΞiXi′ + E,
i=1
ここで, µ は q0 次元未知ベクトル, Xi は p × qi 既知行列 (詳しくは後述; i = 0, 1, . . . , r), Ai は各行が各個体のいくつかの説明変数のグループからなる n × ki 既知 (説明変数) 行列 (i = 1, . . . , r), Ξi は

p × qi 未知行列 (i = 1, . . . , r) である (他の記号は GMANOVA モデルで用いている記号と同様の意味
r∑である). ここで本講演では, n −ki − 1 > 0 で各 Ai が中心化 (つまり, 1′n Ai = 0ki , 0ki は ki 次元
i=1

ゼロベクトル) とし, rank(Xi ) = qi , rank(Ai ) = ki とした. さらに, E の各行は独立で E[E] = 0n 0′p ,Cov[vec(E)] = Σ ⊗ In , Σ は rank(Σ) = p の p × p 未知行列とした.
本講演では, GMANOVA モデルと経時変動の対応を踏まえ, 経時測定データに対する分析において提案したモデルの解釈について, 以下のように報告した. GMANOVA モデルと同様に, このモデル)(
r∑Ai Ξi Xi′ が Y の経時変動を表している. さらに, 平均的な経時変動を表においても 1n µ′ X0′ +i=1

す 1n µ′ X0′ と各説明変数のグループ Ai に対する経時変動を表す Ai Ξi Xi′ (i = 1, . . . , r) と分割して見ると, それぞれが異なる Xi (i = 0, . . . , r) を用いていることが分かる. ここで, GMANOVA モデルで

1

の経時変動の推定と同様に考えると, 例えば Xi の j 行目を (1, tj , . . . , tqj i −1 ) とすることが, 経時変動を測定時点 t1 , . . . , tp の (qi − 1) 次多項式で推定することに対応していることを報告した. これらをまとめると, 平均部分や各説明変数のグループごとに 異なる次数の多項式など を用いた経時変動の推定ができることが分かる. これにより, 説明変数のあるグループについての知見などを生かした経時変動が推定できる.
しかし, 従来の Kollo ら (2005) で提案されているモデルでの推定の際には, Xi や Ai について非常に厳しい仮定を置き, E の各行が従う確率分布が正規分布と仮定し, さらに反復計算も必要となる形となっている. 経時測定データの分析を考えた際には, 特に Xi や Ai に対する仮定をそのまま用いる事は, 説明変数に関する知見を生かした経時変動の推定などができなくなり, さらに用いる多項式にも非常に厳しい制約を課すこととなる. そこで本講演では, これらの仮定を置かずに以下の残差平方
和 (RSS) を最小にすることで得られる最小二乗推定量 (LSE) がどのような形になるかを報告した;

RSS(µ, Ξ1 , . . . , Ξr |Σ)

{(

def.

= tr

Y −1n µ′ X0′ −

r
∑

)

(

Ai Ξi Xi′ Σ−1 Y −1n µ′ X0′ −i=1

つまり, µ や Ξℓ の LSE µ̂, Ξ̂ℓ が次のように得られることを報告した;µ̂ =(X0′ Σ−1 X0 )−1 X0′ Σ−1 Y ′ 1n /n,








r
∑

)′ }
Ai Ξi Xi′

.

i=1

∑



Ξ̂ℓ=(A′ℓ Aℓ )−1 A′ℓ Y −Ai Ξ̂i Xi′  Σ−1 Xℓ (Xℓ′ Σ−1 Xℓ )−1 (ℓ = 1, . . . , r).





i̸=ℓ
しかしながら, Ξ̂ℓ においては, 他の Ξ̂i (i ̸= ℓ) が必要なため, Ξ̂i 全ての初期値や更新のための反復計算が必要である. そこで本講演では, Lütkepohl (1996) などにある vec 作用素と ⊗ (クロネッカー積)の関係, 直行射影行列の性質などを用いて, 次のように反復計算が不要な形での推定量を報告した;)(
r∑vec(Ξ̂ℓ ) = (Q′ℓ Qℓ )−1 Q′ℓ (Σ−1/2 ⊗ In )vec Y −Ai Ξ̂i Xi′ (ℓ = 1, . . . , r),
r
i>ℓ
∑
′
ここで,
Ai Ξ̂i Xi = O (適当な大きさのゼロ行列) としている. さらに Q1 = Σ−1/2 X1 ⊗ A1 ,Qj =

(

)

i>r
j−1

Inp −

∑

PQi

(Σ−1/2 Xj ⊗ Aj ) (j > 2), PQi = Qi (Q′i Qi )−1 Qi である (ここで, (Q′i Qi )−1 の存在

i=1は仮定する). この推定量では, まず vec(Ξ̂r ) が得られ, これを用いる事で vec(Ξ̂r−1 ) が得られる. これを繰り返す事で, r 個全ての推定量が反復計算をせずに得られる. この推定量を用いる事で, Kollo ら(2005) のモデルでの問題点 (Xi や Ai への厳しい仮定, E の各行への正規性の仮定, Ξ̂i を得るために反復計算が必要) を全て避けた新たな推定量が得られた.
さらに本講演では, 推定量の不偏性や Σ の不偏推定量についても報告した.
引用文献:[1] Kollo, T. & von Rosen, D. (2005). Advanced Multivariate Statistics with Matrices, Springer.
[2] Lütkepohl, H. (1996). Handbook of Matrices, John Wiley & Sons.
[3] Potthoﬀ, R. F. & Roy, S. N. (1964). A generalized multivariate analysis of variance model usefulespecially for growth curve problems. Biometrika, 51, 313–326.

2

Generalized Multinomial Distributions with Scalable Variance
Nobuaki Hoshino∗
September 2019

Abstract
In this talk we generalize the multinomial distribution to keep closure under the collapseof cells, motivated by a practice of statistical disclosure control. The resulting family of
distributions is characterized by Bell polynomials, and its parameter spaces are derived.
This family is shown to have the same mariginal first moment as that of the multinomial
distribution, but its second moment can be overdispersed. Also shown is that an infinite
dimensionl distribution of this family exists by Kolmogorov’s extension theorem. In this limit
the marginal distribution of positive frequencies is reduced to a simple random partitioning
distribution of a positve integer, which is useful to deal with a sparse contingency table.
Keywords: Discrete multivariate distribution, Overdispersion, Sparsity∗

School of Economics, Kanazawa University,
mail: hoshino@kenroku.kanazawa-u.ac.jp

Kakuma-machi,

1Kanazawa 920-1192,

Japan.

E-