科研費シンポジウム「融合する統計科学」報告書
本シンポジウムは科学研究費・基盤研究（A) 「大規模複雑データの理論と方法論の総合的研究」 （研究代表者：青嶋 誠（筑波大学）, 課題番号：15H01678）の助成により開催した。

■日時：2018 年 11 月 30 日～12 月 2 日■場所：金沢大学サテライトプラザ■開催責任者：星野伸明（金沢大学）＜

日時

プログラムと目次 ＞

演題

講演者

所属

頁

The Stirling and Eulerian

Xiaoling

早稲田大学

1

numbers in the Edo PeriodDou

Ewens 抽出公式に対する正佃 康司

東京大学 大3

11 月 30 日
13:00
13:40

規近似とポアソン近似

学院総合文化研究科

14:20

Posteriorsomesampling

from

間野 修平

non-exchangeable

統計数理研究5所

priors
15:00

休憩

15:10

統計的形状モデルを利用し広瀬 修

金沢大学理工た大規模点群レジストレー

7

研究域

ションとその高速化
15:50

Kolmogorov-Smirnov

test

Rizky

based on kernel estimation

Fauzi

院数理学府

石井 晶

東京理科大学

16:30

休憩

16:40

Tests for high-dimensional
covariancematrices

Reza

and

strongly

911

理工学部 情correlation matrices underthe

九州大学大学

報科学科

spikedeigenvalue model
17:20

癌の高次元遺伝子解析の諸
問題（1）
[i]

新村 秀一

成蹊大学

13

18:00終了

12 月 1 日
9:30

A Robust-filtering Method

国友 直人

for Small Sample Economic

明治大学政治15

経済学部

Time Series
10:10

ApplicationsofDistance

松井 宗也

Correlation to Time Series
10:50

割引因子ポアソン-ガンマ状

南山大学経営

17

学部
入江 薫

態空間モデルによる計数時

東京大学経済19

学部

系列モデルの逐次分析
11:30

昼休み

13:00

Testingforchangesin西埜 晴久

income inequality in Japan

広島大学社会20

科学研究科・
経済学部

13:40

Particle

Filtering

Non-linear

14:20

for

稗田 尚弥

State-Space

大学院工学研Models for Wind Speeds

究科経営工学and Directions

専攻

The Dantzig selector for a

藤森 洸

幹理工学部応processes

用数理学科

休憩

15:10

多様な関数を用いた経時変

永井 勇動の罰則付推定法

中京大学 国

22

24

26

際教養学部

傾向性仮説と変化点モデル

広津 千尋の様々な応用
16:30

早稲田大学基linear model of diffusion
15:00

15:50

東京理科大学

明星大学連携28

研究センターCharacterizationsofindicator

forfunctions青木 敏

神戸大学大学30

院理学研究科

fractional factorial designs
17:10

休憩

17:20

Pooling

incomplete

samples による統計解析[ii]

布能 英一郎

関東学院大学経済学部

32

18:00

Pitman's

Closeness

張 元宗

目白大学

34

野村 玲実

東京理科大学

36

Domination in Predictive
Density Estimation for Two
OrderedNormalMeans

Under α-Divergence Loss
18:40

終了

19:00

懇親会

12 月 2 日
9:30

10:00

一般欠測データの下での２標本問題における多変量正大学院理学研規母集団の同等性検定

究科 M1

単調欠測データが一様共分

佐伯 悠一郎

東京理科大学散構造を持つ場合の平行性

大学院理学研

仮説検定と水準差の信頼区

究科 M２

38

間
10:40

2-step 単調欠測データのも

川崎 玉恵とでの部分平均ベクトルの東京理科大学40

理学部

検定に対する修正尤度比検
定統計量
11:20

Recent cylindrical models
and their application totree data set

12:00

終了

[iii]

阿部 俊弘

南山大学理工
学部

42

The Stirling and Eulerian numbers in the Edo Period
Xiaoling Dou ∗, Hsien-Kuei Hwang †Similar to the recurrence relation satisfied by the binomial coefficients
( ) (
) (
)
n
n−1
n−1
=
+
,
k
k−1
k
[ ]
{ }
the Stirling numbers of the first kind nk and those of the second kind nk can be computed by
the recurrences:
[ ] [
]
[
]
n
n−1
n−1
=
+ (n − 1)
k
k−1
k
{ } {
}
{
}
n
n−1
n−1
=
+k
,
k
k−1
k
⟨ ⟩
subject to proper boundary conditions. Similarly, the Eulerian numbers nk and the second
order version ⟪nk⟫ satisfy the patterns
⟨ ⟩
⟨
⟩
⟨
⟩
n
n−1
n−1
= (n − k)+ (k + 1)
k
k−1
k
n
n−1
n−1
⟪ ⟫ = (2n − k − 1) ⟪⟫ + (k + 1) ⟪
⟫.
k
k−1
k
Over the last five centuries or so, these numbers emerged naturally and were studied extensively in a large number of diverse areas, ranging from finite calculus and series summations to
combinatorial structures and computer algorithms, and from statistics to spline interpolations,to name a few.
While the history of the developments of these numbers in the West has been largely andfactually clarified, that in the East has remained mostly obscure. In this work, we aim to provide
more historical materials during the Edo Period concerning these numbers, and to speciallyshed further light on their evolution (including introduction and use) in the Wasan History.
The following table summarize our findings on Stirling and Eulerian numbers in the Wasan
History
during{n} the Edo Period, where we also list the closely connected Bell numbers Belln :=∑0≤k≤n k :{Belln }n≥1 = {1, 2, 5, 15, 52, 203, 877, 4140, 21147, 115975, 678570, · · · }.
∑ [ ]Note that k nk = n!, and that a table of these numbers, as well as a procedure (the recurrence)of computing them, already appear in Arima’s 1769 book Shūki Sanpo.
∗†

Faculty of Science and Engineering, Waseda University, Japan
Institute of Statistical Science, Academia Sinica, Taiwan-1-Numbers
[n]
{nk }
⟨nk ⟩
k

⟪nk⟫
Bell numbers

West
James Stirling (1730)
James Stirling (1730)
Leonhard Euler (1736)
Jekuthiel Ginsburg (1928)
G. Dobiński (1877)
Eric Bell (1938)

East
SEKI Takakazu (1640–1708)
SAKA Masanobu (1782)
MATSUNAGA Yoshisuke (1694–1744)
MATSUNAGA Yoshisuke (1694–1744)
MATSUNAGA Yoshisuke (1694–1744)
ARIMA Yoriyuki (1769)

Interestingly, unlike the early developments of these numbers in the West, which are mostly
computational and algebraic, those carried out by the Wasankas already are not merely computational but also were motivated by combinatorial problems, adding another rich dimension to
the diversity and usefulness of these numbers. We will present the combinatorial connectionsof these numbers to certain games frequently played during this Period.

References[1] Y. Arima (1769). Shūki Sanpo.
[2] E.T. Bell (1938). The iterated exponential integers. Annals of Mathematics. 39, 539–557.
∑∑ m ∑ nm[3] G. Dobiński (1877). Summirung der Reihe
nmn! nn!
für m = 1, 2, 3, 4, 5, . . ..
n!
Grunert’s Archiv. 61, 333–336.
[4] L. Euler (1736). Euler Archive, E055.
[5] J. Ginsburg (1928). Note on Stirling’s Numbers, Amer. Math. Monthly 35, 77–80.
[6] J. Stirling (1730). Methodus Differentialis.
[7] Y. Matsunaga (1694–1744). Sanpo Zenkei.
[8] M. Saka (1782). Sanpō Gakkai.
[9] T. Seki (1712). Katsuyo Sanpo. (M. Araki, Y. Otaka).

-2-

Ewens 抽出公式に対する正規近似とポアソン近似
佃康司（東京大・総合文化）

Ewens 抽出公式は集団遺伝学の文脈で Warren John Ewens が 1972 年に導入した確率分
布 (Ewens, 1972) であり，多くの分野で議論されてきた．各分野における Ewens 抽出公式の役割については，Crane (2016) に解説がある．特に，確率の問題として古くから考えられてきたランダム置換と関連しているほか，ノンパラメトリック統計におけるベイズ法で事前分布として用いられる Dirichlet 過程からの標本確率分割が従う確率分布として統計学や機械学習といった分野で盛んに議論されてきた．
集団遺伝学における標本抽出理論において，標本（サイズを n(∈ N) とする）から計算される統計量として標本対立遺伝子分割における要素カウント (C1n , . . . , Cnn ) に興味があり，

Ewens 抽出公式はその確率を与える．ここで，要素カウントにおける Cjn (j ∈ N) はサイズ
n の標本において j 回の観測があったクラスの数を表している．Ewens 抽出公式には，母集団（を定常分布として生成する確率過程）のスケールされた突然変異率 θ と標本サイズn(∈ N) という二つのパラメータが存在する．これまで要素カウントや θ の完備十分統計量∑nである種類の数（Kn (= j=1 Cjn ) で定義される）などの統計量に対して多くの近似が示されてきた．中でも，特によく議論されてきた近似に Kn の正規近似とポアソン近似および(C1n , C2n , . . .) のポアソン近似がある．例えば，Arratia and Tavaré (1992) を参照．
これまでの Ewens 抽出公式に対する近似の多くは，二つのパラメータのうちの片方を固定して他方を無限大に発散するという極限をとる漸近論の設定を考えて導かれてきた．一方，両方のパラメータを同時に大きくするような漸近論も自然であり，Feng (2007) の 4 節で大偏差原理が議論されているものの，これまでよく議論されてきたような他の漸近的性質についてはあまり議論されてこなかった．講演では，Tsukuda (2017, 2018) に基づいて，要素カ
ウント (C1n , Cn2 , . . .) と種類の数 Kn について，n と θ がともに大きくなるような漸近論の設定のもとで講演者が導いた結果を紹介した．

-3-

参考文献
Arratia, R.; Tavaré, S. (1992). Limit theorems for combinatorial structures via discrete
process approximations. Random Structures Algorithms 3, no.3, 321–345.
Crane, H. (2016). The ubiquitous Ewens sampling formula. Statist. Sci. 31, no.1, 1–19.
Ewens, W.J. (1972). The sampling theory of selectively neutral alleles. Theoret. Population Biology 3, 87–112; erratum, (1972) Theoret. Population Biology 3, 240; erratum,(1972) Theoret. Population Biology 3, 376.
Feng, S. (2007). Large deviations associated with Poisson–Dirichlet distribution andEwens sampling formula. Ann. Appl. Probab. 17, no. 5–6, 1570–1595.
Tsukuda, K. (2017). Estimating the large mutation parameter of the Ewens sampling
formula. J. Appl. Probab. 54, no. 1, 42–54; correction, J. Appl. Probab. 55, no. 3, toappear.
Tsukuda, K. (2018). On Poisson approximations for the Ewens sampling formula whenthe mutation parameter grows with the sample size. Ann. Appl. Probab. to appear.

-4-

Posterior Sampling from some Non-Exchangeable Priors1
Shuhei Mano2
The Institute of Statistical Mathematics
January 10, 2019
Let the Dirichlet process with the measure θµ, where µ is a probability measure, be denotedby DP. Let F ∼ DP(θ; µ). It is well known that P(X1 ∈ ·) = µ(·) andP(Xn+1 ∈ ·|X1 , ..., Xn ) =Here, Λn (X1 , ..., Xn ) := n−1θnµ(·) +
Λn (X1 , ..., Xn )(·).
θ+n
θ+n

∑n

i=1 δXi is the empirical distribution.

In Bayesian context, F iscalled a prior process, and the posterior distribution is called the prediction rule. The sequential
sampling scheme is well known as the Blackwell-MacQueen urn scheme (1973), or the Chinese
restaurant process. The prediction rule induces measures on partitions. Let the j-th firstly
appear value of (X1 , ..., Xn ) be Xj∗ , j ∈ {1, 2, ..., k}. Then, Pn,k ∋ (n1 , ..., nk ), nj := #{i; Xi =Xj∗ } is a integer partition of a positive integer n. The distribution of the multiplicities of integers(c1 , ..., cn ), ci := #{j; nj = i}, isn! ∏P(C1 = c1 , ..., Cn = cn ) =
(θ)n
n

i=1

( )ci
1
θ
,
i
ci !

where (θ)n := θ(θ + 1) · · · (θ + n − 1). This measure on partitions is called the Ewens sampling
formula (1972; Antoniak 1974, Sibuya 1993). A generalization, (exchangeable) Gibbs partitions,are commonly used to characterize prior processes.
Definition 1 (Pitman 2006; M 2018). Gibbs partition is the probability measure on partitions
λ ⊢ n ∈ N := {1, 2, ...} of the form
P(C = c) =n (∏vn,l(c)w i )c i 1n!
.
Bn (v, w)
i!
ci !
i=1

Here, l(c) = c1 + · · · + cn is the length, and the normalization constant is written as
Bn (v, w) =n∑vn,k Bn,k (w),Bn,k (w) =∑c∈Pn,k

k=1

n!

n (
∏
wi )ci 1
i=1

i!

ci !

.

The latter is known as the partial Bell polynomial.
Examples of Gibbs partitions are presented in Proceedings and in M (2018). The characterizations are given by Gnedin & Pitman (2005). A Gibbs partition is infinitely exchangeable iﬀ
1
2

Based on ongoing work with Jaeyong Lee at Seoul National University.
E-mail: smano@ism.ac.jp-5-wi = (1 − α)i−1 , α < 1. Without exchangeability, the prediction rule is not available anymore(Pitman 1995, Lee et al. 2013, Proceedings) and we need another sampling scheme than theBlackwell-MacQueen urn. For a Gibbs partition, the length l(c) = k is the suﬃcient statistic for parameters (vn,k ) and the conditional distribution is an A-hypergeometric distribution(Takayama et al. 2018):
P(C = c|AC = b) =xc
1,
ZA (b; x) c!

c ∈ Nn−k+1,0

where the A-hypergeometric polynomial ZA (b; x) is defined by()()
0 1 2 ··· n − k
n−k
wi
A=
, b=
, xi =,i!
1 1 1 ···
1k
and n!ZA (b; x) = Bn,k (w). M (2017) obtained a direct sampler for A-hypergeometric distributions, and it provides a direct sampler from a Gibbs partition without exchangeability. The
sampler can be used as an alternative to the Blackwell-MacQueen urn scheme.
A typical Bayesian semiparametric setting is data (Y1 , ..., Yn ) derived from a hierarchicalmodel (MacEachern 1994; Escober & West 1995)
ind

Yi |Xi , σ ∼ N(Yi ; Xi , σ), σ ∼ π(σ),

ind

Xi |P ∼ P, P ∼ DP(θ; µ).

The marginal with the Dirichlet process, (X1 , ..., Xn ) ∼ P(X1 , ..., Xn ), is sampled directly by
the Blackwell-MacQueen urn scheme. Let us call the sampling from the marginal with a priorprocess prior sampling. To draw from posterior π(Y, σ|X) (we will call posterior sampling), a
Gibbs sampler is used, where we iterative draw values from conditional distributions
Xi |(X−i , σ, Y ),i ∈ [n],σ|(X, Y ).

The prediction rule gives the updating ruleP(Xi = ·|X−i , σ, Y ) ∝ N(Yi ; Xi = ·, σ)θµ(·) +k∑

N(Yi ; Xj ∗ , σ)nj δXj ∗ (·).

(1)

j=1

Let us consider use of a non-exchangeable Gibbs partition. Thanks to the direct sampler
from a Gibbs partition introduced above, the prior sampling is straightforward. For posterior
sampling, a possibility is use of an independent Metropolis–Hastings algorithm with the priorsampling. This cause poor acceptance, since the acceptance is determined based on the vector(X1 , ..., Xn ). Improvement of the acceptance with aid of n-exchangeability was discussed in the
talk.
Remark 1. After the talk, we found that mixture modeling with non-exchangeable Gibbs
partitions have been extensively discussed, where a non-exchangeable Gibbs partition with vn,k =1 is called a product partition model (Quintana & Iglesias 2003, JRSS B 65: 557-574). It seems
that a Gibbs sampler with updating rule similar to (1) has been used. The updating rule is notconsistent with infinite exchangeability.

-6-統計的形状モデルを利用した大規模点群レジストレーションとその高速化
広瀬 修†
金沢大学理工研究域†

1. はじめに
点群位置合わせ問題は，物体の形状を表現する
２つの点集合間に対し，１つの点群をもう一方の点群に移す写像を求める問題である．点群位置合わせ問題は想定する写像の種類に応じて剛体変換と非剛体変換のものに分類され，最近では非剛体変換に基づいた点群位置合わせが，その柔軟さのため非常に活発に研究されている．
Coherent point drift (CPD) は非剛体変換に基づいた点群レジストレーション手法の代表的な手法である [Myronenko 2010]．CPD の成功の要因としてまず挙げられるのが，外れ値への耐性である．ここで，外れ値とは点群によって表現される形状とは無関係に存在する点とする．
CPD は点群位置合わせ問題を混合確率分布の推定問題として定義する．その際，混合分布の構成
分布の１つとして外れ値の分布を明示的に与えることが，外れ値への耐性の主要な要因である．
CPD のもう１つの成功要因として挙げられるのが，非 剛 体 換される点群に対する「変位場の滑らかさ」である．変位場の滑らかさとは，非剛体変換される点群を構成する任意の点の変位と，その他の点の変位が，その距離が近ければ近いほど相関するとした仮定である．この仮定は非常に自然な仮定であるため，CPD は多くの点群位置合わせ問題において精度の高い位置合わせ結果を与える．一方で，変位場の滑らかさの仮定が適切ではない場合，CPD は容易に位置合わせに失敗する．例えば，人間の手の形状マッチングを行う場合，人差し指と中指を構成する点は比較的近くに位置するが，その動きは逆相関する傾向があり，このような場合には変位場の滑らかさの仮定だけでは不十分であるためである．
この問題を克服する方法の１つとして挙げられるのが，教師あり学習に基づく方法である．
もし人差し指と中指の動きが逆相関する傾向にあることを事前に知っていれば，その知識を位置合わせアルゴリズムに組み込むことにより,高精度の位置合わせが期待できるからである．
今回, 新たに開発した教師あり学習法に基づいた点群位置合わせ手法に対する位置合わせ性能-7-の評価について報告する．開発手法は外れ値の分布を構成分布の１つとして持つ混合分布に基づいているため，外れ値への優れた耐性を有す
る．また，物体の形状変化モデルに，訓練デー
タから得られる事前知識を組み込むため，変位場の滑らかさの問題を自然に解決することができる．講演では大規模な点群データに対応するための高速化についても議論を行う．

2. 実験
開発した位置合わせ手法と代表的な点群位置合
わ せ 手 法 で あ る CPD と Thin Plate Spline
Robust Point Matching (TPS-RPM) [Chui 03]の性能の比較を行った．使用したデータは IMM
hand データである[Stegmann 2002]．このデー
タは４０種類の人間の手の画像に対し，手と背景の境界部分に５６個の特徴点を人手で打点したものである．各々の特徴点は４０種類の画像で対応関係がとられるように打点されている．
図１が IMM hand データ中の点群番号６に対するCPD, TPS-RPM および開発手法の位置合わせ結果を表す. 図１の１段目が推定されるべき正解の手の形状を表し，２段目が位置合わせの対象となるデータで，手法の頑健性を検証するため正解データに４種類の改変を施したものである．
左から順に (1)点の複製，(2)欠損, (3)外れ値の付加，(4)回転の改変を表す．図１の３段目の赤色で示された点が最適化の初期形状，すなわ
ち，平均形状を表す．３〜５段目がそれぞれ提案手法，CPD，TPS-RPM の適用結果を表す．デー
タ (1)，(2)，(3)，(4)の全てで提案手法が正解と ほ ぼ 同 一 の 形 状 の 推 定 に 成 功 し た ． CPD と
TPS-RPM については，全てのデータで親指以外の指が細くなる現象が見られた．これは，変位場の滑らかさの欠点の１つを表している．ある
指を形成する点が，近隣に存在する他の指を形成する点ともその変位が相関することが原因である. データ(1) に対しては全ての手法で概ね正解の形状を空いてしたが，CPD と TPS-RPM には
指が細くなる現象が見られた. データ(2) に対して提案手法は欠損領域の点群の推定に成功した．一方で CPD と TPS-RPM は親指を含む全ての指が短く描画された．

図 1 ． 提 案 手 法 ， CPD ， TPS-RPM の 比 較 ． １ 段 目 の 図 は 推 定 さ れ る べ き 正 解 の 形 状 を 表 す ． ２
段目の図は位置合わせの対象となるデータ．３段目の赤の点群が最適化の初期形状を表す．手法 の 頑 健 性 を 検 証 す る た め 正 解 デ ー タ に ４ 種 類 の 改 変 を 施 し た ． 左 か ら 順 に (1) 点 の 複 製 ，(2) 欠 損 ， (3) 外 れ 値 の 付 加 ， (4) 回 転 の 改 変 を 表 す ． ３ 〜 ５ 段目がそれぞれ提 案 手 法 ， CPD ，
TPS-RPM の 適 用 結 果 を 表 す .

-8-

Kolmogorov-Smirnov test based onkernel estimation
Rizky Reza Fauzi, Graduate School of Mathematics Kyushu University
Maesono Yoshihiko, Faculty of Mathematics Kyushu University
1. Boundary-free kernel distribution function estimators
Let X1 , X2 , ..., Xn be independently and identically distributed random variables with
an absolutely continuous distribution function FX and a density fX . The classicalnonparametricestimator of FX has been the empirical distribution function Fn (x) =∑ni=1 I(Xi ≤ x)/n(I(·): indicator function). Let K be a kernel function and h > 0
be a bandwidth.
∫ ∞We assume that K is a symmetric (about 0) continuous nonnegative
function with −∞ K(v)dv = 1, as well as h → 0 and nh → ∞ when n → ∞. Since
distribution function is actually an integral of density function, this kernel density estimator gave an idea to define a kernel distribution function estimator. Nadaraya (1964)
defined it as()
nx − Xi1∑b
FX (x) =W, x ∈ R,n i=1h∫v
where W (v) = −∞ K(w)dw.
Under the condition that fX has one continuous derivative fX′ , it has been proved
that, as n → ∞,
∫ ∞
h2 ′
b
Bias[FX (x)] =f (x)v 2 K(v)dv + o(h2 ),
2 X
−∞
( )
2h
h
1
b
V ar[FX (x)] =FX (x)[1 − FX (x)] − r1 fX (x) + onnn∫∞where r1 = −∞ vK(v)W (v)dv. It is easy to show that r1 is a nonnegative number.
If we deal with R+ or unit interval, the standard kernel distribution function estimator will suﬀer the so called boundary bias problem. To solve this problem, we utilise afunction g which bijectively transform the support A of the random variable under consideration into R, then doing the usual standard kernel distribution function estimationof Y = g(X). Here we propose an estimator[]ng(x) − g(Xi )1∑eW
FX (x) =, x ∈ A.
n i=1h
Just as the standard one, we have a bias and a variance∫ ∞h2eBias[FX (x)] =c(x)v 2 K(v)dv + o(h2 )
2
−∞
( )
f
(x)
1
2h
h
X
V ar[FeX (x)] =FX (x)[1 − FX (x)] −r1 + o,′nn g (x)n-9-

where r1 =∫∞−∞vK(v)W (v)dv andc(x) =fX′ (x)fX (x)g ′′ (x)−.
[g ′ (x)]2[g ′ (x)]3

2.Boundary-free smoothed Kolmogorov-Smirnov type test
The Kolmogorov-Smirnov (KS) test is the most popular GOF test used in practice.
Unfortunately, it lacks of smoothness that can lead to smaller power at the tails, which
is important in many practical applications. It is natural if one uses the naive kernel
distribution function estimator in place of the empirical distribution function. Thus,
instead of the standard KS statistic
Dn =

sup

−∞<x<∞|Fn (x) − FX (x)|being used to test whether random variable X having FX as its distribution, we canreformulate by smoothing it to
b=
D

sup

−∞<x<∞|FbX (x) − FX (x)|,where Fb is the kernel distribution function estimator.
However, a new problem is raising when the support of the random variable that weare dealing with is not the entire real line, i.e. boundary problem. To overcome thisproblem, we propose to use our estimator in section 1 to substitute empirical distributionfunction in standard KS statistic. Therefore, we define the boundary-free smoothed KS
type test as
e = sup |FeX (x) − FX (x)|,D−∞<x<∞

where FeX is our boundary-free kernel distribution function estimator. We compare anumerical study by calculating simulated power of our proposed test with n = 50, andthen we compare it with the result of the standard KS test.
The asymptotic behaviours of our proposed test statistic are stated in the following
theorems.
Theorem 1 Let X be a random variable with distribution function FX supported on aset A. If FeX is the proposed boundary-free kernel distribution function estimator, then
)
(
1
e
.
sup |FX (x) − FX (x)| = op √n−∞<x<∞
Theorem 2 Let X be a random variable with distribution function FX supported on ae is the proposed boundary-free smoothed KS-type statistic, thenset A. If D
√
]
[
∞
√
2π ∑
−(2i − 1)2 π 2elim Pr( nD ≤ x) =
exp
.
n→∞
x i=1
8x2

-10-

Tests for high-dimensional covariance matrices and correlation matrices
under the strongly spiked eigenvalue model
Aki Ishiia , Kazuyoshi Yatab and Makoto Aoshimaba1

Department of Information Sciences, Tokyo University of SciencebInstitute of Mathematics, University of Tsukuba

Introduction

we consider the equality test of covariance matrices when the data dimension is much larger than the sample
size. Suppose we have two classes πi , i = 1, 2. We define independent d × ni data matrices, X i =[xi1 , ..., xini ], i = 1, 2, for πi , i = 1, 2. We assume that xij , j = 1, ..., ni , are independent and identically
distributed (i.i.d.) as a d-dimensional distribution with a mean vector µi and covariance matrix Σi . We
assume ni ≥ 4, i = 1, 2. The eigen-decomposition of Σi is given by Σi = H i Λi H Ti , where Λi =diag(λ1(i) , ..., λd(i) ) having λ1(i) ≥ · · · ≥ λd(i) (≥ 0) and H i = [h1(i) , ..., hd(i) ] is an orthogonal matrix ofthe corresponding eigenvectors. We assume λ2(i) > 0 for i = 1, 2, and λ1(i) s are of multiplicity one in thesense that lim inf d→∞ λ1(i) /λ2(i) > 1 for i = 1, 2. We consider the equality test of covariance matrices asfollows:
H0 : Σ1 = Σ2 vs. H1 : Σ1 ̸= Σ2 .
(1.1)
Aoshima and Yata [1] gave a test procedure based on the quantity of tr(Σ1 −Σ2 ). They also discussed sample
size determination so as to have a prespecified size and power simultaneously. Li and Chen [5] considered thetest problem by using the quantity of tr{(Σ1 − Σ2 )2 }. The above references discussed asymptotic propertiesof their test procedures when d → ∞ and ni → ∞ under the following eigenvalue condition:
λ21(i)tr(Σ2i )→ 0 as d → ∞ for i = 1, 2.

(1.2)

Aoshima and Yata [2] called (1.2) the “non-strongly spiked eigenvalue (NSSE) model”. On the other hand,
Ishii, Yata and Aoshima [3] investigated asymptotic properties of the first principal component and considered the test problem (1.1) when d → ∞ while ni s are fixed under the following eigenvalue condition:∑d2
s=2 λs(i)(1.3)= o(1) as d → ∞ for i = 1, 2.
λ21(i)
Note that (1.3) implies the conditions that λ2(i) /λ1(i) → 0 and λ21(i) /tr(Σ2i ) → 1 as d → ∞. The condition(1.3) is generalized as{ λ2(A-ii)

lim infd→∞1(i)tr(Σ2i )}> 0 for i = 1 and 2.

Aoshima and Yata [2] called (A-ii) the “strongly spiked eigenvalue (SSE) model” and showed that highdimensional data often have the SSE model.

2

A new test procedure under the SSE modelLet ∆ = ||Σ1 − Σ2 ||2F = tr{(Σ1 − Σ2 )2 }. Under (1.3), it holds that as d → ∞∆ = (λ1(1) − λ1(2) )2 + 2λ1(1) λ1(2) {1 − (hT1(1) h1(2) )2 } + o(λ21(1) + λ21(2) ).

(2.1)

We consider (2.1) as a starting point to handle the SSE model (A-ii). We give a test statistic based on (2.1)
and show that it holds an asymptotic null distribution even when (1.3) is not met. Let X i = [x̄i , ..., x̄i ] and

-11-

∑ni
−1
x̄i = n−1
j=1 xij for i = 1, 2. We define the sample covariance matrix as S ini = (ni − 1) (X i −i
X i )(X i − X i )T for i = 1, 2. We denote the dual matrix of S ini by S iD and define its eigen-decomposition
as follows:n∑i −1
S iD = (ni − 1)−1 (X i − X i )T (X i − X i ) =λ̂s(i) ûs(i) ûTs(i) ,
(2.2)
s=1

where λ̂1(i) ≥ · · · ≥ λ̂ni −1(i) (≥ 0) and ûs(i) denotes a unit eigenvector corresponding to the eigenvalue
λ̂s(i) . Note that S ini and S iD share non-zero eigenvalues. If one uses the noise reduction (NR) method given
by Yata and Aoshima [7], λ1(i) and h1(i) are estimated byλ̃1(i) = λ̂1(i) −tr(S iD ) − λ̂1(i)and h̃1(i) = {(ni − 1)λ̃1(i) }−1/2 (X i − X i )û1(i) .
ni − 2(2.3)

for i = 1, 2, where û1(i) is given in (2.2). Let δi = tr(Σ2i ) − λ21(i) , i = 1, 2. If one applies the NR method
to estimating ∆ in (2.1) straightforwardly, the estimation of ∆ includes huge noise for high-dimensional
1/2
data. In this talk, we evaluate the amount of the noise and consider bias correction. Let η = δ1 /λ1(1) +1/2

δ2 /λ1(2) . By using the cross-data-matrix (CDM) method given by Yata and Aoshima [6] and the NR
method, we estimate η by η̂. We provide the following new test statistic:T(λ̃1(1) − λ̃1(2) )2 + 2λ̃1(1) λ̃1(2) {1 − min{1, (h̃1(1) h̃1(2) )2 }}1+η̂
TNR =
.
∑2
2
i=1 2λ̃1(i) /(ni − 1)
Proposition 2.1 (Ishii, Yata and Aoshima [4]). Under (A-ii), H0 and some regularity conditions, it holds
that TNR ⇒ χ21 as min(n1 , n2 , d) → ∞.
We gave another test statistic under (A-ii) by dividing high-dimensional eigenspaces into the first eigenspace
and the others. We gave asymptotic null distribution and the power of the test statistic. We demonstrated
the new test procedure by using actual microarray data sets. We also considered a correlation test for highdimensional data under the SSE model.

References[1] Aoshima, M., Yata, K., 2011. Two-stage procedures for high-dimensional data. Sequential Anal. (Editor’s special invited paper) 30, 356-399.
[2] Aoshima, M., Yata, K., 2018. Two-sample tests for high-dimension, strongly spiked eigenvalue models.
Stat. Sin. 28, 43-62.
[3] Ishii, A., Yata, K., Aoshima, M., 2016. Asymptotic properties of the first principal component andequality tests of covariance matrices in high-dimension, low-sample-size context. J. Stat. Plan. Inference 170, 186-199.
[4] Ishii, A., Yata, K., Aoshima, M., 2018. Equality tests of high-dimensional covariance matrices under
the strongly spiked eigenvalue model. J. Stat. Plan. Inference, revised.
[5] Li, J., Chen, S.X., 2012. Two sample tests for high-dimensional covariance matrices. Ann. Statist. 40,
908-940.
[6] Yata, K., Aoshima, M., 2010. Effective PCA for high-dimension, low-sample-size data with singularvalue decomposition of cross data matrix. J. Multivariate Anal. 101, 2060-2077.
[7] Yata, K. and Aoshima, M. (2012). Effective PCA for high-dimension, low-sample-size data with noise
reduction via geometric representations. J. Multivariate Anal. 105, 193-215.

-12-

癌の高次元遺伝子解析の諸問題（1）
成蹊大学
新村秀一

名誉教授

1. 癌の高次元遺伝子解析を取り巻く諸問題について
1999 年から 2004 年まで、米国の 6 研究グループが医学誌に論文を掲載し、研究に用いた microarray データをインターネット上に公開した。ハーバード大学医学部の Golub ら（1999）はサイエンスに論文を発表し、「顕微鏡などで細胞や遺伝子の生物学的な変異から癌遺伝子を見つける限界に言及し、1970 年から
microarray による体系的な研究が必要であり研究しているが、大きな成果を未だ得ていない」と真摯に述べている。彼らは独自の手法を開発し、それらを Self-organizing maps(SOM)や LOO や生存時間解析等の統計手法と組み合わせている。これらのデータは、例えば癌と健常者の 2 群 100 例を、1 万個の遺伝子の発現量で判別する問題であり、最も判別に適した遺伝子をがん遺伝子と特定する研究である。しかし、幾つかの論文では SVM の利用は述べられているが明確な結果の紹介はない。そして残念なことに NIH が乳がんを除いて、この種の研究に疑問を呈したレポートを出し、医学的な研究が終わったようだ。
2．新しい判別理論と Matryoshka Feature Selection Method(新理論 2)による癌の遺伝子解析の成功
筆者は 2015 年 10 月 25 日に富山市で開催された科研費シンポジュームで石井博士の発表で上記の 6 種のデータが公開されていることを知った。彼女から 28 日にメールで microarray を入手できる HP を知り入手した。そして 12 月 22 日までの 56 日間で癌の遺伝子解析を完成させた。
Fact3：6 種の microarray すべてが LSD すなわち最小誤分類数 MNM=0 である。
Fact4：全ての microarray は、多数の線形分離可能な Small Matryoshuka(SM)とそれ以外の雑音空間に分割できる。最近まで、2 群が各 SM の部分空間で完全に分かれて識別できるので、SM に含まれる遺伝子が癌遺伝子であり、癌の Basic Gene Set(BGS)と考えていた。そしてその和集合と各 BGS は信号空間であると定義した。すなわち、
1） 分散共分散行列に基づく判別分析は LSD を正しく判別できないので、癌の遺伝子解析に全く役に立たない。
2） 改定 IP-OLDF（RIP）と改定 LP-OLDF とハードマージン最大化 SVM（H-SVM）だけが microarrays が LSDであることがわかる。しかしなぜ多くの研究者は SVM を利用しているのに、この重要な信号を見過ごしたかという問題が起きる。
3） RIP と改定 LP-OLDF だけが、microarray を多数の SM と雑音空間に分割できる。これが今回の発表のテーマである。
そして、SM は小標本であるので、簡単に分析できる。しかし、2 群が各 SM で完全に分かれているにもかかわらず、主成分分析、クラスター分析、一元配置による分散分析などで線形分離可能な事実を示さない。
ロジスティック回帰は全ての SM が NM=0 であるが、分散共分散行列による LDF や QDF は NM=0 にならないものが多い。一方、RIP、改定 LP-OLDF と H-SVM は、2 つの SV で-1 以下に class1、1 以上に class2 の症例を正しく判別し、判別スコア(Discriminant Score, DS)の範囲に対する比率 RatioSV が 5％以上になるものが多い。そこで遺伝子の変わりに各 SM に含まれる遺伝子の総合特性値である RIP、改定 LP-OLDF と H-SVMの DS である RipDSs、LpDSs と HsvmDSs を変数とするデータを作成した。これらを PCA やクラスター分析で分析すると、2 群は完全に分かれる。
以上から SM で、癌と健常者の 2 群は完全に分かれているが、通常の統計手法ではそれがわからない。しかし、DS で作成したデータでは LSD の事実が簡単にわかる。すなわち RIP、改定 LP-OLDF と H-SVM の判別
スコアが、高次元 microarray の信号でないかと考えるに至った。これに関する詳細は、広島のシンポジュームで報告を予定している。
3．統計的判別分析の問題
Golub 以前にも他の遺伝子データが公開されていて、統計研究者やパターン認識などの工学者研究者は、質が高く無償で提供されたこの高次元データ(small n large p)を、格好の研究テーマとしてとらえて研究を行ってきた。多くの研究論文には”Feature Selection Methods（統計的にはモデル選択とか変数選択）”
とか”Filtering”という用語を含んだものが多い。そして医学論文には見られない、次の 3 つの困難を指摘する論文もある。

-13-

1） 高次元データ(small n large p；n<<p)の困難さ：これは、例えば 100 症例から 1 万次元の分散共分散行列を構築することが端的な事例である。2000 年以前に、国際会議で発表もあったがいつのまにかなくなった。2015 年 11 月に六本木で行われた JMP の Discovery Summit で JMP の創業者の Sall 博士が特異値分解を用いた横長データの Fisher の LDF を発表した。無償で 1 か月借り受け、6 種の microarrayを分析したが、誤分類数（NM）は高かった。すなわち、分散共分散に基づく LDF は LSD を正しく判別できず、癌の遺伝子解析に役に立たない点である。また、数理計画法では全く問題にならずむしろ largen small p の方が、制約式が増えて計算時間がかかる。すなわち、この困難は 2 変数の相関でもって 1
万変数が関係づけられる統計に限定された困難である。数理計画法は変数間の関係が小さい。
2） NP-hard:1 万変数の判別分析で、適切な部分モデルを選ぶことは困難である。この困難の真の意味を考えていない側面がある。即ち統計的判別関数や 2 次計画法で定式化された H-SVM は、定義域で唯一の最適な判別関数が求まる。部分空間にも最適解がある場合、モデル選択で部分モデルの中から最適解を見つける必要がある。さらに問題なのは LSD では、最小次元の BGS を含むすべてのモデルが最適解になることである。Feature Selection で MNM=0 という基準で最適化モデルを探さない限り、これ等の研究は無意味である。
3） 信号と雑音の分離：この問題は的を得ているが、「信号」の定義がはっきりしない。癌の Microarray データにおける信号は MNM=0 である。2017 年 1 月の金沢における科研費シンポジュームで、私が初めて癌の遺伝子解析に成功したという説明に、青嶋氏より我々の方が先行している旨の意見があった。よく考えてみると、2015 年の富山で石井氏が microarray データで PCA を行うと、第 1 固有値だけがスパイ
ク上に大きな固有値になり、一般的な常識で考えられない結果になるという話を聞いた。しかし、私は
microarray が簡単に入手でき、これまで判別分析で 4 つの問題を解決してきたが、「未解決の 5 番目の癌の遺伝子解析」を解決していないことに気づいた。癌と健常あるいは異なった 2 種の癌が、遺伝子空間で完全に分かれて 2 つの球に布置しているというのが青嶋と矢田らの結論である。このことは、私の「2 群は microarray で MNM=0 であり、それが多くの MNM=0 である SM と MNM が 1 以上の雑音空間に分割できる」という驚く結果を統計的に検証した研究であることに気づいた。
4．癌の遺伝子解析から癌の遺伝子診断
癌の遺伝子解析は、判別分析の 4 つの問題と応用問題として 6 種の microarray で全ての SM を求めることができ Springer の本で解決した問題 5(Shinmura, 2016d)を指す。豊富な実証研究の成果である。
2016 年になって、SM は n 個以下の遺伝子で構成された小標本であるので、統計手法で簡単に分析できると考えた。しかしロジスティック回帰だけが NM=0 になり、PCA やクラスター分析では２つの class が線形分離可能な事実を示さなかった。それも仕方がないと当初は考えたが、RIP、改定 LP-OLDF と H-SVM は、MNM
あるいは NM が 0 である。そこで判別スコアの範囲に対して 2 つの SV 間の距離 2 の比を RatioSV として求めると 6 種の最大値は[11.67%, 38.93%]と異常に大きい。これに反して、Alon の 130 個の BGS は 1％未満である。そして、RIP の判別スコアを遺伝子の代わりに変数として用いたデータを作成した（RipDSs デー
タ）。これを PCA で分析すると健常症例が第 1 主成分で負のある値以下に、癌症例が正のある値以上に布置することが分かった。そして各 SM で求めた RipDS と PCA で求めた総合化された RipDS を癌の悪性度指標と呼び、医学の素人ながら癌の遺伝子診断の突破口を開いたのではないかと考えた。それらの成果をまとめて、2017 年に Amazon から Kindle 版として出版した。予約注文で 600 部以上がダウンロードされ幸先が良いと思ったがその後が続かない。きっと Research Gate で出版案内したので、癌の遺伝子解析に関連した研究者が NIH の敗北宣言後も 600 人程度は細々といると考えられる。9 月時点で RG の Read 数が 11 万を超えたが癌の遺伝子解析関連の Read 数は少ない。また 7 月末にラスベガスで開催された Biocomp18 で 8 月 3日（金）に開催ホテルの Luxor で朝 3 時に目が覚め空港に行く 7 時まで,初めて約 140 にいた Following とFollower の所属を調べた。帰国後 RG がダウンし復旧後に 1391 人に増えた。ひょっとして Biocomp18 で遺伝子関連の専門家に注目されたかと調べてみたが、40 人程度しかいないようである。また、彼らが癌の遺伝子関連の Draft を特に読んでいる事実も得られなかった。癌の遺伝子診断の成果は、専門家に検証してもらわなければ意味がないので、袋小路に入っている。
今後の課題として、青嶋・矢田らの結果に対して、筆者の研究アプローチで具体的な幾つかの事実で同じことを示していることを示す予定でいる。
[1] Aoshima M, Yata K (2017) Two-sample tests for high-dimension, strongly spiked eigenvalue models, Statistica Sinica
Preprint No.ss-2016-0063R2:1-31-14-

A Robust-filtering Method for Small Sample
Economic Time Series
Naoto Kunitomo
In collaboration with Seisho Sato 1

October 25, 2018
Key Words
Non-stationary economic time series, Errors-variables models, Measurement Error,
trend and seasonality, Robust-filtering, SIML, Fourier-Inversion, Large Dimension.

Summary
We develop a new filtering method to estimate the hidden states of random variables and to handle multiple time series data, and small sample economic time series
in particular. Kunitomo and Sato (2017), and Kunitomo, Sato and Kurisu (2018)
have developed the separating information maximum likelihood (SIML) method for
estimating the non-stationary errors-in-variables models. They have discussed theasymptotic properties and finite sample properties of the estimation of unknown
parameters. We utilize their results to solve the filtering problem of hidden randomvariables, which gives a powerful new method of handling macro-economic timeseries.
Kitagawa (2010) has discussed the standard statistical filtering methods alreadyknown including the Kalman-filtering and the particle-filtering methods. Since (i)these methods depend on the underlying distributions such as the Gaussian distributions for the Kalman-filtering and (ii) the procedures essentially depend on thedimension of state variables, there may be some diﬃculty to extend to the highdimension cases even when it is fixed, say 100. On the other hand , we can expectthat our method has some merits when we handle small sample economic times series with non-stationarity and seasonality with many variables because our method
does not depend on the specific distributions as well as the dimensions of randomvariables. See Kunitomo, Awaya and Kurisu (2017) for a comparison of small sample
properties of the ML and SIML methods.
The most important feature of the present procedure is that it can be applicableto small sample time series data with large dimension. Also our new method has asolid mathematical and statistical foundation, which is related to Doob (1953) and
Hannan (1971).
1

We owe Akihiko Takahashi for a suggestion on our filtering method.

-15-

References[1] Brockwell, P. and R. Davis (1990), Time Series : Theory and Methods, 2nd
Edition, Wiley.
[2] Doob, J.L. (1953), Stochastic Processes, Wiley.
[3] Hannan, E.J. (1970), Multiple Time Series, Wiley.
[4] Kitagawa, G. (2010), Introduction to Time Series Analysis, CRC Press.
[5] Kunitomo, N. (2018), ”Common Seasonal Factors,” in preparation.
[6] Kunitomo, N. , Sato and D. Kurisu (2018), Separating Information Maximum
Likelihood Estimation for High Frequency Financial Data, Springer.
[7] Kunitomo, N. and S. Sato (2017), ”Trend, Seasonality and Economic Time
Series : the Non-stationary Errors-in-variables Models,” SDS-4, MIMS, Meiji
University, http://www.mims.meiji.ac.jp/publications/2017-ds.
[8] Kunitomo, N., N. Awaya and D. Kurisu (2017), ”Some Properties of Estimation Methods for Structural Relationships in Non-stationary Errors-in-Variables
Models,” SDS-3, MIMS, Meiji University.

-16-

Applications of Distance Correlation to Time
Series : 距離の相関係数の時系列解析への応用
南山大学

松井 宗也

研究内容
近年、「距離の相関係数」(Distance Correlation) という 2 つの確率ベクトルの依存関係を測る指標が注目されている。ここで 2 つのベクトルの次元は任意で異なっていてもよい。2005 年に Székely, G.j. により提案されて以来盛んに研究されている（詳しい定義等は [Székely et al.(2007)] や [Székely and Rizzo(2009)] 参照）。注目すべき特徴は、ケンドールやスピアマンの順位相関係数といった従来の指標では捉えることの難しい非線形な相関も検知できることである。そもそも 2 つの確率変数が独立であるとは、その同時特性関数がそれぞれの周辺特性関数の積として書けることと同値である。「距離の相関係数」の基本的なアイデアは、同時特性関数と独立な場合のそれ（各々の周辺関数の積）の距離をみて相関を測るものである。2 つの確率ベクトル X ∈ Rp と Y ∈ Rq（次元は任意で異なっていてもよい）の同時特性関数
を φX,Y (s, t) とおき、またそれぞれの周辺特性関数を φX (s)、φY (s) とおく（ここでp, q ∈ N）。また重み付き測度（重み関数（正）かけるルベーグ測度等）を µ(s, t) とおく。すると「距離の相関係数」は∫|φX,Y (s, t) − φX (s)φY (t)|2 µ(ds, dt)
Rp+qと重み付けした L2 距離で定義される。この形を見れば統計量は独立ならばその時に限り 0 となることがわかる。特に Székely,G.j. らは具体的に重み関数を µ(ds, dt) =cp,q |s|−α−p |t|−α−q ds dt とおくことで統計量をT (X, Y ; µ) = E[|X − X ′ |α |Y − Y ′ |α ] + E[|X − X ′ |α ] E[Y − Y ′ |α ]−2 E[|X − X ′ |α |Y − Y ′′ |α ]と陽な形に表した。ここで (X ′ , Y ′ ) は (X, Y ) の iid コピーで、Y ′′ は全ての変数と独立な Y のコピーである。
本研究はこの「距離の相関係数」を 1 次元と多次元の定常な時系列へ応用するものである。主結果は論文 [Davis et al.(2018)] にまとめてある。主な目的は
１．2 つの時系列が独立かどうかを検定する、あるいは「距離の相関」で 2 つの時系
列の相互依存関係を測る。
２．タイムラグをとった系列の「距離の相関」をみることで、１つの時系列が系列相関を持つかどうか検定する。
ことである。時系列解析では、多次元の系列相関を測る指標として（多次元の）自己相関関数が一般的に知られている。ここで言う相関はピアソンの積率相関係数である。非常に便利な指標であるが、非線形時系列モデルの依存関係を捉えるのが難-17-

しいことが知られている。非線形時系列は計量ファイナンス分野において近年よく研究されていて、例えばその代表例として GARCH(p, q) モデルが挙げられる。このモデルは、単に原系列の自己相関をみても依存関係を検出できず、原系列を変換したも（例えば絶対値や 2 乗を考えたもの）の自己相関をみて初めて依存関係が捉えられる。
本研究のアイデアは、そのピアソン流の相関を「距離の相関係数」に置き換える
ことで依存関係をより良く検出しようというものである。論文では新しく「自己距離の相関関数」という指標を自己相関関数に代わりうるものとして提案した。そして、ミキシング条件 (strong mixing) のもとで、標本「自己距離の相関関数」に関連する漸近論を導出した。まず、緩やかモーメント条件のもとで統計量の一致性を示した。その後いくつか追加的な条件を与え、依存関係がある場合と独立な場合の両方で統計量の漸近分布を理論的に導出した。独立な場合とそうでない場合は収束のオーダーが異なることも示した。検定統計量は特に独立な場合は、独立なカイ２乗
変数の線形無限和で表され、一般に上側確率を求めるのが難しいことが知られている。これに対し、ブートスラップ法を用いた数値実験で漸近分布がうまく近似できることも示した。
２つ目の研究である 1 系列における「自己距離の相関関数」に関しては、具体的なモデルとして AR(p) モデルを考えた。モデルの適合度をみるために、パラメータ推定後に得られる残差に「自己距離の相関関数」を適応した。すると興味深いことに、iid 系列に「自己距離の相関関数」を適応した場合と比較して明らかに異なることもわかった（これは自己相関関数と同様の現象である）。その他、誤差項に裾の重い分布を仮定した場合は Székely, G.j. のオリジナルな定義では漸近論がうまく機能
しないこともわかった。パラメータの推定も考えたもとでは、モーメント制約がより厳しくなるためと考えられる。我々が提案した他の重み関数によるものは、裾の厚い誤差項に関してもうまく対応することも確認できた。
実証研究も行い、アマゾンの株価収益率データや風速データなどで提案した方法の応用を試みた。
キーワード：Auto- and cross-distance correlation function, testing independence,strong mixing, ergodicity, Fourier analysis, U -statistics, AR process, residuals

参考文献[Davis et al.(2018)] Davis, R.A., Matsui, M., Mikosch, T. and Wan, P. (2018) Applications of distance correlation to time series. Bernoulli 24, 3087–3116.
[Matsui et al.(2017)] Matsui, M., Mikosch, T. and Samorodnitsky G. (2017) Distance covariance for stochastic processes. Probab. Math. Statist. 37, 355–372.
[Székely et al.(2007)] Székely, G.J., Rizzo, M.L. and Bakirov, N.K. (2007) Measuring and testing dependence by correlation of distances. Ann. Statist. 35, 2769–2794.
[Székely and Rizzo(2009)] Székely, G.J. and Rizzo, M.L. (2009) Brownian distancecovariance. Ann. Appl. Stat. 3, 1236–1265.

-18-割引因子ポアソン-ガンマ状態空間モデルによる計数時系列データの逐次分析
入江 薫
科研費シンポジウム「融合する統計科学」において、二日目にあたる十二月一日の午前のセッションにて、「割引因子ポアソン-ガンマ状態空間モデルによる計数時系列データの逐次分析」という演題で研究成果を発表した。
発表は予稿の内容に沿って行い、公表済みの論文（Chen et al., 2018）および進行中のプ
ロジェクト（Irie et al., 2018）の二つの研究について発表した。シンポジウムの主旨を踏まえ、提案するモデルや計算手法の詳細は省き、主に使用しているデータについて詳しく紹介した。具体的にはニュースウェブサイトへのアクセス数のデータであるが、(1)これは非負の整数値を取る多変量の時系列データでありモデリングに工夫を要すること、(2)また、高頻度で（毎 30 秒ごとに）次々と観測されるために事後・予測分析を短時間で終える
必要があり、この点でもモデリングおよび計算手法に工夫を要すること、これらの統計学上の問題意識に加えて、(3)応用上でも同種のデータ分析の需要が高まっており、主にマー
ケティングの問題として産業においても関心が高まっていることを述べた。提案モデル・
手法の解説においては、実際のデータに適用した結果を主に事後・予測分布を図示することで説明した。
当日の質疑応答では、主にデータの構造について確認の質問を多く受けた。たとえば、ウェブサイトへのアクセス数をどう定義するかという点については一定の恣意性があり、ここでは既定の時間内に同じ IP アドレスからの連続したアクセスがあった場合に、そのドメ
イン間のフローにカウントすると説明した。そのため、ドメイン A からドメイン A へのフロー（同一ドメイン内のフロー）も存在することや、対象のウェブサイト外への（からの）アクセスも定義できることを指摘した。
また、技術的な点については大きく省略したものの、モデルの一部には Two-way ANOVA というよく知られた統計モデルと同様の構造があり、数理統計学を専門とする先生方から数々の指摘があった。特に、ANOVA のパラメータの解釈は識別のための仮定に大きく依存することや、個別効果や独自効果の平均に制約を入れる場合には外れ値に敏感になる問題が議論された。これらの点については共同研究者とともに認識していたが、応用の研究として普段発表する際にはあまり質問されないため、今回の議論は非常に有意義であった。
発表内容のうち半分は未発表のものであるため、今回得られたコメントをもとに論文の執筆を急ぎ、早くに公表することを目指したい。

-19-

Testing for changes in income inequality in Japan
Haruhisa Nishino
Graduate School of Social Sciences, Hiroshima University.∗

A lot of attention was paid to Japanese income inequality around 2000. Tachibanaki (2005)
and Ohtake (2008) are famous works in this field. These works were originally published inJapanese. The English translations were published later. They agree that the Japanese income
inequality has risen constantly from the 1970s to 2000s. In recent years Piketty (2014) suggested
that the income inequalities in developed countries including Japan became wider using tax
statistics. This paper investigates how much the change in Japanese income inequality after
2000 is.
Ohtake and Saito (1998) argue that inequality in the 1980s and 1990s can be explained
mainly by population ageing. The baby boomer generation in Japan, who was born in 194749, became older in the 1980s and 1990s when the proportion of older workers was larger.
Consequently, the income inequality was getting wider, because some of older workers are atsenior positions and richer; other workers are not rich.
Moriguchi and Saez (2008) and World Inequality Database indicate that top income sharesare increasing in Japan, the U.S., and other developed countries after the 1980s. Pickety (2014)
also uses top income share statistics.

0.31

two-or-more

0.30

0.29

0.28

0.27

0.26

1965

1970

1975

1980

1985

1990

1995

2000

2005

2010

2015

Figure 1: Gini coefficients: Two-or-more-persons households (1963–2014)

∗

1-2-1
Kagamiyama,
HigashiHiroshima,hnishino@hiroshima-u.ac.jp-20-

Hiroshima739-8525

Japan,

Email:

As an example, Figure 1 gives Gini coefficients based on Family Income and Expenditure
Survey (Two-or-more-persons households) in Japan: (1963–2014). It indicates that inequality
in Japan does not seem to be increasing since 2000.
This paper proposes two test statistics for testing change in income inequality by usinggrouped income data in Family Income and Expenditure Survey in Japan. One is a parametric test, which assumes lognormal distribution and follows Nishino and Kakamu (2011). The
other is a nonparametric test proposed by Beach and Davidson (1983). The nonparametric test
does not need to assume a parametric distribution, which is useful. The nonparametric test isavailable only for decile data while the parametric one for quintile data also. That is, more various datasets are available for the parametric test than for nonparametric one. In addition, thenonparametric test needs calculating local variances for each group, which is usually unavailable. In this paper, local variances are extrapolated by assuming a parametric model.
These tests lead to the following interpretation that the inequalities at 1974, 1986, 1998
become wider statistically significantly for two-or-more persons households and for workers’
households in common. It indicates that inequality is increasing largely in the 1980s and 1990s.
It cannot be judged that income inequality is increasing since the 2000s.

References[1] Beach, C. M. and Davidson, R. (1983) “Distribution-Free Statistical Inference with Lorenz
Curves and Income Shares”, Review of Economic Studies, 50, 723–735.
[2] 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.
[3] Nishino, H. and Kakamu, K. (2011) “Grouped data estimation and testing of Gini coefficients using lognormal distributions,” Sankhya B, 73, 193–210.
[4] Ohtake, F. (2008) “Inequality in Japan”, Asian Economic Policy Review, 3(1), 87–109.
[5] Ohtake, F and Saito, M. (1998) “Population aging and consumption inequality in Japan,”
The Review of Income and Wealth, 44, 361–381.
[6] Piketty, T. (2014) Capital in the Twenty-First Century, translated by A. Goldhammer, Cambridge: Harvard University Press.
[7] Tachibanaki, T. (2005) Confronting Income Inequality in Japan: A Comparative Analysis
of Causes, Consequences, and Reform, Cambridge: The MIT Press.
[8] World Inequality Database: http://wid.world/ (accessed 25 September 2018).

-21-

Particle Filtering for Non-linear State-Space Models
for Wind Speeds and Directions
Naoya Hiedaa , Takayuki Shiohamab,∗aGraduate School of Engineering, Tokyo University of Science
6-3-1 Niijuku, Katsushika, Tokyo, 125-8585 JAPANb
Department of Information and Computer Technology
Tokyo University of Science
6-3-1 Niijuku, Katsushika, Tokyo, 125-8585 JAPAN

Abstract
The state space form is a useful framework for estimating unobserved state variables fromsome given observations. The applications can be found in diverse areas of natural scienceand engineering such as ecology, epidemiology, meteorology and economics and finance. The
wind speeds and directions have complex time series probability structures involving highly
non-Gaussian and nonlinear transition. In this study, we consider a simulation-based inference
using the sequential Monte Carlo methods for computing the posterior distributions for thestate variables given all available observations. We propose an alternative approach thatallows us to extend the methods of importance sampling distributions incorporating withthe class of circular Markov transition densities. The resulting methods are compared withvarious resampling schemes with real data applications.
Keywords: circular data, EM-algorithm, state-space model, particle filter

1. Introduction
Circular (or directional) data refer to data recorded as points for which directions are measured, typically in the fields of biology, geography, medicine, and astronomy. For such data,
which are usually expressed in terms of compass angles or pairs of sine and cosine variables,
the beginning and end of the scale in the domain coincide. Owing to this periodicity, analyzing
circular data is challenging because traditional statistics are not meaningful, and may even be
misleading when the particular definition of the domain is ignored. Recent developments incircular data analysis using the statistical computing software, R, are summarized in Pewsey
et al. (2013). Although most circular data are in the form of time series, little research has
been carried out in the field of circular time series analysis compared with the number ofcircular time series modeling approaches.
In general, three main approaches are used to model circular time series. The first method
is used to obtain circular-valued random variables by wrapping; one example is the wrapped
autoregressive process of Breckling (2012). The second approach is based on a link function∗

Corresponding author
Email addresses: 4417621@ed.tus.ac.jp (Naoya Hieda), shiohama@rs.tus.ac.jp (Takayuki Shiohama)

Preprint submitted to Journal of LATEX Templates

December 3, 2018-22-

that maps a line onto a circular domain, called a linked autoregressive moving average process.
This model was proposed by Fisher & Lee (1994). The last approach specifies the densityof the conditional distribution, including the Markov process of Wehrly & Johnson (1980),
Möbius transformation of Kato (2010), and hidden Markov models of Holzmann et al. (2006).
Abe et al. (2017) studied thecircular Markov process of Wehrly & Johnson (1980) and obtained
theoretical circular autocorrelation structures under simple model assumptions. According totheir results, circular autocorrelations are determined by the mean resultant length of theunderlying circular density of the process. Abe et al. (2018) considered the circular Markov
processes whose concentration parameter could be time-varying.
Many data in directional time series applications display nonlinear features such as heteroskedasticity and a nonlinear relationship between wind direction and speed. These features
become more and more relevant as the length of the observed time series increases and as
the series itself is subject to changes in the dynamic structure. In this paper, we address thecircular process of Wehrly & Johnson (1980), which allows time-varying concentration parameters. The proposed model can incorporate the time-varying autocorrelations of the observed
circular time series. For this purpose, we introduce a simple nonparametric regression modelto the model parameter with time-varying observed exogenous variables, which cause a reasonable fit of the observed time series. The proposed models are then used to illustrate howwind direction and speed are related to the time-varying parameters. For further detail on
the time series analysis of wind direction, see, for example, Breckling (2012), Ailliot et al.
(2006), and Fuentes et al. (2005).
In an applications in meteorology, bivariate data consists of wind speeds and directions
are often modeled by using projected normal distributions (Mardia & Jupp (2009)). However,time series modeling for such a bivariate dataset is not suﬃciently studied in the literature.
The dataset of wind speeds and directions are called cylindrical data because circular winddirection data are observed along with linear wind speeds ones. Lagona et al. (2015) proposeda hidden Markov model for analyzing cylindrical time series, and the proposed model canadequately explain circular–linear correlation, and temporal autocorrelation of the observed
data. The state space modeling using circular random variable is considered in Mazumder& Bhattacharya (2017) and Kurz et al. (2016). In this study, we extend existing circularstate-space models to cope with cylindrical time series.
Sequential Monte Carlo (SMC) methods are the set of simulation-based methods which
provide a convenient and attractive approach to computing posterior distributions. Overthe last few years, there has been a proliferation of scientific papers on SMC methods andtheir applications. Several closely related algorithms, under the names of bootstrap filters,condensation, particle filters, Monte Carlo filters, interacting particle approximations andsurvival of the fittest, have appeared in several research fields. In general, the parameter
estimation in Gaussian and linear state-space model can be done by usual maximum likelihood
estimation (MLE). However, the EM algorithm is used for estimating model parameters inthe state space model, and it turns out be more robust than the direct solution of the MLE.
The recent development of EM algorithms in SMC methods are explained in Kantas et al.
(2015).

-23-

The Dantzig selector for a linear model ofdiﬀusion processes
Kou Fujimori
Waseda University
Let us consider the following model given by the linear stochastic diﬀerential
equation:∫t

Xt = X0 +Θ⊤ ϕ(Xs )ds + σWt ,(1)0where {Wt }t≥0 := {(Wt1 , . . . , Wtp )}t≥0 is a p-dimensional standard Brownian
motion, Θ is a p × p sparse deterministic matrix, σ = diag(σ1 , . . . , σp ) is a p × p
diagonal matrix and ϕ(x) = (ϕ1 (x1 ), . . . , ϕp (xp ))⊤ for x = (x1 , . . . , xp )⊤ ∈ Rp
is a smooth Rp -valued function. We will propose some estimators for the true
values (Θ0 , σ 0 ) of (Θ, σ) based on the observation of {Xt }t≥0 at n+1 equidistanttime points 0 =: tn0 < tn1 < . . . < tnn , under the high-dimensional and sparsesetting, i.e., p ≫ n and the number of nonzero components of the true value Θ0
is relatively small.
To deal with high-dimensional and sparse parameters, various kinds of estimators for regression models have been discussed. One of the most famous
estimation methods is the l1 -penalized method called Lasso proposed originally by Tibshirani (1996), which has been studied for regression models withhigh-dimensional and sparse parameters in various models including the ones ofstochastic processes.
On the other hand, a relatively new estimation procedure called the Dantzig
selector was proposed for linear regression models by Candés and Tao (2007) as
follows.
{}pj⊤β̂D := arg min ∥β∥1 , C := β ∈ R : sup |Z (Y − Zβ)| ≤ λ ,
1≤j≤p

β∈C

where λ ≥ 0 is a tuning parameter. When λ = 0, the Dantzig selector coincides
with the classical estimators such as the LSE in general cases and the MLE in
Gaussian noise cases. For λ > 0, the Dantzig selector searches for the sparsest
β within the given distance of the classical estimators. The Dantzig selector has
been studied well especially for i.i.d. models. For example, Bickel et al. (2009)
showed that the Dantzig selector has some properties similar to Lasso estimator for linear regression models in the sense of the consistency. In addition, as
well as Lasso, the Dantzig selector has variable selection consistency for some
regression models. Fan et al. (2016) showed the varable selection consistency of-24-the Dantzig selector for general single index models by using the irrepresentableconditions which are obtained from the KKT condition of the optimization
problem. The Dantzig selector also has a good potential to be applied for othermodels including the models of stochastic processes. For instance, Antoniadis
et al. (2010) applied this method to estimate regression parameter for Cox’s
proportional hazards model and proved the obtained estimator has the consistency. Fujimori (2017) studied the variable selection consistency of the Dantzig
selector for the proportional hazards model and construct asymptotically normal estimators for the regression parameter and the cumulative baseline hazard
function. Moreover, it is well-known that the Dantzig selector for linear models
has computational advantages since it can be solved by a linear programming,
while Lasso demands a convex program.
This talk dealt with the Dantzig selector for the linear models of stochastic processes (1) to estimate the drift matrix Θ0 and prove the consistency inthe sense of lq norm for every q ∈ [1, ∞] and the variable selection consistency
under some appropriate conditions. Moreover, using the variable selection consistency, construction of a new estimator which has an asymptotic normality
was discussed. We can prove the consistency of the Dantzig selector by thestandard way which is similar to Bickel et al. (2009). However, since dealing
with the KKT conditions of the Dantzig selector for our model is more diﬃcult
due to the complicated structure of the model than those of i.i.d. models, it
may be hard to obtain the same results as Fan et al. (2016) concerning thevariable selection consistency. Therefore, the another type of variable selection
consistency by using a thresholding method was proposed in this talk. This talkis based on the paper Fujimori (2018)

References
Antoniadis, A., Fryzlewicz, P. and Letué, F. The Dantzig selector in Cox’s
proportional hazards model. Scand. J. Stat. 37, no.4, 531-552. (2010).
Bickel, P. J., Ritov, Y. and Tsybakov, A. B. Simultaneous analysis of lasso and
Dantzig selector. Ann. Statist. 37, no. 4, 1705-1732. (2009).
Candés, E. and Tao, T. The Dantzig selector: statistical estimation when p ismuch larger than n. Ann. Statist. 35, no.6, 2313-2351. (2007).
Fan, Y., Gai, Y. and Zhu, L. Asymtotics of Dantzig selector for a general singleindex model. J. Syst. Sci. Complex. 29, no.4, 1123-1144. (2016).
Fujimori, K. Cox’s proportional hazards model with a high-dimensional andsparse regression parameter. arXiv:1710.10416[math.ST]. (2017).
Fujimori, K. The Dantzig selector for a linear model of diﬀusion processes. To
appear in Stat. Inference Stoch. Process. (2018).
Tibshirani, R. Regression shrinkage and selection via the Lasso. J. Roy. Statist.
Soc. Ser. B 58, no.1, 267-288. (1996).

-25-多様な関数を用いた経時変動の罰則付推定法
中京大学 国際教養学部 永井 勇
各個体に対して経時的に測定することで得られるデータは経時測定データと呼ばれ, 多くの分野で収集され分析がされている. この経時測定データの分析では, 経時変動と呼ばれるデータに潜む変動を上手く捉えることが一つの目的である. 本講演では, i 番目の個体における第 j 時点 tj での測定値を yi (tj ) と表し, t1 , . . . , tp は全ての個体で共通とした.
このとき, 未知の q 個の係数 m1 , . . . , mq を用いて m1 + m2 tj + · · · + mq tjq−1 とすると, 経時変動を測定時点の (q − 1) 次の多項式で推定できる. ここで m = (m1 , . . . , mq )′ とすると,この項は (1, tj , . . . , tq−1
)m となる. この項を平均的な経時変動を表す項とした. 一方で, 各j
個体の性別などからなる k 個の説明変数の各変数に対しても同様に(q − 1) 次の多項式で経()
時変動を推定することを考えると, 未知の係数 ξij を用いて

q
∑
i=1

ξ1i ti−1
j ,...,

q∑i=1

ξki tji−1 の形でモデル化できた. ここで, (i, j) 成分が ξij からなる k × q 未知行列を Ξ とすると, この項は
(1, tj , . . . , tq−1
)Ξ′ となり, これが i 番目の個体の説明変数 ai (k 次元ベクトル) に関する経時変
j
動を表現する項とした. この二つの項の表現より, X の j 行目を (1, tj , . . . , tjq−1 ) とし, i 番目の個体の経時測定データ yi = (yi (t1 ), . . . , yi (tp ))′ の経時変動全体を測定時点の (q −1) 次多項式で捉えられるとすると, 誤差ベクトル εi を用いて, 各個体に対して yi′ = Xm+XΞ′ ai +εiとモデル化できた. ここで, εi は E[εi ] = 0p (p 次元ゼロベクトル), Cov[εi ] = Σ (p × p 未知
正定値行列), εi ⊥⊥ εj (j ̸= i) とした. さらにここで, n を個体数, 1n を全成分が 1 の n 次元ベクトル, i 行目が i 番目の経時測定データ yi′ からなる n × p 既知行列 Y , i 行目が i 番目の個体の説明変数 a′i からなる n × k 既知行列 A (ただし A′ 1n = 0k ), i 行目が ε′i からなる n × p
誤差行列 E を用いると, 前述の各個体に対するモデルは次の形でまとめることができた;
Y = 1n m′ X ′ + AΞX ′ + E.
このモデルで q 次元未知ベクトル m と k × q 未知行列 Ξ を推定することで, Y の経時変動全体を測定時点 t1 , . . . , tp の (q − 1) 次多項式で捉えることができる. また, このモデルはPothoﬀ and Roy (1964) により提案された GMANOVA (一般化多変量分散分析) モデルと同じ形である. さらに, X の (i, j) 成分を tj に関する柔軟な関数に置き換えることで, それらの関数の重み付和により経時変動全体を捉えることも考えられる. このとき, 従来の推定法では Y へ過剰に適合してしまうが, これを回避する推定法を Nagai (2011) で提案した.
GMANOVA モデルを用いた場合, Y の経時変動全体を推定する際に用いる関数は, 平均的な経時変動の項 (= 1n m′ X ′ ) も 各説明変数に関する経時変動の項 (= AΞX ′ ) も同じ関数系 (例えば, 同じ次数の多項式) を用いた形に限られる. そこで本講演では, 平均や説明変数により異なる関数系を用いた推定を考えた. つまり例えば, 平均的な経時変動は (q0 − 1)次多項式で, いくつかの説明変数 A1 に関する経時変動は (q1 − 1) 次多項式で, 残りの説明
変数 A2 に関する経時変動は (q2 − 1) 次多項式で推定することを考えた. このとき, 上記の
GMANOVA モデルと同様に考えると, 次のような Extended GMANOVA モデル (Kollo &von Rosen (2005) Deﬁnition 4.1.3 を少し変形したモデル) による推定を考えることになった;
Y = 1nµ′X0′ + A1Ξ1X1′ + A2Ξ2X2′ + E,
ここで µ は q0 次元未知ベクトル, Xi は j 行目が (1, tj , . . . , tqj i −1 ) などからなる p × qi 既知行
列 (i = 0, 1, 2; rank(Xi ) = qi ≤ p), Ai は n × ki 既知行列 (i = 1, 2; rank(Ai ) = ki ≤ n, た
だし A′i 1n = 0ki ), Ξi は ki × qi 未知行列 (i = 1, 2), E = (ε1 , . . . , εn )′ は n × p 誤差行列であ-26-

る. また, 元の Kollo and von Rosen (2005) での定義や推定などでは 1n , A1 , A2 や X0 , X1 ,X2 に対応する部分などに階層構造を仮定しているが, 本講演ではこの構造を仮定しなかった. この µ, Ξ1 , Ξ2 を推定することで, 各項に対して (qi − 1) 次多項式などを用いて経時変動全体が推定できる. 例えば, 平均的な経時変動は 3 次の, 性別に関する経時変動は 2 次の,他の変数に関する経時変動は 4 次の多項式を用いた形で経時変動全体の推定ができる.
しかしながら, 経時変動全体が複雑な場合, 測定時点の多項式を複数組み合わせて推定しても経時変動全体を捉えることができない場合がある. そこで, 多項式だけを用いるのではなく, Nagai (2011) と同様に柔軟な関数を用いて経時変動全体を推定することが考えられる. そこで本講演では, 上記のモデルにおいて, X0 や X1 を測定時点の (q0 − 1) 次および
(q1 − 1) 次の多項式を用いる項, X2 を柔軟な関数を用いる項として考えた. これにより, 平均的な経時変動やいくつかの説明変数 A1 に対しては測定時点の多項式を用い, 別の説明変
数 A2 に対しては柔軟な関数を用いて経時変動が推定できる. しかし, 多項式の部分を単に柔軟な関数に置き換えて従来の推定法で未知の µ や Ξ1 , Ξ2 を推定すると, 目的である経時変動全体ではなく経時測定データ Y へ過剰に適合してしまうという問題が起きる. 本講演
では, Nagai (2011) と同様に多変量一般化リッジ回帰による推定法 (Yanagihara, Nagai and
Satoh, 2009) を用いて, この問題を回避することを考えた.
しかし, 従来の推定法では Ξi の推定量を得るために Ξj (i ̸= j) が必要となるため, 柔軟な関数を用いて推定する項 (= A2 Ξ2 X2′ ) だけに着目し, Ξ2 の推定量の部分にのみ罰則を付けることが難しい. そこでまず本講演では, vec 作用素などの性質 (例えば Lütkepohl (1996) 参
照) を用いて, Ξ1 の推定量を用いずに Ξ2 の推定量を得る手法を提案した. その結果, Nagai(2011) の元となった Hoerl and Kennard (1970) により提案された手法を用いることで, 柔軟な関数を用いる部分の推定量に罰則を付けることが可能であることを示した.
次に, 過剰適合を回避するために導入した罰則パラメータの最適化について考えた. そのために, Mallows (1973) や Yanagihara and Satoh (2010) などと同様の手法で, Cp 型情報量
i.i.d.
規準を構築し, さらに n − k1 − k2 − p − 2 > 0 と εi ∼ Np (0p , Σ) を仮定した下での Cp 型情報量規準のバイアスを補正した規準も構築した. そして最後に, これらの情報量規準を工夫して展開することで情報量規準を最小にする罰則パラメータが陽に得られることを示した.
数値実験による比較などについては当日報告した.
引用文献:[1] Hoerl, A. E. & Kennard R. W. (1970). Ridge regression: biased estimation for
nonorthogonal problems. Technometrics, 12, 55–67.
[2] Kollo, T. & von Rosen, D. (2005). Advanced Multivariate Statistics with Matrices,
Springer.
[3] Lütkepohl, H. (1996). Handbook of Matrices, John Wiley & Sons.
[4] Mallows, C. L. (1973). Some comments on Cp . Technometrics, 15, 661–675.
[5] Nagai, I. (2011). Modiﬁed Cp criterion for optimizing ridge and smooth parameters in
the MGR estimator for the nonparametric GMANOVA model. Open J. Stat., 1, 1–14.
[6] Potthoﬀ, R. F. & Roy, S. N. (1964). A generalized multivariate analysis of variancemodel useful especially for growth curve problems. Biometrika, 51, 313–326.
[7] Yanagihara, H., Nagai, I. & Satoh, K. (2009). A bias-corrected Cp criterion for optimizing ridge parameters in multivariate generalized ridge regression. Jpn. J. Appl.
Stat., 38, 151–172 (in Japanese).
[8] Yanagihara, H. & Satoh, K. (2010). An unbiased Cp criterion for multivariate ridge
regression. J. Multivariate Anal., 101, 1226–1238.

-27-

．序論

傾向性仮説と変化点モデルの様々な応用
広津千尋（明星大学連携研究センター
広津千尋（明星大学連携研究センター）（明星大学連携研究センター）

1

用量反応解析においては，通常, 厳密な反応曲線を想定することが難しいために単調性，凸性，Ｓ字性のような形状制約がよく設定される．それらのうち正規分布モデルに対する単調仮説については， isotonic regression がよく知られている . しかしながら，それはBartholomew によってやや直観的に導入され，とくにこのような制約のある母数空間に対する最適性は自明ではない. さらに制約付き最小二乗法は, 計算及び分布論が複雑で, 凸性,
Ｓ字性問題, 正規分布以外の確率モデル, あるいは 2 元配置交互作用問題への拡張には困難が伴う. 一方，著者等のアプローチは Hirotsu (1982)で導かれた一般制約仮説に対する検定の完全類を基にしており, その意味での最適性を持っている. それによると単調性, 凸性,
Ｓ字性仮説それぞれに対し, 単純, 2 重, 3 重累積和に基づく単調増大な統計量が示唆される.
本稿ではそのうち規準化最大対比，および規準化二乗和を用いる方法について論ずる.
一方, これらの形状制約は変化点モデルと密接な関係がある. すなわち, 単調性, 凸性,
Ｓ字性仮説はそれぞれ , 段差変化点 , スロープ変化点 , 変曲点モデルと対応する (Hirotsu
and Marumo, 2002). 例えば, 段差変化点モデルを表す対比は単調対比の典型であり, 逆にすべての単調対比は段差変化点対比の一意正係数線形結合で表される. 例えばＰＭＤＡでは副作用自発報告が収集され，その経年変化が解析されている.その場合, 単調増加傾向を
いち早く検出すると同時に, 増加傾向の生じた時点を推測することは応用上も大変有意義である. 同様のことが凸性, Ｓ字性問題についても示される. 最大対比統計量は, これら変化点モデルに対する efficient score 検定を与える. すなわち, 従来統計学の二つの異なる流れの中で研究されてきた制約仮説と変化点問題を統合的に扱う事が出来る．
以上のように，累積和に基づく方法は様々な確率モデル，様々な実質科学上の問題を理論的に一貫した方法で扱い，効率の良い計算アルゴリズムを与えることが出来る．正規分
布，2 項分布，Poisson 分布についてはこれまでいろいろな機会に発表しているので，本論では最近の研究である独立な2 × 2表の系列への応用について述べる．重要な例としてケー
ス・コントロール研究がある．

．問題の概要2

二つの 2 項分布母集団を 層上で比較する問題を考える．オッズ比の一様性が仮定出来る場合は例えば Gart (1970)により共通オッズ比に関する推論を行えば良い．一様性検定については最尤法や，Breslow & Day 検定が良く知られている他，Zelen (1971)の正確検定も利用出来る．
正確推論に関しては Agresti (1992)のサーベイが参考になる．本論はとくに層

然な順序がある場合に，単調性，および凸性を対立仮説とする検定について論ずる．
これらは上記の論文では論じられていないが，単調性に対する最大対比検定は Hirotsu 他(2001)の特殊ケースに当たり，その特殊ケースはまた太田他(2003)で詳しく論じられている．
に自

-28-

礎統計量の規準化に用いられている漸近分散を正確分散に置き換える改良を行い，Fleiss 他(2003)のロジット線形モデルによる交互作用解析と比較する．その結果は良く整合する．次にオッズ比に関する凸性仮説検定の定式化を行い，新たに最大対比正確検定を提案する．これら二つの問題について累積χ 乗検定も与える．

本論ではまずこの検定について，基．定式化
第層のデータを ，その確率を= 1で表す，ただし， = 1, 2 は二つの 2 項
分布母集団， = 1, 2は成否の反応，そして = 1, ⋯ , は層を表す．この時， を並べたベ
クトルをとして，確率分布は  = ∏ ∏  ! ∏  / !!$で与えられる．ここ= ) ++ +, +- を想定した上，周辺和 および を与えた
で 対数線形模型：log
相似検定を考えると，条件付き確率分布は , . = ∏ / - 1  exp -  という簡単な形になる．ただし，は改めて を並べたベクトルとし，- が推論対象である対数オッズ比である．この先，- に関する単調性，および凸性検定は 2 項分布や，Poisson
分布等の 1 母数指数分布族に対する方法(Hirotsu，2017)にならい，カーネルおよび漸化式
遂行の不等式の変更を適切に行うことにより実行出来る．
3





∙




"#

%
"#







∙


"#




%
"#

0#



∙

##

∙

##

##

参考文献
Agresti, A. (1992). A survey of exact inference for contingency tables. Statistical Science 7,131-177.
Fleiss, J. L., Levin, B. and Paik, M. C. (2003). Statistical methods for rates and proportions.
Wiley Series in Probability and Statistics, New York.
Gart, J. J. (1970). Point and interval estimation of the common odds ratio in the combinationof 2 x 2 tables with fixed marginals. Biometrika 57, 471-475.
Hirotsu C. (1982). Use of cumulative efficient scores for testing ordered alternatives indiscrete models. Biometrika 69, 567-577.
Hirotsu, C. (2017). Advanced analysis of variance. Wiley Series in Probability andStatistics, New York.
Hirotsu, C. and Marumo, K. (2002). Change point analysis as a method for isotonic
inference. Scand. J. Statist. 29, 125–138.
Hirotsu, C., Aoki, S., Inada, T. & Kitao, Y. (2001). An exact test for the association between the
disease and alleles at highly polymorphic loci with particular interest in the haplotype
analysis, Biometrics 57, 0769-778.

太田絵里，青木敏，広津千尋(2003). 2 × 2 × 分割表における単調仮説の検定．応用統計学 32，107－．

126

Zelen, M. (1971). The analysis of several 2 x 2 contingency tables. Biometrika 58, 129-137.

-29-

Characterizations of indicator functions of fractional factorial designs
神戸大学大学院理学研究科

青木 敏

1. はじめに実験計画法は計算代数統計における主要な研究対象の一つであり，中でも，Fontana, Pistone and Rogantin (2000, JSPI) により導入された指示関数 (indicator function) による一部実施計画の特徴付けは，古典的な接近法のひとつである．指示関数の構造は，2 水準の一部実施計画に関しては完全に解明されている．例えば，レギュラーな一部実施計画 D1の定義関係 x1 x2 x4 = x1 x3 x5 (= x2 x3 x4 x5 ) = 1 は，その指示関数 f1 の各項に対応する．
また，一部実施計画 D2 の指示関数 f2 の形は複雑であるが，定数項はこの計画の一部実施度 (3/8) を，1 次の項は因子 x4 以外の水準数が同数であることを，2 次の項は因子 4 が他の 3 因子と直交することを，それぞれ表している．
D1
x1
1
1
1
1
−1
−1
−1
−1

x2
1
1
−1
−1
1
1
−1
−1

x3
1
−1
1
−1
1
−1
1
−1

x4
1
1
−1
−1
−1
−1
1
1

x5
1
−1
1
−1
−1
1
−1
1

D2
x 1 x2 x 3 x4
1
1
1
1
1
1 −1 −1
1 −1
1 −1
−1
1 −1 −1−1 −1
1 −1−1 −1 −1
1D3
x1 x2 x 3
0 0 0
0 1 2
0 2 1
1 0 2
1 1 1
1 2 0
2 0 1
2 1 0
2 2 2

1 1
+ (x1 x2 x4 + x1 x3 x5 + x2 x3 x4 x5 )
4 4
3 1
1
1
3
f2 (x1 , x2 , x3 , x4 ) = − x4 + (x1 x2 + x1 x3 − x2 x3 ) + (x1 x3 x4 + x2 x3 x4 ) + x2 x3 x4
8 8
8
8
8
このように，2 水準計画については，分解能や aberration の概念は，指示関数の係数と直接的な関係がある．一方で多水準計画に関しては，水準を有理数とする標準的な設定においては，多水準，あるいは水準数が因子ごとに不均一な場合などの一般的な計画に対する指示関数の構造は解明されていない．例えば計画 D3 は，定義関係が x1 +x2 +x3 = 0 (mod 3)のレギュラーな一部実施計画であるが，その性質を指示関数 f3 から読み取るのは容易ではない．
f1 (x1 , x2 , x3 , x4 , x5 ) =311f3 (x1 , x2 , x3 ) = 1 − (x1 + x2 + x3 ) + (x21 + x22 + x23 ) + (x1 x2 + x1 x3 + x2 x3 )
2
2
4
3
3
27
+ x1 x2 x3 + (x21 x2 + x21 x3 + x1 x22 + x22 x3 + x1 x23 + x2 x23 ) − (x21 x22 + x21 x23 + x22 x23 )
2
4
4
9
9
33 2
− (x1 x2 x3 + x1 x22 x3 + x1 x2 x23 ) + (x21 x22 x3 + x21 x2 x23 + x1 x22 x23 ) − x21 x22 x23
4
2
4

-30-

2. 指示関数の性質
水準を有理数とする一般の一部実施計画について，その性質を指示関数の係数に結びつける方法を説明する．x1 , . . . , xn を n 個の因子とする．j = 1, . . . , n について，因子 xj の水準の集合を Aj ⊂ Q とおく．Q は有理数の集合である．Aj の要素数を rj = #Aj とし，
rj ≥ 2 と仮定する（j = 1, . . . , n）．D = A1 × · · · × An を，因子 x1 , . . . , xn の完全実施計
画とよび，その部分集合 F ⊂ D を一部実施計画という．計画 F の指示関数は，x1 , . . . , xnの多項式として一意的な表現f (x1 , . . . , xn ) =θa xa ,
(1)a∈Lをもつ．ただし，xa =naj
j=1 xj ,L = {a = (a1 , . . . , an ) ∈ Zn≥0 : 0 ≤ aj ≤ rj − 1, j = 1, . . . , n}であり，Z≥0 は非負整数の集合である．F の指示関数は，F 上で 1，D \ F で 0 をとる
観測値ベクトル y に対する D 上の補間多項式に他ならない．
多項式 f ∈ Q[x1 , . . . , xn ] が，ある一部実施計画 F ⊂ D の指示関数であることは，f と2
f が D 上で交絡することと同値である．従って，式 (1) で表される f が，ある計画の指示関数であるとき，2θa xa =
θ a xamod I(D)a∈La∈Lθa1 θa2 xa1 +a2 mod I(D)=a1 ∈L a2 ∈Lが成り立つ．ただし，I(D) は，完全実施計画 D 上で 0 となる Q[x1 , . . . , xn ] の多項式の集合で，
D の計画イデアルとよばれる．
 
 I(D) の被約グレブナー基底 G に関する
a1 +a2θa1 θa2 xの標準表示を r =μa xa と書く．多項式 f が，ある計画の指示
a1 ∈L a2 ∈L

a∈L関数であることは，代数方程式系 θa = μa , a ∈ L を満足することと同値である．この代数方程式系に，サイズ，直交性などの条件に対応する制約を付け加える．このようにして構築された代数方程式系の解集合は，与えられた性質をもつすべての一部実施計画に対応する．それを同値類に分類することで，与えられた性質をもつ一部実施計画の分類を得る
ことができる．

3. 23 × 3, 24 × 3 計画の一部実施計画の分類
計算代数ソフトウェア Macaulay2 により，23 × 3 計画の一部実施計画のすべての直交計
画，および，24 × 3 計画の一部実施計画のすべての強度 3 の直交計画を列挙し，同値類に分類した．いずれも，3 つの同値類に分類されることが分かった．詳細は，Aoki (2018),
arXiv:1810.08417 を参照されたい．

-31-

Pooling incomplete samples による統計解析
関東学院大経済 布能 英一郎 (Eiichiro Funo)

1. Introduction
Asano(1965) は、On estimating multinomial probabilities by pooling incomplete samplesという paper を発表した (Ann. Inst. Stat. Math., 17, 1-17)。Pooling incomplete samplesという名前は、これに由来する。後に、Johnson, Kotz, Balakrishnan(1997) は、分布論を網羅した書籍 Discrete Multivariate Distributions (Wiley) にて、第 35 章 Multinomial
disributions の中で、Incomplete and modified multinomial disributions という節を設け(第 10 節, 71-72 ページ)、ここで Asano の成果を解説している。それによると
X = (X1 , X2 , · · · , Xk ) は、multinomial(NY = (Y1 , · · · , Ym (m < k)x ; p1 , · · · , pk ) に従い、

p1
pm 
は、X とは独立で、multinomial Ny ; m
, · · · , mに従うものとする。このとき、
i=1 pii=1 pipi の最尤推定量 p̂i は
xi + yixip̂i =(i ≤ m, )p̂i =(i > m)NyNx
Nx 1 + m
j=1 xjである。
更に、次のことも知られている。

(i) p̂i は pi の minimum variance unbiased estimator である (Asano, 1965)(ii) p̂i は、自乗損失下で許容的 (布能, 2000)(iii) Pooling incomplete samples による多項分布からの 2 標本問題において、Total information = Within information + Between information (以後、これを Kullback 情報
量の直和分解と略記) が成り立つ。(布能, 2011)
注 Kullback(1959, 1968) は、2 標本問題において、仮説 H1 を「2 つの母集団は同じ」、
仮説 H2 を「2 つの母集団は異なる」に設定したとき、Between Information を H2 の下における「(最良) 推定量とパラメータの間の Kullback 距離」、Within Information
を「H1 の下での推定量と、仮説 H2 の下での推定量の間の Kullback 距離」、Total
Information を「H2 の下での推定量と、H1 の下でのパラメータ間の Kullback 距離」と定めた。そして、(通常の) 多項分布において、Kullback 情報量の直和分解 が成り立つことを示した。

-32-

2. その後の発展と課題
Pooling incomplete samples によるサンプリングは、多項分布の場合に限らず、一般の分布の下で考えることが出来る。
「カテゴリーが減少する場合、減少後のカテゴリーには、元の生起確率を比例配分したものを割り当てる」ものであるから、離散分布にはこの考え方をそのまま当てはめられる。ところが、多項分布の場合に成り立っていることが、他の分布のとき、あるいは、多項分布の場合でも、Asano が設定した Pooling incomplete samplesの枠組から外れると、成り立たないことが生じる。
2.1 unbiasedness Asano が示した方法は、Pooling incomplete samples を伴う他の分布に適用できる場合もあれば、うまくゆかない場合もある。
2.2 許容性の証明 Pooling incomplete samples を伴う多項分布の場合、Stepwise Bayes
法は、の許容性の証明法として「相性が良い」と言える。それは、Brown(1981) が「sample
space が有限、母数空間がコンパクトのとき、任意の許容的推定量に対して、sequence ofpriors(あるいは、1 つの prior) が存在し、これに対する Stepwise Bayes 解 (あるいは、1つの prior に対する Bayes 解) となっている」という定理を示しているからである。他方、
Stepwise Bayes 法は、母数空間がコンパクトでないときには、必ずしもうまくゆくとは限らない。
2.3 2 標本問題における Kullback 情報量の直和分解 多項分布以外の Pooling incomplete
samples では、成り立たない場合の方が多いと言わざるを得ない。
たとえば、Pooling incomplete samples を伴う負の多項分布の場合、「負の多項分布で自然に考えられる Pooling incomplete samples 」と思える確率モデルで計算したところ、
Kullback 情報量の直和分解が成り立たかった。他方、負の多項分布 = 負の二項分布 × 多
項分布 なので、この「多項分布」の部分で Pooling incomplete samples によるモデルを設定したところ、Kullback 情報量の直和分解が成立した。
更に、パラメータを分離して Kullback 情報量の直和分解が成立するような確率モデルの構築も行った。
今後の課題 カテゴリーが減少した際、「Kullback 情報量の直和分解が成立」するような
Pooling incomplete samples 風の確率モデルを探すことも必要だが、Kullback 情報量の直和分解が成立しない原因を突き止めることも重要なことである。願わくば、Total informationと Within information + Between information の差が、何らかの意味を持つ量であることを突き止められたらと思う。
文献 文中で十分な検索情報を与えたもの以外に、[1] 布能英一郎 (2000). Missing, Poolingを伴う離散サンプリングの統計的推測. 京都大学数理解析研究所講究録 No 1161, 141158. [2] 布能英一郎 (2011). On estimating multivariate discrete probabilities by pooledincomplete samples and related topics. 京都大学数理解析研究所講究録 No 1758.
[3]

Kullback, S.(1968). Information Theory and Statistics, Revised edition. Dover.

-33-

Pitman’s 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)

1 Introduction
We consider Pitman closeness domination in predictive density estimation problems when theunderlying loss metric is α-divergence {D(α)}, a loss introduced by Csiszàr (1967).The underlying distributions considered are normal, including the distribution of the observables, thedistribution 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. 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 topredict a 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 the σi2 , i = 1, 2 are
unknown and/or satisfy the additional order restriction, σ12 ≤ σ22 .
The class of α−divergence losses is given byZ

Dα {p̂(ỹ|y), p(ỹ|ψ)} =fαp̂(ỹ|y)p(ỹ|ψ)dỹ,p(ỹ|ψ)

(1)

where, for −1 ≤ α ≤ 1
fα (z) =
4(1+α)/2 ), 1−α2 (1 − zz log z, − log z,|α| < 1α=1
α = −1.

(2)

Here KL loss corresponds to α = −1, and α = 1 is sometimes referred to as reverse KL loss.
Chang and Strawderman (2014, JMVA, Theorem 2.1) have derived the general form of Dαloss for the case of normal models and shown that it is a concave monotone function of quadratic
loss and is also a function of the variances (observed, predicand, and plug-in).
An alternative criterion to evaluate the goodness of estimators was introduced by Pitman(1937) as follows:
Let T1 and T2 be two estimators of θ. Then T1 is closer to θ than T2 if Pitman nearness(P N ) of T1 compared to T2
P Nθ (T1 , T2 ) = P {|T1 − θ| < |T2 − θ|} > 1/2.
For the case, when the estimators are equal with positive probability, Nayak (1990),Gupta andSingh (1992) defined the modified Pitman nearness (M P N ) of T1 compared to T2 . Setting
M P Nθ (T1 , T2 ) = P {|T1 − θ| < |T2 − θ||T1 6= T2 } = P {|T1 − θ| < |T2 − θ|, T1 6= T2 }/P {T1 6= T2 }T1 is closer to θ than T2 if M P Nθ (T1 , T2 ) > 1/2.
Here is a brief review of some of the relevant literature for the problem of estimating the
P i
P imean. Let X̄i = nj=1
Xij /ni , s2i = nj=1
(Xij − X̄i )2 /(ni − 1) be the unbiased estimators of µiand σi2 , respectively, based on samples of size ni from two normal populations, X1 ∼ N (µ1 , σ12 )and X2 ∼ N (µ2 , σ22 ) respectively. When the variances are unknown, the unbiased estimatorµ̂GD =

n1 s22n2 s21X̄
+X̄2
1n1 s22 + n2 s21
n1 s22 + n2 s21

was proposed by Graybill and Deal (1959)and they gave a necessary and sufficient condition on
n1 and n2 for µ̂GD to have a smaller variance than both X̄1 and X̄2 .
When estimating the ordered means, Oono and Shinozaki (2005, JSPI) proposed truncatedestimators of µi , i = 1, 2,
GD
µ̂OS
},
1 = min{X̄1 , µ̂

GD
µ̂OS
},
2 = max{X̄2 , µ̂

(3)and showed that µ̂OSdominates the X̄i in terms of MSE if and only if MSE of µ̂GD is not largerithan that of X̄i to estimate µi when µ1 = µ2 .
-34-

When there are order restrictions given on both means and variances, µ1 ≤ µ2 , σ12 ≤ σ22 ,Chang, Oono and Shinozaki (2012, JSPI)have proposed
 OS, µ̂i 

=
µ̂CS
i




if s21 ≤ s22

2
1
X̄1 + n1n+n
X̄2 , if s21 > s22
min X̄1 , n1n+n
2
2

(4)

OS
CS
OS
They show that µ̂CS
2 stochastically dominates µ̂2 , but µ̂1 cannot dominate µ̂1 even in termof MSE when µ2 − µ1 is sufficient large. µ̂CS
is also Pitman closer to µ2 than µ̂OS22 .

2 Pitman closeness in predicting density function under the D(α) Loss MetricIn this section we will establish Pitman closeness results under the {D(α)} loss metric for certain
predictive density estimation problems involving two normal populations when the means areordered. Here we only discuss the unknown variance cases.
Definition 1. Given two predictive density estimates fˆ1 (ỹ|x) and fˆ2 (ỹ|x) of a density f (ỹ|ψ)based on data x from a distributions X ∼ g(X|ψ), ψ ∈ Ω, fˆ2 (ỹ|x) is closer to f (ỹ|ψ) than fˆ1 (ỹ|x)with respect to the D(α) metric under the modified Pitman closeness criterion, if ∀ψ ∈ Ω,Pψ {Dα (fˆ2 (ỹ|x), f (ỹ|ψ)) < Dα (fˆ1 (ỹ|x), f (ỹ|ψ))|fˆ2 (ỹ|x) =6 fˆ1 (ỹ|x)} ≥ 1/2,
with strict inequality for some ψ ∈ Ω.

2.1 Case when variances are unknown and unrestricted
It is assumed that µ1 ≤ µ2 and that no restriction is given on unknown σi2 . It is shown thatplug-in predictive density with µ̂OSi , (3), is Pitman closer to the true predictive density thanplug-in predictive density with X̄i under Dα loss. We wish to predict the density of a future
observation from a normal population with mean µi and unknown variance σ 2 , i.e. we wish topredict the densityf (ỹ) ∼ N (µi , σ 2 ).
Let
X̄i ∼ N (µi , σi2 /ni ),Si2 ∼ σi2 χ2ni −1 , i = 1, 2

are independent, where µ1 ≤ µ2 and it is desired to predict the density of a future independent
variable Yi ∼ N (µi , aσ 2 ), i = 1, 2, where a > 0 is known. We have the following main result.
c2
Theorem 1. The plug-in predictive density estimate fˆiOS (ỹ) ∼ N (µ̂OSi , aσ ) is Pitman closerX̄c2 ) for all µ ≤ µ and σ 2 , i = 1, 2 under the D(α) metricto f (ỹ|µi , aσ 2 ) than fˆi i (ỹ) ∼ N (X̄i , aσ12ic2 if and only if µ̂GD is Pitman closer to µ than X̄ forfor all −1 ≤ α ≤ 1 and every estimator σ
iall σ12 and σ22 when µ1 = µ2 = µ.

2.2 Case when variances are ordered
We give Pitman closeness domination results in predictive density estimation when unknown
variances σ12 and σ22 satisfy the order restriction σ12 ≤ σ22 .
First we estimate the density of a future observation from a normal population with mean
µ2 and variance σ 2 = aσ22 , where a is known, i.e. we estimate the densityf (ỹ) ∼ N (µ2 , aσ22 ).
c2
Theorem 2. The plug-in predictive density estimate fˆCS (ỹ) ∼ N (µ̂CS
2 , aσ2 ) is Pitman closer
c2
to f (ỹ|µ2 , aσ22 ) than fˆOS (ỹ) ∼ N (µ̂OS
2 , aσ2 ) under the D(α) metric for all −1 ≤ α ≤ 1 and forc2 .
any estimator σ
2
Next we consider estimating the predictive density with smaller variance σ12 , N (µ1 , aσ12 ).
b12 ) can not be Pitman
Theorem 3. The plug-in predictive density estimate fˆCS (ỹ) ∼ N (µ̂CS
1 , aσc2
closer to f (ỹ|µ1 , aσ12 ) than fˆOS (ỹ) ∼ N (µ̂OS
1 , aσ1 ) when µ2 − µ1 is sufficient large, under thec2 .
{D(α)} metric for all −1 ≤ α ≤ 1 and for any estimator σ
1-35-

一般欠測データの下での 2 標本問題における多変量正規母集団の同等性検定
東京理科大学・理
東京理科大学・理
東京理科大学・理

1

野村 玲実
八木 文香
瀬尾 隆

はじめに
一般欠測データの下で, 2 つの多変量正規母集団の同等性検定問題について考える．すなわち， µ(ℓ)

と Σ(ℓ) , ℓ = 1, 2 をそれぞれ第 ℓ 母集団の平均ベクトルと分散共分散行列とするとき, この検定問題
は H0 : µ(1) = µ(2) , Σ(1) = Σ(2) である平均ベクトルと分散共分散行列の同時検定を考えることと同じであり，完全データの下では，尤度比検定統計量とその修正尤度比検定統計量が与えられている(Muirhead (1982), Srivastava (2002) など参照). また 1 変量の場合については，Zhang et al. (2012)
で議論され, 正確な尤度比検定が与えられている. 一方, 欠測データの下での同時検定については, 単調型欠測の場合の議論があり, Hao and Krishnamoorthy (2001) や Hosoya and Seo (2015, 2016) などで尤度比検定統計量の導出とその帰無分布に対する近似上側パーセント点などが与えられている.
ただし, Hao and Krishnamoorthy (2001) や Hosoya and Seo (2015) では 2-step 単調欠測データの下での 1 標本問題における同時検定を議論しており, Hosoya and Seo (2016) では 2-step 単調欠測デー
タの下での多標本問題について議論している. このように単調欠測データの下での同時検定についての議論はあるが, 単調でない一般欠測データの下での同時検定の議論はない. そこで本報告では一般
欠測データの下での 2 標本問題に対する同時検定問題, すなわち, 2 つの多変量正規母集団の同等性検定問題を考え, その尤度比検定統計量を導出し, 尤度比検定を与える.

2

最尤推定量と反復法
第 ℓ 母集団 (ℓ = 1, 2) に対して互いに独立で欠測値を含む p 次元観測ベクトルがそれぞれ n(ℓ) 個あるとする. このとき H0 : µ(1) = µ(2) , Σ(1) = Σ(2) である仮説検定問題を考える．尤度比検定統計量を与えるため, 最尤推定量の導出について考える．同じ欠測パターンをもつ観測ベクトルを１つのグループとし, そのグループの総数を K (ℓ) とする. また, 第 k 番目のグループにおける観測ベクトル∑K (ℓ) (ℓ)(ℓ)の数を nk (k = 1, 2, . . . , K (ℓ) ) とする．ここに k=1 nk = n(ℓ) , ℓ = 1, 2 である．このとき, k 番目
(ℓ)

(ℓ)

のグループに属する観測ベクトル xkj (k = 1, 2, . . . , K (ℓ) , j = 1, 2, . . . , nk , ℓ = 1, 2) の欠測部分を除
(ℓ)

(ℓ)

(ℓ)

(ℓ)

(ℓ)

くような変換行列 Bk を xkj にかけ，データの変換を行う．すなわち z kj = Bk xkj と変換する(ℓ)(Srivastava and Carter (1986) 参照). この z kj に対して, 正規性を仮定し, 対数尤度関数を微分する
ことによって, µ(ℓ) の最尤推定量は (ℓ)
−1  (ℓ)

K
K
∑
∑
′
′
(ℓ) (ℓ) b (ℓ)−1 (ℓ) (ℓ) (ℓ) b (ℓ)−1 (ℓ) b (ℓ) = µnk Bk ΛBk
nk Bk Λz̄ k , ℓ = 1, 2
k
k
k=1

k=1

b (ℓ) はとなり, Σ(ℓ) の最尤推定量 Σ
(ℓ)
K
∑

k=1

′

(ℓ) (ℓ) b (ℓ)−1 (ℓ)
nk Bk Λ
Bk =
k

(ℓ)
K
∑

′

′

(ℓ) b (ℓ)−1 b (ℓ) b (ℓ) b (ℓ)−1 (ℓ)
Bk ΛUk Uk Λk
Bk
k

k=1

-36-

を満たす解となる. ただし
(ℓ)

b (ℓ)
Λ
k

(ℓ) b (ℓ) (ℓ)′
Bk Σ
Bk ,

=

(ℓ)
Zk

(ℓ) ′
b (ℓ) = Z(ℓ) − B(ℓ) µUen(ℓ) ,
k
k
k b
k

(ℓ)
(ℓ)
(ℓ)
(z k1 , z k2 , . . . , z (ℓ) ),
kn

=

(ℓ)
z̄ k =

k

nk
1 ∑
(ℓ)

nk

(ℓ)

z kα ,

α=1′

en(ℓ) = (1, 1, . . . , 1) .
k

e も求めることがe Σ
同様に帰無仮説 (H0 : µ(1) = µ(2) , Σ(1) = Σ(2) ) の下での µ や Σ の最尤推定量 µ,b (ℓ) や Σe の正確な値は求めることが困難なため, Srivastava and Carter (1986) のアイデアできる. Σ
b (ℓ) を例として以下に手順を示す. まず, Σ , µ を第 i 近による反復法を用いて数値的に求める．Σ
i
i
(ℓ)

(ℓ)

(ℓ)

(ℓ)

似値とし, ∆i を分散共分散行列の第 i 更新量とする. 初期値 Σ0 は完全データ部分から計算した分
(ℓ)

(ℓ)

(ℓ)

散共分散行列の最尤推定量とする. このとき, 反復式は Σi+1 = Σi + ∆i+1 と書くことができる. こ
こに vec(∆i+1 ) = (Qi )−1 vec(Ei ) で
(ℓ)

(ℓ)
Qi =

(ℓ)
K
∑

(ℓ)

(ℓ)

(ℓ) (ℓ)
(ℓ)
(ℓ)
(ℓ)
(ℓ)
(ℓ)
(nk Dik ⊗ Dik − Dik ⊗ Fik − Fik ⊗ Dik ),

(ℓ)
Ei =

k=1
(ℓ)

(ℓ)
K
∑

(ℓ)

(ℓ)

(ℓ)

(nk Dik − Fik )

k=1
(ℓ)

であり, Dik , Fik は
(ℓ)′

(ℓ)

(ℓ)′

(ℓ)

(ℓ)′

(ℓ)

(ℓ)′

(ℓ)

(ℓ)

= Bk (Bk Σi Bk )−1 Bk ,(ℓ)

= Bk (Bk Σi Bk )−1 Uik Uik (Bk Σi Bk )−1 Bk

Dik

Fik

(ℓ)

(ℓ)

(ℓ)′

(ℓ)

(ℓ)

(ℓ)′

(ℓ)

である. ただし, Uik = Zk − Bk µi e′ (ℓ) である. よって尤度比は
(ℓ)

(ℓ)

(ℓ)

(ℓ)

nk

(1)
K
∏

λ

=

b | 2 nk
|Λ
k
(1) 1

(1)

(2)
K
∏

k=1

k=1

(1)
K
∏

(2)
K
∏

e | 2 nk
|Λ
k
(1) 1

(1)

k=1

b | 2 nk
|Λ
k
(2) 1

(2)

e | 2 nk
|Λ
k
(2) 1

(2)

k=1

′

e (ℓ) である. 尤度比検定統計量 −2 log λ の極限分布は自由度
e (ℓ) = B(ℓ) ΣBで与えられる. ただし, Λkkkp(p + 3)/2 のカイ二乗分布であり, その近似精度をモンテカルロ・シミュレーションにより評価した.

参考文献[1] Hao, J. and Krishnamoorthy, K. (2001). Inferences on a normal covariance matrix andgeneralized variance with monotone missing data. Journal of Multivariate Analysis, 78,
62-82.
[2] Hosoya, M. and Seo, T. (2015). Simulaneous testing of the mean vector and the covariancematrix with two-step monotone missing data. SUT Journal of Mathematics, 51, 83-98.
[3] Hosoya, M. and Seo, T. (2016). On the likelihood ratio test for the equality of multivariate
normal populations with two-step monotone missing data. Journal of Statistical Theory andPractice, 10, 673-692.
[4] Muirhead, R. J. (1982). Aspects of Multivariate Statistical Theory. Willey, New York.
[5] Srivastava, M. S. (2002). Methods of Multivariate Statistics. Willey, New York.
[6] Srivastava, M. S. and Carter, E. M. (1986). The maximum likelihood method for nonresponse in sample surveys. Survey Methodology, 12, 61-72.
[7] Zhang, L., Xu, X. and Chen, G. (2012). The exact likelihood ration test for equality of two
normal populations. The American Statistician, 66, 180-184.

-37-単調欠測データが一様共分散構造を持つ場合の平行性仮説検定と水準差の信頼区間
東京理科大学・理
東京理科大学・理
東京理科大学・理
防衛大学校・数学

佐伯 悠一郎
八木 文香
瀬尾 隆
百武 弘登

はじめに

1

本研究ではプロフィール分析における平行性仮説検定とその水準差について考える．平行性仮説検定問題は，単調欠測データの下で µ(1) , µ(2) , · · · , µ(g) を g 個の母集団の平均ベクトルとすると，次のように表される．
H : µ(ℓ) − µ(g) = γ (ℓ) 1p , ℓ = 1, 2, . . . , g − 1 vs. A : not H,
ただし，γ (ℓ) は水準差と呼ばれる未知パラメータであり，1p = (1, 1, · · · , 1)′ とする．平行性仮説検定
は Srivastava (1987) 等で完全データに対して尤度比検定統計量と漸近分布を用いた近似パーセント点が与えている．また欠測データの下での議論については Onozawa, Takahashi and Seo (2013) で 2-step
単調欠測データに対し，2 群におけるホテリングの T 2 型検定統計量を用いた検定と，多群に対する尤
度比検定統計量を与えている．どちらも検定統計量の漸近分布と近似上側パーセント点を提案している．一方，Sawatsuhashi and Hyakutake (2016) では観測ベクトルにランダム効果モデルと一様共分散構造を仮定した場合の 2 群における，一般の単調欠測データの下でのプロフィール分析が議論され，検定統計量とその正確な分布が得られている．一様共分散構造とは繰り返し測定データに対して仮定される分散共分散行列の構造の一つである．本報告では Sawatsuhashi and Hyakutake (2016) を基に，一様構造を持つ一般の単調欠測データに対して，多群の場合の平行性仮説検定と，2 群かつ 2-step 単調欠
測データに対し水準差の信頼区間について考える．

平行性仮説検定

2

(ℓ)

(ℓ)

(ℓ)

k-step 単調欠測データとして，xij ∼ Npi (µi , Σi ) (i = 1, 2, . . . , k , j = 1, 2, . . . , ni ) を第 ℓ 群の(ℓ)

第 i ステップ内 j 番目にある観測ベクトルとする．ただし，p = p1 > p2 > · · · > pk である．また µi
は µ(ℓ) の第 1 から第 pi 成分で構成される (pi × 1) ベクトルであり，Σi は Σ の pi 次主座小行列とする．
本報告では分散共分散行列の構造をとし，尤度比検定統計量の正確な分布を導出した．
{}
Σ = σ 2 (1 − ρ)I p + ρ1p 1′p , ただし，− 1/(p − 1) < ρ < 1.
観測ベクトルを C i 1pi −1 = 0 , C i C ′i = I pi を満たす (pi − 1) × pi の対比行列 C i を用いて変換する．
具体的な行列 C i については Sawatsuhashi and Hyakutake (2016) などを参照されたい．変換後の観
(ℓ)
(ℓ)
(ℓ)
(ℓ)
(ℓ)
測ベクトルは y ij = C i xij ∼ Npi −1 (η i , τ0 I pi ) となる．ここで，η i = C i µi , τ0 = σ 2 (1 − ρ)∑i(ℓ)
である．さらに pk−i+1 =h=1 qh となる qm (m = 1, 2, . . . , k) に対し，y ij を (q1 − 1) 次元，q2(ℓ)

次元 , . . . , qk 次元ごとに区切り z ms とする．m は縦に分割後のブロック番号であり，s はサンプル番∑k−m+1 (ℓ)(ℓ)(ℓ)
ni とする．µ(ℓ) についても同様に分割し ξ m とすると号を意味する．ただし，Mm =
i=1
(ℓ)

(ℓ)

(ℓ)

(ℓ)

z 1s ∼ Nq1 −1 (ξ 1 , τ0 I q1 −1 ) , z ms ∼ Nqm (ξ m , τ0 I qm ), m = 2, 3, . . . , k となる．このことから，一般のk-step の場合の尤度関数を与えることができるので H : ξ (1) = ξ (2) = · · · = ξ (g) vs. A : not H の下で
(ℓ)
(ℓ)
∑g
∑g ∑Mm(ℓ)(ℓ)(ℓ) ∑ m(ℓ)(ℓ)
の尤度比は Mm = ℓ=1 Mm , z m = (1/Mm ) M
s=1 z ms , z m = (1/Mm )s=1 z ms とするとℓ=1(
λk = 1 +gk ∑∑m=1 ℓ=1

(

)(

g M
k ∑
m
∑
) ∑
( (ℓ)
)′ ( (ℓ)))−(1/2)νk(ℓ)
z ms − z (ℓ)
z
−
z
m
ms
m
(ℓ)

′
(ℓ)
z m − z (ℓ)
z m − z (ℓ)
Mm
m
m /

m=1 ℓ=1 s=1-38-となる．また，νk = M1 (q1 − 1) + M2 q2 + · · · + Mk qk である．さらにこの尤度比を用いた検定統計量
Tk の分布について)
νk − g(p − 1) ( − ν2kTk =λk − 1 ∼ F(g−1)(p−1), νk −g(p−1)(g − 1)(p − 1)
であることが示される．補足として，いくつかのパラメータ設定でのモンテカルロ・シミュレーションを行ったところ，数値的にも F 分布に一致することを確認している．

3

水準差の信頼区間
(ℓ)

(ℓ)

(ℓ)

観測ベクトル xij (i = 1, 2 , j = 1, 2, . . . , ni , ℓ = 1, 2) は Npi (µi , Σi ) に従うとし，平行性(1)µ − µ(2) = γ1p が認められているとする．欠測データにおけるパラメータ ρ, σ 2 の MLE(最尤推定量)は容易に求められないので，これらを既知とみなし最後に推定量を代入する．n∗i = (ni ni )/(ni +
(ℓ)
(2)
(ℓ)
(ℓ) ∑ni(ℓ)(1)(2)
ni ), xi = (1/ni ) j=1
xij , ui = xi − xi とおくと ρ が既知の場合の γ の MLE は γ
b =
(ℓ) (
∑2
∑
∑
ni
(ℓ)
(ℓ)
2
∗ −1 ′∗ −1
i=1 ni τ1i 1pi ui /i=1 ni τ1i pi である．ただし，τ1i = 1+(pi −1)ρ である．ここで V i =j=1 xij −∑∑(ℓ) )( (ℓ)(ℓ) )′(ℓ)(ℓ)
xi
xij −xi, Vi = 2ℓ=1 V i , S i = (1/ni −2)V i , ni = 2ℓ=1 ni と定義し，Siotani, Hayakawa
and Fujikoshi ∑(1985) で与えられている完全データの下でのMLEを基に，パラメータ ρ と σ 2 の推定
∑
∑
∑
σ2, σ
e2 = 2i=1 tr(V i )/ 2i=1 (ni − 2)pi を提案す
τ2i − tr(S i ))/ 2i=1 {pi (pi − 1)}e
量として ρe = 2i=1 (be2 は共に MLE ではないが，どちらも 2-step 単調欠測データのる．ただし τb2i = 1′pi S i 1pi とする．ρe, σ
全ての観測値を使用した統計量である．Jensen の不等式とこれらの推定量を用いて γ の 100α% 近似信
頼区間
( 2
( 2
)− 1
)− 1
2
2
∑
∑
n∗i p2in∗i p2i
γ
e−
<
γ
<
γ
e
+
τb2i t2ni −2 ( α2 )τb2i t2ni −2 ( α2 )(1) (2)

i=1

∑
−1 ∗
−1 ∗ ′
ni p i ,e1i
ni 1pi ui / 2i=1 τe1i
i=1 τ

∑2

が得られた．ただし，γe=は自由度 ni − 2 の t 分布の上側 100(α/2)% 点である．

4

(1)

i=1

τe1i = 1 + (pi − 1)eρ であり，tni −2 (α/2)

シミュレーション結果

前節で導出した近似信頼区間を C1 とし，完全データ部分のみを使用した信頼区間を C2 とすると，
11100α% 信頼区間 C2 は γ
b1 − tn1 −2 (α/2)((n∗1 p21 )/b
τ21 )− 2 < γ < γ
b1 + tn1 −2 (α/2)(n∗1 p21 /b
τ21 )− 2 である．
有意水準 5%，繰り返し回数を 1,000,000 回とし，被覆確率を比較した．結果として C1 の方が被覆確率が高く，かつその確率は 0.95 を大きく上回らない，及び狭い区間であることが確認できた．

参考文献[1] Onozawa, M., Takahashi, S. and Seo, T. (2013). Tests for profile analysis based on two-stepmonotone missing data. Discussiones Mathematicae Probability and Statistics, 33, 171-190.
[2] Sawatsuhashi, K. and Hyakutake, H. (2016). Profile analysis for random eﬀects model in two
sample problem with monotone missing. Far East Journal of Theoretical Statistics, 52, 235-251.
[3] Siotani, M., Hayakawa, T. and Fujikoshi, Y. (1985). Modern Multivariate Statistical Analysis:
A Graduate Course and Handbook, American Sciences Press, Inc., Columbus, Ohio.
[4] Srivastava, M. S. (1987). Profile analysis of several groups. Communications in Statistics -Theory
and Methods, 16, 909-926.

-39-

2-step 単調欠測データのもとでの部分平均ベクトルの検定に対する修正尤度比検定統計量

東京理科大・理
東京理科大・理

1

川崎 玉恵
瀬尾 隆

はじめに

１標本，２標本問題における部分平均ベクトルの仮説検定問題を 2-step 単調欠測データのもとで議論した．部分平均ベクトルの仮説検定問題について，欠測値を含んでいない完全デー
タの場合は，尤度比を用いた検定統計量が１標本問題では Rao (1949) や Giri (1964) などで議論されており，これらは Rao の U 統計量と呼ばれている．欠測値を含むデータに対しては，Kawasaki and Seo (2016) が１標本問題における 2-step 単調欠測データのもとでの最尤
推定量と尤度比検定統計量を導出している．
本報告では，2-step 単調欠測データのもとでの１標本，２標本問題における部分平均ベク
トルの仮説検定問題について，１標本問題では Kawasaki and Seo (2016) で与えた尤度比検
定統計量を，２標本問題では最尤推定量を導出し，尤度比検定統計量を与えたのち，検定統計量の帰無分布を大標本漸近枠組みの下で漸近展開することで，それぞれの修正尤度比検定
統計量を与えた．さらにこれらの精度について，モンテカルロ・シミュレーションにより数値的評価を行った．

2

１標本問題における部分平均ベクトルの仮説検定問題

x1 , x2 , . . . , xN1 が Np (µ, Σ) に，xN1 +1 , xN1 +2 , . . . , xN が Np1 +p2 (µ(12) , Σ(12)(12) ) に従う互
いに独立な確率ベクトルとする．ただし，µ と Σ はそれぞれ次のよう分割される．




(
)
(
)
Σ11 Σ12 Σ13µ10
Σ(12)(12) Σ(12)3
µ(12)



 Σ
µ =  µ2  =, Σ =  21 Σ22 Σ23  =(1)
Σ3(12)
Σ33µ3µ3
Σ31 Σ32 Σ33
また，xj をそれぞれ p1 × 1, p2 × 1, p3 × 1 の確率ベクトルとし，xj = (x′1j , x′2j , x′3j )′ =(x′(12)j , x′3j )′ と表す．ただし，xij , i = 1, 2, 3, j = 1, 2, . . . , N1 は pi × 1，p = p1 + p2 + p3 であり，N2 = N − N1 である．このとき以下の仮説検定問題を考えた．

H0 : µ(23) = µ(23)0 given µ1 = µ10 vs. H1 : µ(23) ̸= µ(23)0 given µ1 = µ10
この問題における尤度比検定統計量 −2 log λ は Kawasaki and Seo (2016) で与えられている．ただし)−N1 /2()(
e (12)(12) | −N/2 |Φe 33 ||Φλ=b 11 | · |Ψb 33 |
b 22 ||Ψ|Ψであり，記号の説明については Kawasaki and Seo (2016) を参照されたい．この尤度比検定
統計量における帰無分布の漸近展開を以下の大標本漸近枠組み

δi =

Ni
→ positive constants as N1 , N2 → ∞, i = 1, 2
N-40-

の下で行い，修正尤度比検定統計量を導出した．ただし，摂動展開は以下を用いる．






V
V
V
zF 1
x
F
11
F
12
F
13
F
1


 √

 √ 
z F = z F 2  = N1 xF 2  , VF = VF 21 VF 22 VF 23  = N1 − 1(SF − Ip ),xF 3
VF 31 VF 32 VF 33
zF 3
(
(
)
)
(
)
√
√
xL1
z L1
VF 11 VF 12
zL =
= N2
, VL == N2 − 1(SL − Ip(12) )z L2
xL2
VF 21 VF 22
また，xF = (x′F 1 , x′F 2 , x′F 3 )′ と SF はそれぞれ N1 × p 部分の，xL = (x′L1 , x′L2 )′ と SL はそれぞれ N2 × (p1 + p2 ) 部分の標本平均ベクトルと不偏標本分散共分散行列とする．
以上から，尤度比検定統計量 −2 log λ の特性関数は{}(p2 + p3 )cit(−2 log λ)−(p2 +p3 )/2−1E[e] = (1 − 2it)1+{(1 − 2it) − 1} + o(N −1 )2Nと展開することができた．ただし
{
}
1
1
1
1
c=
p2 (2p1 + p2 + 2) + p3 (p + 3) − p3 (p3 + 2) .
p2 + p3 2
δ1
δ1
よって，−2 log λ の帰無分布に対する漸近展開はPr(−2 log λ ≤ x) = Gp2 +p3 (x) +(p2 + p3 )c{Gp2 +p3 +2 (x) − Gp2 +p3 (x)} + o(N −1 )2Nのように与えられた．ただし，Gf (x) は自由度 f の χ2 分布の分布関数とする．さらに，バートレット修正係数を ρ = 1 − c/N と提案した．

3

２標本問題における部分平均ベクトルの仮説検定問題

2-step 単調欠測データのもとでの２標本問題における部分平均ベクトルの仮説検定問題
(1)

(2)

(1)

(2)

(1)

(2)

(1)

(2)

H0 : µ(23) = µ(23) given µ1 = µ1 vs. H1 : µ(23) ̸= µ(23) given µ1 = µ1
(i)

(2)

(i)

を考えた．ただし，µ1 , µ(23) の分割は１標本問題のときの (1) 式と同様であり，i = 1, 2 は母集団を表す添え字とする.
仮説検定問題 (2) における尤度比検定統計量を導出し，１標本問題と同様に大標本漸近枠組みのもとで漸近展開を行い，修正尤度比検定統計量の導出をおこなった．また，100 万回のモンテカルロ・シミュレーションにより本報告で与えた修正尤度比検定統計量についての数値的評価を行った．

参考文献[1] Giri, N. C. (1964). On the likelihood ratio test of a normal multivariate testing problem. Ann. Math. Stat., 35, 181–189.
[2] Kawasaki, T. and Seo, T. (2016). A test for subvector of mean vector with two-stepmonotone missing data. SUT J. Math., 52, 21–39.
[3] Rao, C. R. (1949). On some problems arising out of discrimination with multiple
characters. Sankhyā, 9, 343–364.

-41-

Recent cylindrical models and their application to tree data set
Toshihiro Abe1 & Ichiro Ken Shimatani2
1
Nanzan University, Nagoya, Japan
2
The Institute of Statistical Mathematics, Tokyo, Japan

1

Introduction

If the objective variable is circular, we need multivariate versions of circular probability distributions, and in the two-dimensional case, there are two types. The first is when two objective
variables are circular; then, the probability distribution is defined on the torus. The second
type is when one is circular and the other is linear; then, the probability distribution should be
defined on S 1 × R, or if linear variables take only non-negative values, the domain is S 1 × [0, ∞).
In both cases, the probability distribution is defined on a “cylinder” and such distributions arecalled “cylindrical distributions.” The latter type of data are commonly seen, for example, invelocity data, which consist of a speed and a direction. Often, these show specific correlations.
In such cases, correlations are present between the linear and circular variables in the form thatthe variances over the circular variables become high when speeds are slow. In this talk, by
focusing on cylindrical distributions on S 1 × [0, ∞), we introduce those distributions that aremathematically tractable, namely the probability density function (pdf) can be expressed byelementary functions or, at most, well-known special functions, and provide examples of theirapplications.

2

Cylindrical distributions

Let (Θ, L) be a pair of circular and linear random variables, respectively. When we observecylindrical data, if the linear components increase, the degree of concentration around a certain
direction also tends to increase. This section introduces models that satisfy these conditions.
Johnson and Wehrly (1978) proposed a cylindrical distribution whose joint pdf is given by√ν 2 − κ2exp{−νl + κl cos(θ − µ)}, (θ, l) ∈ [−π, π) × [0, ∞),(1)
fJW (θ, l; µ, ν, κ) =2π
where 0 ≤ κ < ν and −π ≤ µ < π. The distribution is obtained by maximizing the (Shannon)
entropy with E(L), E(L cos Θ) and E(L sin Θ) constant (see, Mardia, 1975, p. 352). The
functional form of the density indicates that as the length part (= l) increases, the concentration
around the location (= µ) also increases. Clearly, the mode of the density always occurs at l = 0.
As a more flexible symmetric cylindrical distribution, we show a special case of the density
in Abe and Ley (2017), whose joint pdf is given byfW M (θ, l; µ, κ, α, β) =αβ αlα−1 exp [−(βl)α {1 − tanh(κ) cos(θ − µ)}] ,2π cosh(κ)(2)

where (θ, l) ∈ [−π, π) × [0, ∞), α > 0 and β > 0 are the linear scale and shape parameters, and−π ≤ µ < π and κ ≥ 0 are the circular location and concentration parameters, respectively. We-42-

call this cylindrical distribution the Weibull–von Mises distribution. The density (2) is basically
proportional to the product of the Weibull and von Mises (replacing κ with κl) distributions,fW M (θ, l) ∝ lα−1 exp [−(βl)α ] × exp {(βl)α tanh(κ) cos(θ − µ)} ,divided by the normalizing constant.
As in Abe and Ley (2017), another cylindrical distribution is obtained by replacing the
Weibull distribution with the Gamma distribution. The pdf is given byfGM (θ, l; µ, κ, α, β) = C −1 lα−1 exp [−βl {1 − tanh(κ) cos(θ − µ)}] ,(3)

where (θ, l) ∈ [−π, π) × [0, ∞), α, β, γ, κ > 0, and −π ≤ µ < π. The normalizing constant C isgiven by∫ π∫ ∞
2πΓ(α)(cosh(κ))α Pα−1 (cosh(κ)).
lα−1 exp[−βl{1 − tanh(κ) cos(θ − µ)}]dldθ =βα−π 0We call this distribution the Gamma–von Mises distribution.

2.1

Parameter estimation

The most preferable orientation of asymmetry for focal tree i is written as


ni
∑
−
→
→
µi (a) = arg a S +CI ji −x−j xi  ,

(4)

j=1

−
→
where ni is the number of competitors, the summation covers all the competitors, and S =(0, −1). The parameter a quantifies the relative importance of solar radiation and competitorswith a given size and distance. By substituting µ = µi (a) into the symmetric cylindrical distributions and letting (θ1 , l1 ), . . . , (θn , ln ) be a sample of n independent and identically distributed
couples of circular and linear observations drawn from a cylindrical distribution. Then, the likelihood function can be expressed, in the cases of the basic three symmetric distributions, namely
the Johnson–Wehrly, Weibull–von Mises, and Gamma–von Mises distributions, respectively as
LJW (µ, κ, ν) =LW M (µ, κ, α, β) =

n
∏
i=1
n
∏

fJW (li , θi ; µ, κ, ν),(5)fW M (li , θi ; µ, κ, α, β),(6)fGM (li , θi ; µ, κ, α, β).

(7)

i=1

LGM (µ, κ, α, β) =

n
∏
i=1

The likelihood functions of the other distributions can be written in the same way. In general, it
is diﬃcult to give closed-form expressions for the maximum likelihood estimates (MLEs); hence,
numerical methods should be used to find the solutions.

References[1] Abe, T. & Ley, C. (2017). A tractable, parsimonious and flexible model for cylindrical
data, with applications. Econometrics and Statistics, 4, 91–104.

-43-