科研費シンポジウム「多様な分野における統計科学に関する理論と方法論の革新的展開」（Innovative development of theory and methodology on statisticalscience in various ﬁelds）
日時：2021 年 9 月 3 日（金）〜9 月 5 日（日）(Date: Friday, 3, – Sunday, 5, September, 2021)
場所：新潟大学駅南キャンパスときめいと 講義室 A,B(Place: Niigata University, Station South Campus Tokimate, Lecture
Room A, B)（TEL: 025-248-8141）
科学研究費・基盤研究（A) （課題番号：20H00576）「大規模複雑データの理論と方法論の革新的展開」
研究代表者：青嶋 誠（筑波大学）
Makoto Aoshima (University of Tsukuba)
開催責任者：蛭川 潤一（新潟大学）
Junichi Hirukawa (Niigata Universiry)Program
Friday, 3, September
Reception 13:00-13:15
Opening 13:15-13:20 Junichi HIRUKAWA (Niigata University)
Afternoon Session I (in Japanese) 13:20-14:45
Chair 宮田 庸一 (高崎経済大学)
1. 13:20-14:00 太刀岡 勇気
デンソーアイティーラボラトリ
深層学習モデルの内部状態クラスタリングへのノイズ掃き出し法とクロスデータ行列法の利用

2. 14:05-14:45 劉慶豊（Qingfeng Liu）
小樽商科大学（Otaru University of Commerce）

Machine CollaborationCoﬀee Break 14:45-15:00
Afternoon Session II (in English) 15:00-17:10
Chair Hiroshi Shiraishi (Keio University)
3. 15:00-15:40 張元宗（Chang Yuan-Tsung）
目白大学（Mejiro University）

Simultaneous estimation of multiplicative Poisson means in two-way contingency tables
4. 15:45-16:25 栗栖 大輔（Daisuke Kurisu）
東京工業大学（Tokyo Institute of Technology）
非定常な関数時系列データの統計分析（Statistical analysis of nonstationary functional time series）

5. 16:30-17:10 明石 郁哉 (Fumiya Akashi)
Faculty of Economics, University of Tokyo
Statistical inference for nonstationary heavy-tailed time series models by L1 approach

Saturday, 4, September
Morning Session (in Japanese) 9:20-11:30
Chair Muneya Matsui (Nanzan University)
6. 9:20-10:00 桃﨑 智隆 1 ・長 光司 1 ・中川 智之 2 ・富澤 貞男 21

東京理科大学大学院 理工学研究科

2

東京理科大学 理工学部

ベイズ法を用いた分割表における尺度の推定

7. 10:05-10:45 永井 勇
中京大学 教養教育研究院
高次元小標本データにおける多変量線形回帰モデルでの罰則を用いない推定法

8. 10:50-11:30 阿部俊弘
法政大学 経済学部

Cauchy 型分布の単純な EM アルゴリズムについて
Lunch 11:30-13:20
Afternoon Session I (in Japanese) 13:20-14:00
Chair 阿部俊弘 (法政大学)
9. 13:20-14:00 松井 宗也 (Muneya Matsui)
南山大学 (Nanzan University)
ランダムフィールドにおける距離の共分散 (Distance covariance for random ﬁelds)

Coﬀee Break 14:00-15:00

Afternoon Session II (in English) 15:00-16:25
Chair Fumiya Akashi (University of Tokyo)
10. 15:00-15:40 澁木 涼太郎（Shibuki Ryotaro）, 中村 知繁（Nakamura Tomoshige）, 白石 博
（Shiraishi Hiroshi）
慶應義塾大学理工学研究科 (Graduate School of Science and Technology, Keio University),
慶應義塾大学理工学部 (Faculty of Science and Technology, Keio University),
慶應義塾大学理工学部 (Faculty of Science and Technology, Keio University)
ランダムフォレストを用いた時系列分位点回帰 (Time Series Quantile Regressions by using Random

Forest)11. 15:45-16:25 竹内 努（Tsutomu T. TAKEUCHI）
名古屋大学素粒子宇宙物理学専攻（Division of Particle and Astrophysical Science, Nagoya University, Japan）
高次元統計学の方法による銀河の分光マップの解析（Analysis of Spatially Resolved Spectral Map of a Galaxy by the High-Dimensional Statistical

Method）
Coﬀee Break 16:25-16:40
Afternoon Session III (Guest speakers session) 16:40-18:10
Chair Junichi Hirukawa (Niigata University)
12. 16:40-17:20 Ngai Hang Chan
The Chinese University of Hong Kong
Statistical Inference for Glaucoma Detection
13. 17:30-18:10 Konstantinos Fokianos
University of Cyprus
Mixtures of Nonlinear Poisson Autoregressions

Sunday, 5, September
Morning Session (in Japanese) 10:00-11:25
Chair 阿部俊弘 (法政大学)
14. 10:00-10:40 宮田 庸一
高崎経済大学 経済学部
ブリッジ推定量におけるチューニングパラメーターの選択について

15. 10:45-11:25 種市 信裕・関谷 祐里・外山 淳
北海道教育大学（札幌校）
・北海道教育大学（釧路校）
・数学利用研究所
3 次元分割表における条件付き独立性検定の変換統計量について
Closing 11:25-11:30 Junichi HIRUKAWA (Niigata University)

深層学習モデルの内部状態クラスタリングへの
ノイズ掃き出し法とクロスデータ行列法の利用
太刀岡 勇気 (デンソーアイティーラボラトリ)

1 はじめに複数の話者やスタイルを制御できる end-to-end 音声合成システムで、タグづけにより話者・スタイルを制御する方法を提案した [1]。ところで統計学の分野では、データの次元数 d と標本数 n の間に d ≫ n の関係がある場合、高次元小標本 (high-dimension, low-sample-size; HDSS) データと呼ばれ、特殊な取り扱いが必要である [2]。埋め込み表現を主成分分析 (PCA) や t-SNE により次元圧縮して可視化することが、埋め込み表
現は 29 程度の次元、話者・スタイル数は 22 ∼ 24 程度であるため、この問題は HDSS のクラスタリング問題
であると考えられる。ここでは、音声合成器内部の埋め込み表現に対して、PCA と、HDSS データを扱う代
表的な 2 手法であるノイズ掃き出し法 (NRM)[3] とクロスデータ行列法 (CDM)[4] での分析を比較した。

2 実験
学習話者 ftak,fhar の 2 人の女性話者について、4(= n) スタイルでの内部表現を PCA した結果を、図 1 に示す。話者によって分離できていることと、h,t,z,k の順に上から並んでいることから、スタイルが同じものは近くに配置されていることが推測される。NRM による結果を図 2 に示す。ラベルによる並び順は異なるがほぼ同傾向といえる。CDM による結果を図 3 に示す。話者毎に集合分割すると、スタイル別にクラスタリング
できることがわかった。
スパイク指数の推定結果を図 4 に示す。点推定の結果、その平均値、中央値を示す。外れ値もあるがそれなりに収束している。下段にはべき乗近似により得られた αs を示している。これより αs > 0.8 であり、γ > 0.2
であれば標本固有値が一致性を持つことがわかり、ここでのは指数 γ = 2/9 ≃ 0.22 であるから問題ない範囲であるといえる。ただし一般化スパイクモデルが仮定している αs の s に関する単調減少性は満たされていない。
3

100

6

ftakh

90

1

5

fharh
fhart

80
70

4

0

ftakz

-1
ftakk

fharz

60
3

50
40

2

30

-2

20

1

fhark

-3

10

-4

0
-4

-2
0
2
1st Principal Component Score

Variance Explained (%)

ftakt

Eigen Value

2nd Principal Component Score

2

4

0
1

2

3 4 5 6 7
Principal Component

8

図 1 Principal component analysis (PCA) of embedded vectors.

1fhark

2.5

90

3.5

80
fharz

1

ftak

3
k

70

0.5
ftakz

0
-0.5 fhart
fharh

2.5

60
50

2

40

1.5

-1

30

Variance Explained (%)

1.5

Eigen Value

2nd Principal Component Score

100

4

2

1
ftakt

-1.5

20
0.5

-2
ftakh

-2.5
-2

10

0

0
2
1st Principal Component Score

0

4

1

2
3
4
5
Principal Component

6

3
fharz
ftakz

1

ftakh
fharh

0

-1

ftak
fhark
fhart
ftak

-2

100

4

90

3.5

80
70

3
Eigen Value

2nd Principal Component Score

2

4.5

60
2.5
50
2
40
1.5

k

t

-3

30

1

20

0.5

10

0
-2

0
2
1st Principal Component Score

Variance Explained (%)

図 2 Noise reduction methodology.

0

4

1
2
3
Principal Component

図 3 Cross-data-matrix methodology.
20
raw
mean
median

15
10
5
0
1

2

3

4

5

6

4

5

6

s

1.6
1.4
1.2
1
0.8
1

2

3
s

図 4 Estimated α (embedded vector, style).

3 まとめ
合成音声器の内部状態を分析する際の問題が HDSS の設定であることに着目し、NRM と CDM を適用した。NRM は通常の PCA よりも若干固有値が割り引かれるものの結果は大差がなかった。CDM は初期の集合の分割によって結果が大きく変わることがわかった。

参考文献[1] 太刀岡勇気，“話者タグによる tacotron2 の話者性制御，” 日本音響学会研究発表会講演論文集 (秋季)，pp.857–858 (2020)．
[2] 青嶋誠，矢田和善，高次元の統計学，共立出版 (2019)．
[3] K. Yata and M. Aoshima, “PCA consistency for the power spiked model in high-dimensional settings,” Journal ofMultivariate Analysis, 122, 334–354 (2013).
[4] K. Yata and M. Aoshima, “High-dimensional inference on covariance structures via the extended cross-data-matrix
methodology,” Journal of Multivariate Analysis, 151, 151–166 (2016).

2

Machine Collaboration∗
Qingfeng Liu† and Yang Feng‡
August 2, 2021
Ensemble learning has emerged and been extensively studied by many in the past few decades.
In general, the idea of ensemble learning is to combine the predictions obtained from different learning methods (hereafter, base machines), or predictions based on different subsamples to improveprediction performance. Bagging, stacking, and boosting are three prominent examples. In bagging(Breiman, 1996)/stacking, base machines first run in parallel and independently, and then the final
prediction is constructed as a simple/weighted average of the predictions from these base machines.
In boosting (Schapire et al., 1998), the base machines work jointly in a top-down manner. In the
aforementioned algorithms, the output from each base machine is fixed after being calculated. Like
human collaboration, an idea that may yield potential improvement is to let the base machines communicate with each other and update their outputs after observing the predictions of the other basemachines. Based on this idea, we propose the Machine Collaboration or MaC learning framework.
Compared with bagging, stacking, and boosting, MaC has the following desirable features. Figure
1 provides the schematic for bagging, stacking, boosting, and MaC. As illustrated, bagging andstacking are parallel & independent, boosting is sequential & top-down, while MaC is circular &interactive. For MaC, pieces of information are passed repeatedly between base machines around a“round table,“ but not one-way or top-down. In this type of scheme, the base machines update theirstructures and/or parameters according to the information received from the other machines. We
demonstrate that MaC can deliver competitive performance when compared with the base machinesor the other ensemble methods.

M1

M1

M2

M2
T

···

MK

Bagging/Stacking

MK

T

MK

+

to
up

M
MK
4
to

+

M2

les

···

up

···

les

···

+

cy
c

+

cy
c

les

· · 3·M

MK

Boosting

M1
T

M2

cy
c

M2

+
to

+

up

M1

Circular & Interactive

les

M1

Sequential & Top-downcycParallel & Independent

up

to

T

Collaboration

Figure 1: Bagging, boosting, and machine collaboration
The main contributions of this work are fourfold. First, we propose a new type of ensemble∗ The authors gratefully acknowledge the support of the Japan Society for the Promotion of Science through
KAKENHI Grant No. JP19K01582 (Liu), the Nomura Foundation for Social Science Grant No. N21-3-E30-010 (Liu)and a National Science Foundation CA-REER Grant No. DMS-1554804 (Feng).
† Corresponding author. Department of Economics, Otaru University of Commerce, Otaru City, Hokkaido, Japan.
Email: qliu@res.otaru-uc.ac.jp.
‡ School of Global Public Health at New York University, NY, NY, USA. Email: yang.feng@nyu.edu.

1learning framework, MaC, which is circular & interactive. The circular & interactive aspect could
be a potential direction for exploring new methods of ensemble learning. Second, we present somedesirable finite statistical properties of MaC. Third, we demonstrate via extensive simulations that
MaC performs better than all individual base machines and the ensemble methods SL and LSBoost. Lastly, in the analysis of real data, we compare MaC with the competing methods on 119
benchmark datasets in the Penn Machine Learning Benchmarks (PMLB) (Olson et al., 2017) for
evaluating and comparing machine learning algorithms. The results of this analysis demonstrate thenotable advantages of MaC for most datasets.

References
Breiman, L., 1996. Bagging predictors. Machine learning 24, 123–140.
Breiman, L., Friedman, J., Stone, C., Olshen, R., 1984. Classification and Regression Trees. The
Wadsworth and Brooks-Cole statistics-probability series, Taylor & Francis.
Cortes, C., Vapnik, V., 1995. Support-vector networks. Machine learning 20, 273–297.
Dasarathy, B.V., Sheela, B.V., 1979. A composite classifier system design: Concepts and methodology. Proceedings of the IEEE 67, 708–713.
Dong, X., Yu, Z., Cao, W., Shi, Y., Ma, Q., 2020. A survey on ensemble learning. Frontiers ofComputer Science 14, 241–258.
Friedman, J.H., 2001. Greedy function approximation: A gradient boosting machine. The Annalsof Statistics 29, 1189 – 1232.
Hastie, T., Tibshirani, R., Friedman, J., 2009. The elements of statistical learning: data mining,inference, and prediction. Springer Science & Business Media.
Ho, T.K., 1995. Random decision forests, in: Proceedings of 3rd international conference on document analysis and recognition, IEEE. pp. 278–282.
van der Laan, M.J., Polley, E.C., Hubbard, A.E., 2007. Super learner. Statistical Applications in
Genetics and Molecular Biology 6, 1544–6115.
Lu, X., Van Roy, B., 2017. Ensemble sampling, in: Guyon, I., Luxburg, U.V., Bengio, S., Wallach,
H., Fergus, R., Vishwanathan, S., Garnett, R. (Eds.), Advances in Neural Information Processing
Systems, Curran Associates, Inc.. pp. 3259–3267.
Mendes-Moreira, J., Soares, C., Jorge, A.M., Sousa, J.F.D., 2012. Ensemble approaches for regression: A survey. Acm computing surveys (csur) 45, 1–40.
Olson, R.S., La Cava, W., Orzechowski, P., Urbanowicz, R.J., Moore, J.H., 2017. Pmlb: a large
benchmark suite for machine learning evaluation and comparison. BioData Mining 10, 36.
Qi, Y., Liu, B., Wang, Y., Pan, G., 2019. Dynamic ensemble modeling approach to nonstationary
neural decoding in brain-computer interfaces, in: Wallach, H., Larochelle, H., Beygelzimer, A.,d'Alché-Buc, F., Fox, E., Garnett, R. (Eds.), Advances in Neural Information Processing Systems,
Curran Associates, Inc.. pp. 6087–6096.
Sagi, O., Rokach, L., 2018. Ensemble learning: A survey. Wiley Interdisciplinary Reviews: Data
Mining and Knowledge Discovery 8, e1249.
Schapire, R.E., 1990. The strength of weak learnability. Machine Learning 5, 197–227.
Schapire, R.E., Freund, Y., Bartlett, P., Lee, W.S., 1998. Boosting the margin: A new explanation
for the effectiveness of voting methods. The annals of statistics 26, 1651–1686.
Yu, Z., Luo, P., Liu, J., Wong, H.S., You, J., Han, G., Zhang, J., 2018. Semi-supervised ensemble
clustering based on selected constraint projection. IEEE Transactions on Knowledge and Data
Engineering 30, 2394–2407.

2

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

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

1

Introduction

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

J
X
j=1

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

i=1

I X
J
X

xij .

i=1 j=1

As discussed in Hara and Takemura (2006) complete sufficient statistics are x1 = (x1+ , . . . , xI+ ) and x2 =(x+1 , . . . , x+J ). The MLE of λij is( xi+ x+j
if x++ 6= 0
L
x++
λ̂M
=
ij
0
if x++ = 0.
They have given a class of improved estimators which shrink the MLE toward the origin under the normalized
squared error loss. The simple one is)(xi+ x+jdHTδij =, i = 1, . . . , I, j = 1, . . . , J,1−x++x++ + dNext section we consider one-dimensional shrinkage to an order statistic or a specified point.

2

One-dimensional shrinkage to an order statistic or a specified pointFirst, we consider one-dimensional shrinkage to an order statistic.
Let x(`)+ be the `-th smallest observation among x1+ , . . . , xI+ . We assume that I ≥ ` + 2 and consider thefollowing estimator which shrinks xi+ toward x(`)+ when xi+ ≥ x(`)+ :()(xi+ − x(`)+ )+x+j(1)
xi+ − ϕ(W ), i = 1, . . . , I, j = 1, . . . , J,
δij =
x++
W +d
where W =PI++
i=1 (xi+ − x(`)+ ) , a = max(0, a) and d is a positive constant. Then we have the following.

Theorem 2.1. Suppose that ϕ(W ) is a non-decreasing function satisfying 0 ≤ ϕ(W ) ≤ 2(I − ` − 1) and that d ≥PI(1)L
sup ϕ(W )/2. Then δij , i = 1, . . . , I improves upon the MLE λM
ij , i = 1, . . . , I under the loss function
i=1 (λ̂ij −2λij ) /λij for any j = 1, . . . , J.
L
Next we consider the estimators shrink λ̂M
to a specified non-negative values, bi .
ij

1Let bi ≥ 0, i = 1, . . . , I be given numbers and we propose the following shrinkage estimator which shrinks xi+
to bi when xi+ ≥ bi :()
x+j(xi+ − bi )+
(2)
δij =
xi+ − ϕ(N, W ), i = 1, . . . , I, j = 1, . . . , J,x++
W + d(N )where W =PI+
i=1 (xi+ − bi ) and N = #{i|xi+ ≥ bi }. Then we have the following.

Theorem 2.2. Suppose that ϕ(N, W ) is a non-decreasing function of W and satisfies 0 ≤ ϕ(N, W ) ≤ 2(N − 1)+(2)
for any 0 ≤ N ≤ I. Suppose that d(N ) ≥ supW ϕ(N, W )/2. Then δij , i = 1, . . . , I improves upon the MLE
P
I
L
2
λ̂M
ij , i = 1, . . . , I under the loss function
i=1 (λ̂ij − λij ) /λij for any j = 1, . . . , J.
It may be noticed that the shrinkage is made only when N ≥ 2.
Remark Theorems 2.1 and 2.2 can be generalized directly to the case of Poisson multiplicative model for a
multi-way contingency tables.
Next section we also consider two-dimensional shrinkage to order statistics or to a specified point.

3

Two-dimensional shrinkage to order statistics.

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

1
x++



(xi+ − x(`)+ )+xi+ − ϕ1 (W1 )W1 + d1

(x+j − x+(m) )+x+j − ϕ2 (W2 ),W2 + d 2
i = 1, . . . , I, j = 1, . . . , J,
(2.4)



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

1
x++

(

(xi+ − bi )+xi+ − ϕ1 (N1 , W1 )W1 + d1 (N1 ))()(x+j − cj )+x+j − ϕ2 (N2 , W2 ),W2 + d2 (N2 )i = 1, . . . , I, j = 1, . . . , J,

where W1 =

PI

PJ

+
j=1 (x+j − cj ) , N1 = #{i|xi+ ≥ bi , i = 1, . . . , I} and N2 = #{j|x+j ≥PJPI
Although it may be natural to put the condition i=1 bi = j=1 cj , we do not need it in the+
i=1 (xi+ − bi ) , W2 =ci , j = 1, . . . , J}.
following.

(2.5)

Theorem 3.2. Suppose that ϕi (Ni , Wi ) is a non-decreasing function of Wi and satisfies 0 ≤ ϕi (Ni , Wi ) ≤ (Ni −1)+
(4)
for any Ni ≥ 0, and that di (Ni ) ≥ (Ni − 1)+ /Ni supWi ϕi (Ni , Wi ), for any Ni ≥ 0, i = 1, 2. Then δij improves
P
P
I
J
L
upon the MLE λ̂Munder the loss function i=1 j=1 (λ̂ij − λij )2 /λij .
ij
It may be noticed that the shrinkage in the i−th dimension is made only when Ni ≥ 2.

2

Statistical analysis of nonstationary functional time series
Daisuke Kurisu
Tokyo Institute of Technology

1 Introduction
In this study, we develop an asymptotic theory for estimating the time-varying characteristics of a
locally stationary functional time series. The notion of a locally stationary functional time series is anextension of that of a locally stationary process introduced by Dahlhaus (1997). For some probability
space (Ω, A, P ) and for p ≥ 1, let Lp (Ω, A, P ) denote the space of real-valued random variables such that∥X∥p = (E[|X|p ])1/p < ∞. Let D ⊂ Rd be a compact set and let L2 = L2 (D) denote a Hilbert spacewith inner product ⟨·, ·⟩ defined by∫⟨x, y⟩ =x(t)y(t)dt, x, y ∈ H = L2 .
D

Further, let LpH = LpH (Ω, A, P ) denote the space of H-valued random functions X such thatp

1/p

(E[∥X∥ ])( [(∫)p/2 ])1/p
2= E
X (t)dt< ∞.

2 Settings
2.1 local stationarity
The definition of an H-valued locally stationary process is as follows: The H-valued stochastic process{Xt,T } in LpH is locally stationary if for each rescaled time point u ∈ [0, 1], there exists an associated(u)H-valued process {Xt } in LpH with the following properties:(u)(i) {Xt }t∈Z is strictly stationary.
(ii) It holds that((u)
∥Xt,T − Xt ∥ ≤

1
t
−u +
T
T

)
(u)

Ut,T a.s.,

(u)

(u)

for all 1 ≤ t ≤ T where {Ut,T } is a process of positive variables satisfying E[(Ut,T )ρ ] < C for someρ > 0, C < ∞ that is independent of u, t, and T .
This definition is a natural extension of the notion of local stationarity for real-valued time series introduced in Dahlhaus (1997).

2.2 dependence structure
(u)

(u)

For each u ∈ [0, 1], we assume that {Xt } is Lp -m-approximable, that is, Xt(u)Xt∈ LpH is of the form

= fu (εt , εt−1 , εt−2 , . . . ),where εj are i.i.d. elements taking values in a measurable space S, and fu is a measurable function
fu : S ∞ → H. Note that Xt

(u)

(m)

is strictly stationary. We also assume that if {εj,k } is an independentcopy of {εj } defined on the same probability space, then letting
(m)

(u)

(m)

Xm,t = fu (εt , εt−1 , . . . , εt−m+1 , εt,t−m , εt,t−m−1 , . . . ),we have∑(u)

vp (Xt(u)− Xm,t ) < ∞.

m≥1

We can show that a wide class of functional time series (e.g., functional AR(1) process and functional
ARCH(1) process) is Lp -m-approximable under some regularity conditions. See Hörmann and Kokoszka(2010) and Horváth and Kokoszka (2012) for details on the properties of Lp -m-approximable randomfunctions.

3 Results
We introduce a kernel-based method to estimate the time-varying covariance operator and the timevarying mean function of a locally stationary functional time series. Subsequently, we derive the convergence rate of the kernel estimator of the covariance operator and associated eigenvalue and eigenfunctions.
We also establish a central limit theorem for the kernel-based locally weighted sample mean. As applications of our results, we discuss the prediction of locally stationary functional time series and methods fortesting the equality of time-varying mean functions in two functional samples. We also discuss a problem
to estimate the number of principal components to be used. To the best of our knowledge, this is the first
paper that develops an asymptotic theory of estimating time-varying characteristics of locally stationary
functional time series based on kernel methods.

References
Dahlhaus, R. (1997). Fitting time series models to nonstationary processes. Ann. Statist. 25 1-37.
Hörmann, S. and Kokoszka, P. (2010). Weakly dependent functional data. Ann. Statist. 38 1845-1884.
Horváth, L. and Kokoszka, P. (2012). Inference for Functional Data with Applications. Springer, NewYork.
Kurisu, D. (2021). On the estimation of locally stationary functional time series. arXiv2105.11873.

2

Statistical inference for nonstationary heavy-tailed
time series models by L1 approach
Fumiya Akashi (University of Tokyo)∗

Abstract
This talk proposes a robust L1 -estimation methods for a time-varying autoregressive (AR) model. First, we construct a robust local linear estimator based on
the self-weighting approach (Ling (2005, Journal of the Royal Statistical Society,
Series B)), and show the asymptotic normality of the proposed estimator. Second,
the generalized empirical likelihood statistic is proposed to test the hypothesis of
the coeﬃcients of the model. The simulation experiments are given in the talk.

1

The model and self-weighted statistics

Suppose that an observed stretch {Y1,T , ..., Yt,T } is generated by the following time-varying
autoregressive (AR) model
Yt,T = β(t/T )⊤ Xt−j,T + ϵt ,(1)

where β(u) = (β1 (u), ..., βp (u))⊤ is a function [0, 1] → Rp , Xt−1,T is a p-dimensional vector
deﬁned as Xt−1,T := (Yt−1,T , ..., Yt−p,T )⊤ , and {ϵt : t ∈ Z} is a sequence of i.i.d. zero-medianrandom variables. We assume some conditions for {ϵt : t ∈ Z}. In particular, the taildistribution of ϵt is assumed to satisfy xα P (ϵt > x) → qC and xα P (ϵt < −x) → (1 − q)C
as x → ∞ with some α > 0, q ∈ [0, 1] and C ≥ 0. In particular, ϵt do not have the inﬁnite
variance when α < 2.
We deﬁne the self-weighed local linear L1 estimator for β(u0 ) (u0 ∈ [0, 1]) as(β̂(u0 ), γ̂) := arg min
β,γ

∗

T
∑

(
K

t=p+1

t − t0Th)
wt−1,T |et−1,T (β, γ)|,This talk is partially based on a joint work with Junichi Hirukawa (Niigata University) and Kon-

stantinos Fokianos (University of Cyprus).

1where K is a kernel function, h is a bandwidth parameter which goes to zero as T → ∞,t0 is an integer satisfying |t0 /T − u0 | < 1/T , γ is a nuisance parameter and
[
(
) ]⊤
t
et−1,T (β, γ) := Yt,T − β +− u0 γ Xt−1,T (β, γ ∈ Rp ).
T
In particular, wt−1,T (t = p + 1, ..., T ) is called the self-weight, which is a measurable
function of Xt−1,T . The self-weighting approach for time series model was originallyproposed by Ling (2005) for a stationary AR-model.
The limit distribution of the proposed estimator is given in the following theorem.
Theorem 1. Under some regularity conditions, we have
(
)
]
√ [
κ2
d
−1
−1
→ N 0p ,Σ(u0 ) Ω(u0 )Σ(u0 )T h β̂(u0 ) − β(u0 ) −4f (0)2(T → ∞),

where Σ(u0 ) and Ω(u0 ) are nonsingular matrices.
Another important issue in statistical inference for inﬁnite variance process is hypothesis testing for the coeﬃcients of the model (1). Let us consider the hypothesis
H : β(u) = β0 (u) (∀u ∈ [0, 1]), and we apply the generalized empirical likelihood (GEL)
test statistic. Let us deﬁne the self-weighted moment function()
gt,T (b) := wt−1,T sign Yt,T − b⊤ Xt−1,T Xt−1,T (b ∈ Rp )and the self-weighted GEL test statistic asrT (β(·)) := 2

sup

T
∑

{
}
ρ λ⊤ gt,T (β(t/T )) ,λ∈Λ̂T (β(·)) t=p+1

where ρ is a score function whose domain is Vρ (⊂ R), and{}Λ̂T (β(·)) := λ : λ ∈ Rp and λ⊤ gt,T (β(t/T )) ∈ Vρ for all t ∈ {1, ..., T } .
Then, we get the following theorem.
p

Theorem 2. Under some regularity conditions, we have rT (β0 (·)) −→ χ2p as T → ∞.
Remarkably, the proposed statistics always converge to pivotal limit distributions regardless of whether the model has ﬁnite or inﬁnite variance. Finite sample performanceof the proposed statistics is also illustrated by some simulation experiments.

References
Ling, S. (2005). Self-weighted least absolute deviation estimation for inﬁnite variance
autoregressive models. Journal of the Royal Statistical Society. Series B: Statistical
Methodology 67 (3), 381–393.
2

ベイズ法を用いた分割表における尺度の推定
桃﨑 智隆 1

長 光司 1
1

中川 智之 2

富澤 貞男 2

東京理科大学大学院 理工学研究科
2

東京理科大学 理工学部

2 元分割表の解析では, 行変数と列変数の間の連関尺度を用いて独立性からの隔たりの程度を測る. 連関尺度には, 例えば, Cramér 係数 (Cramér, 1946) などがある. 一方で, 行と列が同じ分類からなる正方分割表の解析では, 行と列の変数間の独立性に代わり対称性からの隔たりの程度を測る. Tomizawa et al. (1998) は power

divergence (または, diversity index) を用いて, 対称性からの隔たりの程度を測る様々な尺度を提案した.
分割表におけるこれらの尺度は分割表のセル確率の関数として表されるので, 尺度の値は未知の値を取る.
そのため尺度の推定には, 標本比率を用いた尺度のプラグイン推定量が用いられる. サンプル数が十分にあるとき, 標本比率のプラグイン推定量は近似不偏推定量となる. しかし, サンプル数が十分にない場合は, 推定量のバイアスや平均二乗誤差 (MSE) が大きくなってしまう.

Tomizawa et al. (2007) は, バイアスの高次のオーダーを導出し, バイアス補正を行うことで尺度の推定量を改良した. これらの研究では数値実験を用いて, 改良型推定量は標本比率のプラグイン推定量よりも速く尺度の真値に近づくことが示された. しかしながら, 改良型推定量にはいくつかの問題点がある. 確かに, サンプ
ル数が小さい場合でも改良型推定量はバイアスを小さくできているが, MSE も小さくできるとは限らない. さらに, 改良型推定量の値域と尺度の値域が異なる可能性もある. 例えば, Tomizawa et al. (1998) の尺度の値域
は 0 以上 1 以下であるが, 改良型推定量の値域は 0 以上 1 以下を超えている. また, このような場合, 改良型推定量を用いた尺度の推定量の信頼区間の解釈は容易ではない.
本講演では, ベイズ法を用いて, 十分なサンプルサイズがなくてもバイアスや MSE を小さくできる尺度の推定量を提案した.

2 元 r × c 分割表を考える. 確率変数ベクトル n = (n11 , n12 , . . . , n1c , n21 , . . . , nrc )⊤ は多項分布 M (n, p)∑に従うとする. ここで, n = i,j nij , pij は (i, j) セル確率, p = (p11 , p12 , . . . , p1c , p21 , . . . , prc )⊤ , ⊤ は転置を表す.
セル確率ベクトル p がディリクレ事前分布p(p|α) =Γ(rcα) ∏ α−1p(Γ(α))rc i,j ijに従うとすると, p の事後平均は以下のようになる.
(α) ⊤p̂(α) = (p̂11 , p̂12 , . . . , p̂1c , p̂21 , . . . , p̂rc
)
(α)

(α)

ただし
(α)

p̂ij =(α)(α)

nij + αn + rcαであり, Γ(·) はガンマ関数である. ここで, α = 0 のときは, 事後平均 p̂(α) は標本比率 p̂ = n−1 n に対応している.

関数 f (·) を分割表の尺度とする. ここで, f (·) は p において少なくとも 4 回微分可能であるとする. 尺度

f (p) の値は, f (p) において p を p̂ に置き換えた推定量 f (p̂) を用いて推定できる. また, f (p) において p をp̂(α) に置き換えた推定量 f (p̂(α) ) で推定することも可能である. 本講演では, ディリクレパラメータ α を選択する方法の 1 つとして, f (p̂(α) ) の MSE を最小化する α を考えた. そのため, f (p̂(α) ) の MSE を漸近的に評価することで, f (p̂(α) ) の MSE を最小化する α を導出した. すなわち, 以下のようなディリクパラメータの導出を行った.

α̃ = argmin lim n2 MSE[f (p̂(α) )]n→∞α

推定量 f (p̂(α) ) の MSE は

MSE[f (p̂(α) )] =)1 (A1 α2 − 2A2 α + (the terms independent of α) + o(n−2 )2nとなる. ここで

[(

)(

)
]
∂f (p)⊤
A1 =tr(rcp − 1rc )(rcp − 1rc ) ,∂p⊤) [( 2]()∂f (p)∂ f (p)1⊤⊤
tr(diag(p) − pp )A2 = (rcp − 1rc )
2
∂p
∂p∂p⊤
)(
[(
)
]
∂f (p)∂f (p)⊤+ rc tr(diag(p) − pp )
∂p
∂p⊤
[(
)
( 2
)
]
∂f (p)∂ f (p)⊤⊤+ tr(rcp − 1rc )(diag(p) − pp ) ,∂p∂p∂p⊤∂f (p)∂p

1rc は全ての要素が 1 の rc × 1 ベクトル, diag(p) は主対角線の要素が p の対角行列である. よって, f (p̂(α) )の MSE を最小化するディリクレパラメータ α はα̃ = argmin lim n2 MSE[f (p̂(α) )] =
α

n→∞

A2
A1

(1)

となる. 以上より, 尺度 f (·) の推定量として f (p̂(α̃) ) を提案した.
本講演では, 数値実験を用いて提案推定量が多くの場面において標本比率のプラグイン推定量や改良型推定
量よりもバイアスや MSE を小さくできることを示した. さらに, ディリクレパラメータに α = 1(一様事前分布) や α = 1/2(Jeﬀreys 事前分布) を用いたり, Fienberg and Holland (1973) の方法で選択したディリクレパ
ラメータを用いたりした場合における事後平均のプラグイン推定量よりも, 提案手法によって選択されるディ
リクレパラメータを用いた場合における事後平均のプラグイン推定量の方が, バイアスや MSE を小さくできることも明らかになった. また, 導出したディリクレパラメータを用いることで, 尺度の信用区間をモンテカル
ロシミュレーションにより構成可能であることも示した.

参考文献
Cramér, H. (1946). Mathematical Methods of Statistics. Princeton, N.J., Princeton Univ. Press.
Fienberg, S. E. and Holland, P. W. (1973). Simultaneous estimation of multinomial cell probabilities.
Journal of the American Statistical Association, 68(343), 683–691.
Tomizawa, S., Miyamato, N., and Ohba, N. (2007). Improved approximate unbiased estimators of measureof asymmetry for square contingency tables. Advances and Applications in Statistics, 7(1), 47–63.
Tomizawa, S., Seo, T., and Yamamoto, H. (1998). Power-divergence-type measure of departure fromsymmetry for square contingency tables that have nominal categories. Journal of Applied Statistics,
25(3), 387–398.

Cauchy 型分布の単純な EM アルゴリズムについて
法政大学 経済学部
阿部 俊弘
対称分布の非対称化として, Azzalini (1985) に始まる skew-symmetric 分布が良く知られている. これの発展については, Azzalini (2005) や Azzalini & Capitanio (1999) にまとめられている.
Lin et al. (2007) では確率表現を活用して歪正規分布の EM アルゴリズムを与えている. この歪正規分布の多変量化はいくつのものが提案されているが, Sahu et al. (2003) の多変量歪正規分布についての EM アルゴリズムは Lin (2009) により与えられている.
Lin et al. (2007) や Lin (2009) による EM アルゴリズムでは, skew パラメータに対する推定式が陽に解けないため, 数値的手法が必要となる. これに対して, Abe et al. (2021) では, 歪正規分布の確率表現に overparameter を導入することにより, EM アルゴリズムの陽的表現を導出し, さらに, Chen et al. (2014) と比較し, 推定手法として優れていることを示した. 歪正規分布は様々な良い性質がある一方で, 歪正規分布の Fisher information matrix は skew パラメータが 0 の周りで,特異になる分布であることが指摘されている (Azzalini, 2013; Hallin & Ley, 2014).
この講演では, まず初めに対称な Cauchy 型分布の EM アルゴリズムを統一的な視点で与え, 関連したモデルについても EM アルゴリズムを与える: 対数 Cauchy 分布は裾の重い分布の一例であり, その名前から明らかなように, Z が Cauchy 分布に従うときに X = exp(Z) が対数 Cauchy 分
布に従う. 対数 Cauchy 分布は, Cauchy 分布と同様に一切の (非自明) モーメントが無限大になり,しばしば非常に裾の重い分布 (super-heavy tailed distribution) とされる. 対数 Cauchy 分布のパ
ラメータ推定に関して, 標本の自然対数をとったものの中央値は, µ のロバストな推定量になり, 標
本の自然対数をとったものの中央絶対偏差は σ のロバストな推定量になる. 対数 Cauchy 分布は,有意な外れ値または極値が発生するようなある種の生存過程や生物種の豊富度パターンのモデル化に用いることができる. 例としては, ヒト免疫不全ウイルスの感染から発症までの時間が挙げられる.
上記の Cauchy 分布と同様にして, 密度関数の形から t 分布の EM アルゴリズムも与えていく.
Liu & Rubin (1995) でも t 分布の確率表現から EM アルゴリズムの陽的表現を導出しており, 彼らのものとは別の観点のものであるが, M-step での解は本質的に彼らのものと同じものを導出する.
上で述べた歪正規分布にはいくつかの別形があるが, t 分布の歪分布化についても様々なものが提案されている. skew-symmetric 分布の pdf については有名なもので次のような表現がある (seeAzzalini, 2005): まず, G を x = 0 の周りで対称な分布の cdf とし, f0 を x = 0 の周りで対称な分布の pdf とする. このとき, skew-symmetric 分布の pdf

2G(λx)f0 (x)において x 7→

x−µ
とすると,
σ

) (
)
(
x−µ
2
x−µ
f0
G λ
σ
σ
σ

となる. Azzalini & Capitanio (2003) では, 歪正規分布に scale mixture をすることにより, pdf が
√
(
) (
)
ν+1
x−µ
2
x−µ
fACST (x; ν) = FT λ
; 1, ν + 1 fT
;ν
(1)
σ
σ
(x − µ)2 /σ 2 + νσ1である skew-t 分布を導出している. ここで, fT , FT はそれぞれ, 自由度 ν > 0 の t 分布の pdf と cdf()− ν+1
2x2)
fT (x; ν) = √1+,
ννB 2 , 21
1(ν1
FT (x; σ, ν) :=
σ

(

∫ x
−∞

fT

)
t
; ν dtσ(2)であり, B(·, ·) は Beta 関数である.

References
Abe, T., Fujisawa, H., Kawashima, T. & Ley, C. (2021). EM algorithm using overparameterization for multivariate skew-normal distribution. Econometrics and Statistics, 19, 151–168.
Azzalini, A. (1985). A class of distributions which includes the normal ones. Scandinavian Journal of Statistics 12(2), 171–178.
Azzalini, A. (2005). The skew-normal distribution and related multivariate families. Scandinavian Journal of Statistics 32 (2), 159–200.
Azzalini, A. (2013). The Skew-Normal and Related Families. vol. 3. Cambridge University Press.
Azzalini, A. & Capitanio, A. (1999). Statistical applications of the multivariate skew normal distribution. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 61(3),579–602.
Azzalini, A. & Capitanio, A. (2003). Distributions generated by perturbation of symmetry withemphasis on a multivariate skew t-distribution. Journal of the Royal Statistical Society: Series
B (Statistical Methodology), 65(2), 367–389.
Chen, L., Pourahmadi, M. & Maadooliat, M. (2014). Regularized multivariate regression models
with skew-t error distributions. Journal of Statistical Planning and Inference 149, 125–139.
Hallin, M. & Ley, C. (2014). Skew-symmetric distributions and Fisher information: the double
sin of the skew-normal. Bernoulli 20, 1432–1453.
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.
Liu, C. & Rubin, D. B. (1995). ML estimation of the t distribution using EM and its extensions,
ECM and ECME. Statistica Sinica, 19–39.
Ma, Y. & Genton, M. G. (2004). Flexible class of skew-symmetric distributions. Scandinavian
Journal of Statistics, 31(3), 459–468.
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.

2

Distance covariance for random ﬁelds
Nanzan University

Muneya Matsui

Abstract
We study an independence test based on distance correlation for random ﬁelds (X, Y).
We consider the situations when (X, Y) is observed on a lattice with equidistant grid
sizes and when (X, Y) is observed at random locations. We provide asymptotic theory
for the sample distance correlation in both situations and show bootstrap consistency.
The latter fact allows one to build a test for independence of X and Y based on the considered discretizations of these ﬁelds. We illustrate the performance of the bootstrap
test in a simulation study involving fractional Brownian and inﬁnite variance stableﬁelds. The independence test is applied to Japanese meteorological data, which areobserved over the entire area of Japan.

Introduction
It is well known that two q- and r-dimensional random vectors X and Y, respectively, are independent if and only if their joint characteristic function factorizes, i.e.,



 

φX,Y (s, t) = E exp(i s⊤ X + it⊤ Y) = E exp(i s⊤ X) E exp(i t⊤ Y) = φX (s) φY (t) , s ∈ Rq , t ∈ Rr .
Since this identity should hold with any (s, t) ∈ Rq × Rr , in general it is recommended to use aweighted L2 -distance between φX,Y and φX φY for the test of independence.
We speciﬁcally call the weighted L2 -distance as distance covariance between X and Y if theweight function is given asZ
2
φX,Y (s, t) − φX (s)φY (t) |s|−(q+β) |t|−(r+β) ds dt , β ∈ (0, 2),
T β (X, Y) = cq cr
Rq+r

R
where the constants cd for d ≥ 1 are chosen such that cd Rd (1 − cos(s′ x)) |x|−(d+β) dx = |s|β . The
quantity T β (X, Y) is ﬁnite under suitable moment conditions on X, Y. The corresponding distancecorrelation is given by
T β (X, Y)
Rβ (X, Y) = p.
pT β (X, X) T β (Y, Y)
Of course, X and Y are independent if and only if Rβ (X, Y) = T β (X, Y) = 0. Thanks to the choiceof the weight function, T β (X, Y) has an explicit form: assuming that (Xi , Yi ), i = 1, 2, . . . , are iid
copies of (X, Y), we haveT β (X, Y) = E[|X1 − X2 |β |Y1 − Y2 |β ] + E[|X1 − X2 |β ]E[|Y1 − Y2 |β ] − 2 E[|X1 − X2 |β |Y1 − Y3 |β ] ,and Rβ (cX, cY) = Rβ (X, Y) for c ∈ R, i.e., Rβ is scale-invariant.
With the ﬁnal form, the distance covariance for vectors is easily extended into that for random ﬁelds.
Let (X, Y) be a pair of random ﬁelds on B ⊂ Rd and deﬁne the norm ∥ f ∥2 =R

−1
2
|B| B f (u)du 1/2 . The distance covariance T β (X, Y), β ∈ (0, 2), between two random ﬁelds
X, Y on some bounded Borel set B ⊂ Rd of ﬁnite positive Lebesgue measure is deﬁned by


β
β
β 
β
T β (X, Y) = E ∥X1 − X2 ∥2 ∥Y1 − Y2 ∥2 + E ∥X1 − X2 ∥2 E ∥Y1 − Y2 ∥2ββ− 2 E ∥X1 − X2 ∥2 ∥Y1 − Y3 ∥2 ,where (Xi , Yi ), i = 1, 2, . . . , are iid copies of (X, Y), and the distance correlation Rβ (X, Y) is deﬁned
correspondingly. We apply the statistic Rβ (X, Y) to an independence for random ﬁelds.

Usually a sample of whole paths of (X, Y) is rarely at our disposal, and we observe (X, Y) on alattice with equidistant grid sizes or observe (X, Y) at random locations. Accordingly asymptotic
behavior of the test statistic depends on each of sampling schemes. We provide asymptotic theory
for the sample distance correlation in both situations and show bootstrap consistency. The latter
fact allows one to build a test for independence of X and Y based on the considered discretizationsof these ﬁelds. We illustrate the performance of the bootstrap test in a simulation study involving fractional Brownian and inﬁnite variance stable ﬁelds. The independence test is applied toJapanese meteorological data, which are observed over the entire area of Japan.

Main contents
1. Distance covariance for random ﬁelds on a lattice in [0, 1]d .
We provide asymptotic theory for the lattice case on B = [0, 1]d for increasing intensity p whichis the total (deterministic) number of grid points. We show that the discretized sample distancecorrelation converges to non-discretized whole path sample distance correlation as pn → ∞ depending on the sample size n → ∞. Therefor they have the same asymptotic distribution. We
show the bootstrap consistency for this sampling scheme.
2. Distance covariance for random ﬁelds at random locations.
We deﬁne T β (X (p) , Y (p) ) and Rβ (X (p) , Y (p) ) for non-lattice based discretizations X (p) , Y (p) of X, Y
on B. We choose a random number N p of random locations (Ui ) where the processes X, Y are oba.s.
served. Typically, these locations are uniformly distributed on B and N p → ∞ as p increases withthe sample size n to inﬁnity. We show that the sample distance correlation at random locations
well approximate the whole path sample distance correlation as pn → ∞. Moreover we showthat n T n,β (X (p) , Y (p) ) conditional on (N (p) ) p>0 has the same weak limit as n T n,β (X, Y). We provide
consistency of a suitable bootstrap procedure in this case.
3. A Monte Carlo study.
We conduct a Monte Carlo study of the ﬁnite sample behavior of the sample distance correlation
for β = 1. We consider a fractional Brownian sheet and, as a heavy-tailed alternative, we choose
symmetric 1.8-stable Lévy sheets. The observations are given either on a lattice or at random locations in [0, 1]2 . In addition, we illustrate the performance of the bootstrap procedure for the test
for independence based on distance correlation in the cases of ﬁxed locations on a lattice and of
randomly scattered locations. We focus on independent pairs X, Y, Brownian or 1.8-stable sheets.
4. An application to Japanese meteorological data.
We apply our results to Japanese meteorological data. We choose the 3 most fundamental factors: temperature (temp), precipitation (prec) and wind speed (wind) which have been observed
for a long time and over a wide range of Japan. We focus on the distance correlations for thepairs (prec & temp), (prec & wind), (temp & wind) and conduct the bootstrap tests for pair-wise
independence.

References[1] Dehling, H.G., Matsui, M., Mikosch, T., Samorodnitsky, G. and Tafakori, L. (2020) Distance covariance for discretized stochastic processes. Bernoulli 26, 2758–2789.
[2] Matsui, M., Mikosch, T. and Samorodnitsky, G. (2017) Distance covariance for stochastic processes.
Probab. Math. Statist. 37, 355–372.
[3] Matsui, M., Mikosch, T., Roozegar, R. and Tafakori, L. (2021) Distance covariance for random ﬁelds.
arXiv:2107.03162.

2

Time Series Quantile Regressions by using Random Forest(ランダムフォレストを用いた時系列分位点回帰)
Ryotaro Shibuki∗ , Tomoshige Nakamura† , Hiroshi Shiraishi‡

In this paper, we discuss an estimation procedure of conditional quantile by using random forests in timeseries setting. Our study is an extension of the quantile random forest (QRF) by Meinshausen (2006), thegenerarized random forest (GRF) by Athey et. al (2019) and random forest in time series setting by Davis andNielsen (2020).
Model

Let (εt )t≥1 be a sequence of i.i.d. random variables with E[εt ] = 0 and E[ε2t ] < ∞, and ﬁx an integerp ≥ 1. Given a measurable function g : Rp → R, deﬁne the process (Yt )t≥1 recursively by
Yt = g(Xt ) + εt ,Xt = (Yt−1 , . . . , Yt−p ).

(1)

In addition to the initial data η = (Y0 , Y−1 , · · · , Y1−p ), suppose that we have T observations Y1 , . . . , YT fromthe model (1) available and that we group them in input-output pairs, DT = {(X1 , Y1 ), . . . , (XT , YT )}. For
each ﬁxed value τ ∈ (0, 1), we seek forest-based (function) estimator of q∗τ : Rp → R deﬁned by a solution of
local estimating equation of the formΨτ (q∗τ , x) := E[ψqτ∗τ (Yt )|Xt = x] = 0,

for all x ∈ Rpt(2)

where ψqτ (y) = τ − 1{y≤q} . Let q∗τ ≡ (q∗τ (x))x∈Rp be the solution of (2) under the model (1). We assume that
there exists q∗τ (x) for all x ∈ Rp and τ ∈ (0, 1).
Double sample

We next deﬁne our random forests following Athey et. al (2019). Our random forests consistsof the double sample trees, which are regression trees based on two subsamples Is and Js from sample DT .
Deﬁnition 1. (Double Sample) Suppose that sample DT is available and the sub-sample size s = s(T ) withs ≤ T is provided. Letjskl s mI
I
J
J
I
J
As := A = A ∪ A ) ⊂ {1, 2, . . . , T } A ∩ A ) = ∅, A =
, A =
2
2
For any A = AI ∪ AJ ) ∈ As , we deﬁne two sub-samples Is and Js by Is = DAI , Js = DAJ where Dθ· ={(Xt , Yt )}t∈A· .
Splitting rule

We next deﬁne splitting rule in order to construct the double-sample regression trees following

Wager and Athey (2018).
Deﬁnition 2. (Splitting rule) Given subsample Js in Deﬁnition 1, we deﬁne a sequence of partitions P0 , P1 , . . .
by starting form P0 = {Rp } and then, for each ℓ ≥ 1, construct Pℓ from Pℓ−1 by replacing one set (parent node)P ∈ Pℓ−1 by (child node) C1 := {x = (x1 , . . . , xp ) ∈ P ⊂ Rp : xξ ≤ ζ} and C2 := {x = (x1 , . . . , xp ) ∈ P ⊂ Rp :xξ > ζ}, where the split direction ξ ∈ {1, . . . , p} is randomly chosen (i.e., random split) .
∗ Graduate School of Science and Technology, Keio University† Faculty of Science and Technology, Keio University‡ Faculty of Science and Technology, Keio University

1

Double-sample regression trees

A given partition Λ of Rp is called “recursive” if Λ = Pℓ for some ℓ ≥ 0,

where P0 , . . . , Pℓ are obtained as above. Note that the splitting rules determining how to choose node, direction
and position of a split may depend on the data DT , the double sampling procedure A ∈ As and the sequencei.i.d.

of independent random splittings ξ = {ξi }i=1,...,ℓ with ξi ∼ Ξ. By using the recursive partition Λ, we deﬁne
our double-sample regression trees.
Deﬁnition 3. Given a recursive partition Λ(A, ξ, DT ) = {L1 , . . . , L|Λ| } and a ﬁxed x ∈ Rd , q ∈ R, τ ∈ (0, 1),
our double-sample regression tree T (q, x; A, ξ, DT ) is deﬁned by
X 1{Xt ∈LA (x)}ψqτ (Yt )T (q, x; A, ξ, DT ) =|LA (x)|It∈A

where LA (x) = {Xt : Xt ∈ L̃(x)} ∩ Is and L̃(x) ∈ Λ(A, ξ, DT ) is a leaf containing x (i.e., x ∈ L̃A (x)).
Random Forests According to Wager and Athey (2018), the predictor T deﬁned by Deﬁnition 3 is called“k- PNN predictor” if the assumption (A-3) is satisﬁed. Then, we deﬁne our random forests following Athey
et. al (2019).
Deﬁnition 4. For a ﬁxed x ∈ Rd , q ∈ R, τ ∈ (0, 1), our forest score is deﬁned byΨτT (q, x) :=TX
1 X
T (q, x; A, ξ, DT ) =:αt (x)ψqτ (Yt )
|As |
t=1
A∈As

where αt (x) = |A1s |P

Quantile estimatorA∈As αA,t (x) and αA,t (x) =1{Xt ∈LA (x)}.
|LA (x)|By using the above random forest, we can deﬁne an estimator of the conditional quantile

as follows.
Deﬁnition 5. For each τ ∈ (0, 1) and given Xt = x, we deﬁne an estimator of q∗τ (x) by
( T
)
X
τ
τ
q̂T (x) ∈ arg min
αt (x)ψq (Yt )
.
q∈R

t=1

2

Theorem 1. Under some regularity conditions and subsample size s(T ) satisﬁes s(T )/T → 0 and s(T ) → ∞as T → ∞. For each τ ∈ (0, 1) and x ∈ Rp , any sequence of estimators q̂Tτ (x) converges in probability to q∗τ (x),that is,pq̂Tτ (x) → q∗τ (x)as T → ∞.

参考文献
An, H. Z., and F. C. Huang. (1996). The geometrical ergodicity of nonlinear autoregressive models. Statist.
Sinical, 6(4), 943-956.
Athey, P., Tibshirani, J., and Wager, S. (2019). Generalized random forests. Annals of Statistics, 47(2), 11481178.
Davis, Richard A. and Nielsen, Mikkel S. (2020). Modeling of time series using random forests: theoreticaldevelopments. Electronic Journal of Statistics, 14(2), 3644-3671.
Meinshausen, N. (2006). Quantile regression forests. Journal of Machine Learning Research (JMLR), 7, 983-999.
Wager, S. and Athey, P. (2018). Estimation and inference of heterogeneous treatment eﬀects using random
forests. Journal of the American Statistical Association, 113(523), 1228-1242.
2

Analysis of Spatially Resolved Spectral Map of a Galaxy by theHigh-Dimensional Statistical Method
Tsutomu T. TAKEUCHI1,2, Kai T. KONO1, Kazuyoshi YATA, Makoto AOSHIMA3, Aki ISHII4,
Kento EGASHIRA5, Kouichiro NAKANISHI6, Suchetha COORAY1,†, Kotaro KOHNO7
1. Division of Particle and Astrophysical Science, Nagoya University, Nagoya 464-8602, Japan
2. The Research Center for Statistical Machine Learning, the Institute of Statistical Mathematics, Tachikawa, Tokyo 190-8562, Japan
3. Institute of Mathematics, University of Tsukuba, Ibaraki 305-8571, Japan
4. Department of Information Sciences, Faculty of Science and Technology, Tokyo University of Science, Chiba 278-8510, Japan
5. Graduate School of Science and Technology, University of Tsukuba, Ibaraki 305-8571, Japan
6. The ALMA Project, National Astronomical Observatory of Japan, Mitaka, Tokyo 181-8588, Japan
7. Institute of Astronomy, the University of Tokyo, Mitaka, Tokyo 181-0015, Japan† JSPS Research Fellow (DC1)

1. Galaxy Formation and Evolution
A galaxy is a huge agglomeration of stars, interstellar medium (ISM: gas and dust), and dark matter(DM), a complex system with a complicated interaction between each component. There are a few
billions of galaxies in the observable Universe, which delineate the appearance of the visible Universeat optical wavelengths. However, galaxies have formed from a tiny (order of ~ 10-5) fluctuation ofmatter (mainly DM) in the early Universe, when the age of the Universe was only 380,000 yr. The
initial Gaussian fluctuations of DM start to grow by gravitational interactions, finally to form
virialized structures called dark halos. The dark halos approach each other and finally merge to form
larger halos. The formation proceeds from smaller to larger structures. This is the so-calledhierarchical structure formation, currently the most reliable scenario of the structure formation in the
Universe. During the merging of dark halos, the baryonic gas falls into the gravitational potential
wells of DM and is compressed there to form first stars and galaxies. Finally, some galaxies merge
and form larger galaxies. Present-day large galaxies (up to Mbaryon ~ 1012 M☉) have formed in themerger process. Strong merging process is often accompanied by an effective compression of gas,
inducing burst of star formation. Therefore, galaxy evolution is a highly complex process that depends
on the environment of galaxies (number density of ambient galaxies and gas), as well as their internal
processes. The formation and evolution of galaxies are regarded as one of the most importantphenomena in the history of the Universe ranging of 13.8 billion years.

2. Evolution of the Interstellar Medium and Star Formation in Galaxies
2.1 Star formation in the interstellar medium
The most important internal physical process that drives the galaxy evolution is the star formation.
The star formation proceeds as follows: 1) the interstellar medium (ISM) gravitationally contracts, 2)
the ISM cools down and changes its phase as ionized, atomic, and then molecular gas, 3) high-density

molecular gas clumps are formed, and 4) the nuclear fusion ignites. Since stars form in molecular
clouds, it is important to elucidate the physical state of the molecular clouds for the understanding ofgalaxy evolution. Observationally, an empirical relation between the star formation rate (SFR) andthe gas mass density is known (Kennicutt-Schmidt law: Schmidt 1959; Kennicutt 1989, 1998). The
understanding of the fundamental processes in the ISM is also expected to be a clue to explain therelation.
2.2 Spectroscopic observation of the ISM
Direct laboratory experiments are impossible for most of the problems in astronomy. Instead, onlya unique method to obtain the physical information of the ISM in remote objects. Particularly, since
the electromagnetic emission of astronomical molecules are mainly emitted as the radio emission
lines, the astrophysics of the molecular clouds has been developed in radio astronomy. The largest
and most efficient observational facility to observe molecular emission lines is the Atacama Large
Millimeter/Submillimeter Array (ALMA).
Modern astronomical spectroscopic facilities as ALMA can provide tremendous amount of
information on the atoms, molecules, ions, and dust in the ISM. However, spectroscopic observations
are very time-consuming in general, and mapping an extended object by observing many positions
on the sky is not easy. As a result, the number of independent sample size (number of observed
positions) n is much smaller than the wavelength/frequency dimension d, i.e., n ≪ d. Such data arecalled high-dimensional low-sample size (HDLSS) data. Such a problem was regarded as an ill-posed
problem in traditional astronomy and has been simply given up to study further. It has long been
believed that we must make d smaller than n by throwing away the information of the data, in orderto analyze such type of problem. Obviously a method to make a maximal use of the information of
the data is desirable.
2.3 High-dimensional statistics
In a research field other than astronomy, for example in genomics, it is not rare to handle data withthe dimension of the base sequence d ~ 105 while the sample size n ~ 100. A new statistical method
to deal with the HDLSS data has been developed in recent decade, which is referred to as the highdimensional statistical method. The high-dimensional statistics is continuously providing new
findings (e.g., Aoshima 2018).
In this study, we apply the high-dimensional statistics to the spectroscopic mapping data obtainedby ALMA. As mentioned above, the observational cost of spectroscopic mapping is very expensive,and for the ALMA observation, the sample size is at most n ~ 200, while the frequency dimension
is d ~ 2000. Thus, the spectroscopic map of the ALMA is typically HDLSS. In traditional astronomy,
they extracted some emission lines that are already known to be useful and used them to classify
objects. However, ALMA revealed that molecular emissions appear to be significantly different
even between neighboring molecular clouds in the same galaxy (Fig. 1), and the limitation of thetraditional method was clarified. The excessively huge amount of information obscures the relation
between the evolution of the molecular cloud and star formation. The high-dimensional statistics is
expected to make it possible to classify the molecular clouds by using all the emission lines, which
has been considered to be forlorn.

Figure 1: The spectral map of the central region of a starburst galaxy NGC 253 by ALMA (adopted from
Ando et al. 2017). Left: spectra of star forming regions, Right: corresponding regions in NGC 253.

3. Application of the High-Dimensional PCA
The spatial dimension of the map of the central region of NGC 253 is 231, and the spectral
dimension is 2248, namely n = 231 and d = 2248, typical HDLSS data. We applied the highdimensional statistical analysis to the original ALMA spectral mapping data of NGC253 (Fig. 2).

Figure 2: The result of the high-dimensional PCA to the NGC 253. Left: eigenvalue distribution, Right:
distribution of the PC1 and 2.

It is remarkable that the complicated molecular spectra are characterized by first a few PCs. The
PC1 and 2 are particularly dominant, and show a butterfly-like pattern. The physical meaning of thesefeatures are simple. The molecular clouds are coherently rotating around the center of the galaxy,
which causes a systematic Doppler shift of the spectral lines. Thus, PC1 represents the total intensityof the spectral lines, and PC2 shows the Doppler shift (blueshift and redshift) of the lines.
Then we further proceeded the PCA by eliminating the effect of the Doppler shift to explore more
subtle features in the spectra. The result of the PCA to the Doppler-corrected spectra is presented inFig. 3. New PC2 is much smaller than the original one, and there is no butterfly pattern in the PC
distribution. It shows that the systemic rotation is well eliminated. By examining the spatial
distribution of the PCs on the map of the star forming regions, we discovered that the new PC2 andPC3 represent the effect of local expansion of these regions.

Figure 3: The result of the high-dimensional PCA to the Doppler-corrected spectral map of NGC 253. Left:eigenvalue distribution, Right: distribution of the PC1, 2, and 3.

4. Conclusion
Galaxies ubiquitously exist in the present-day Universe, but they have been formed from a tiny
fluctuation of matter in the early Universe. Evolution of galaxies is mainly driven by the starformation, a transition from ISM to stars. Various phases of the ISM are related, and the evolution of
the ISM is a key to complete the understanding of the galaxy evolution. Spectroscopic observations
are of vital importance to extract and interpret the information of matter in galaxies. Spectroscopic
mapping and similar methods are fundamentally important to reveal the ISM physics, but the data are
high-dimensional low sample size. We summarize the conclusions from this study.
1.

We applied the high-dimensional PCA on the NGC 253 spectral map. ALMA mapping data
are typically HDLSS in general, and in this case n = 231 and d = 2228.

2.

Very large variety in the molecular line spectra of NGC253 map can be described only by twoPCs. Each PC consists of ~ 20 elements, much fewer than d. Because these elements may be
a part of same features, the key features may be reduced to several.

3.

The high-dimensional PCA successfully chose two PCs that reproduce the general propertiesof the ALMA spectroscopic map of NGC253.

4.

The controlling feature was HCN(4-3) rotational lines. PC1 describes the total intensity of thelines, and PC2 represents the Doppler shift caused by the systemic rotation.

5.

After correcting the Doppler shift due to the systemic rotation, we could obtain informationon the smaller-scale velocity field described by PC2 (new) and PC3. These are caused by
outflow of starburst regions.

We stress that this result is not obtained by choosing a handful of features by hand, but by making
use of the full information of the high-dimensional data.

Statistical Inference for Glaucoma Detection 1
Ngai Hang Chan
Department of Statistics
The Chinese University of Hong Kong
Shatin, NT, Hong Kongnhchan@sta.cuhk.edu.hk
Abstract
This talk will review some of the statistical techniques used in theanalysis and detection of one of the most commonly encountered eye
diseases, Glaucoma. By means of some of the recent artificial intelligence algorithms, ophthalmologists are employing modern machinelearning technologies such as CNN and its variates to assist in early
detection of eye symptoms in Glaucoma studies. This talk will addresssome of the related issues and discuss the potential of statistical ideas
in these studies. Some examples will be given.

1

September 3–5 , 2021, Time Series Conference, Niigata University, Niigata, Japan.
Research supported in part by grants from HKSAR-RGC-TBS and HKSAR-RGC-GRF.

Mixtures of Nonlinear Poisson Autoregressions
Paul Doukhan1 Konstantinos Fokianos 2 Joseph Rynkiewicz3
1 AGM, UMR 8088, University of Cergy-Pontoisee-mail: doukhan@u-cergy.fr
2 Department of Mathematics & Statistics, University of Cyprus

e-mail: fokianos@ucy.ac.cy
3 SAMM, Université de Paris 1

e-mail: joseph.rynkiewicz@univ-paris1.fr

1Introduction

We consider models for count time series, which allow for the mean process to change according to thevalues of an unobservable discrete random variable. Such models are directly related to Markov switching
AR models (see Hamilton (1994, Ch.22), Cappé et al. (2005) and Frühwirth-Schnatter (2006, Ch. 11-12)).
This class of processes is defined by regime specific models where transition among different regimes
is determined by the state of an unobserved Markov chain. The conditional distribution of the process,given the past and the specific regime, is assumed to be Poisson with a time-varying mean modelled by anon-linear autoregressive infinite order model. The Poisson distributional assumption has been employed
for modeling count time series by several authors including Rydberg and Shephard (2000), Davis et al.
(2003), Ferland et al. (2006) and Fokianos et al. (2009), among others. The aim of this contribution is toextend this framework to Markov switching Poisson processes which can deal with complex dynamics
and take into account several issues found in count data including proper modeling of non-linearities,
overdispersion, multimodal conditional distribution and outliers. In the context of Gaussian time series,
these issues were examined by Le et al. (1996) and Wong and Li (2000), among others.
Count time series analysis is a research topic receiving considerable attention over the last years, see
Kedem and Fokianos (2002, Sec 4 & 5) and the recent edited volume by Davis et al. (2016) for several
additional references. INARCH (INteger ARCH) and INGARCH (INteger GARCH) processes have beenfound useful in applications because they allow for estimation, model assessment and forecasting by
employing existing statistical software. They belong to a broad class of models that falls under the exponential family framework; see Douc et al. (2017) who provide a unified point of view for time series
generalized linear models. This study includes models and associated inference for mixtures of counttime series models, a topic which has not attracted a lot of attention. Related early work was given byAlbert (1991). Subsequently Carvalho and Tanner (2005, 2007) proposed a mixture-of-experts approach tomodel nonlinearities in count time series models. These authors studied maximum likelihood estimation,
investigated identifiability and asymptotic normality of the estimates, and employed the AIC and BIC for
selecting the number of mixture components in the model. However, they imposes strict conditions fordeveloping asymptotic inference and their study focused explicitly on log-linear models for count timeseries. A mixture of linear models was considered by Zhu et al. (2010) but the authors did not provide theoretical evidence about its properties. More recently, Berentsen et al. (2018) applied a Markov-switching
Poisson log-linear model autoregressive model to a study of corporate defaults. For the case of continuous
valued time series, studies of mixture models were given by Jacobs et al. (1991), Le et al. (1996), Francq andRoussignol (1998) and Francq and Zaköian (2001), among others. In addition, Wong and Li (2000, 2001)1

have used a two-component mixture model to extend AR and conditional ARCH models, respectively.
Stability properties of such models have been studied by Saikkonen (2007) and some further work is given
by Kalliovirta et al. (2015). In the context of GARCH models see Francq et al. (2001) while t-distributed
autoregressive models were considered by Wong et al. (2009). Mixtures of AR models, under a Bayesianframework, have been considered by Wood et al. (2011), among others.
The main goals of this work is to fill several gaps in the literature of mixture models for count timeseries models. Precisely,
1. We introduce Markov switching non linear autoregressive Poisson models and study their properties
by considering two cases: (a) For the mixture setup, we show that infinite order models are weakly
dependent. (b) For the Markov switching case, we prove that finite order models are geometricβ-mixing.
2. We study the statistical properties of the maximum likelihood estimator for the case of finite order
autoregression of Markov switching models.
3. We show that the marginal likelihood ratio test for testing the number of hidden regimes converges
to a Gaussian process. This fact implies that the BIC estimator is consistent estimator for selecting
the true number of regimes.
4. We provide a counterexample which shows that estimation of the true number of regimes, under amisspecified GARCH type count time series model, is not possible.
5. We apply the theoretical results to the weekly number of E.coli cases and compare mixture models
to models with interventions.

References
Albert, P. S. (1991). A two–state Markov mixture model for a time series of epileptic seizure counts. Biometrics 47,
1371–1381.
Berentsen, G. D., J. Bulla, A. Maruotti, and B. Støve (2018). Modelling corporate defaults: A Markov-switching Poissonlog-linear autoregressive model. ArXiv e-prints. https://arxiv.org/abs/1804.09252.
Cappé, O., E. Moulines, and T. Rydén (2005). Inference in Hidden Markov Models. New York: Springer.
Carvalho, A. X. and M. A. Tanner (2005). Modeling nonlinear time series with local mixtures of generalized linear
models. Canad. J. Statist. 33, 97–113.
Carvalho, A. X. and M. A. Tanner (2007). Modelling nonlinear count time series with local mixtures of Poisson
autoregressions. Comput. Statist. Data Anal. 51, 5266–5294.
Davis, R. A., W. T. M. Dunsmuir, and S. B. Streett (2003). Observation-driven models for Poisson counts. Biometrika 90,
777–790.
Davis, R. A., S. H. Holan, R. Lund, and N. Ravishanker (Eds.) (2016). Handbook of Discrete-Valued Time Series. Handbooks of Modern Statistical Methods. London: Chapman & Hall/CRC.
Douc, R., K. Fokianos, and E. Moulines (2017). Asymptotic properties of quasi-maximum likelihood estimators in
observation-driven time series models. Electron. J. Stat. 11, 2707–2740.
Ferland, R., A. Latour, and D. Oraichi (2006). Integer–valued GARCH processes. Journal of Time Series Analysis
Analysis 27, 923–942.
Fokianos, K., A. Rahbek, and D. Tjøstheim (2009). Poisson autoregression. Journal of the American Statistical Association 104, 1430–1439.
Francq, C. and M. Roussignol (1998). Ergodicity of autoregressive processes with Markov-switching and consistencyof the maximum-likelihood estimator. Statistics 32, 151–173.

2

Francq, C., M. Roussignol, and J.-M. Zakoïan (2001). Conditional heteroskedasticity driven by hidden Markov chains.
Journal of Time Series Analysis 22, 197–220.
Francq, C. and J.-M. Zaköian (2001). Stationarity of multivariate Markov-switching ARMA models. J. Econometrics 102,
339–364.
Frühwirth-Schnatter, S. (2006). Finite mixture and Markov switching models. Springer Series in Statistics. Springer, NewYork.
Hamilton, J. D. (1994). Time series analysis. Princeton University Press, Princeton, NJ.
Jacobs, R., M. Jordan, S. Nowlan, , and G. Hinton (1991). Adaptive mixtures of local experts. Neural Computation 3,
79–87.
Kalliovirta, L., M. Meitz, and P. Saikkonen (2015). A Gaussian mixture autoregressive model for univariate time series.
J. Time Series Anal. 36, 247–266.
Kedem, B. and K. Fokianos (2002). Regression Models for Time Series Analysis. Hoboken, NJ: Wiley.
Le, N. D., R. D. Martin, and A. E. Raftery (1996). Modelling flat stretches, bursts, and outliers in time series using
mixture transition distribution models. Journal of the American Statistical Association 91, 1504–1515.
Rydberg, T. H. and N. Shephard (2000). A modeling framework for the prices and times of trades on the New Yorkstock exchange. In W. J. Fitzgerlad, R. L. Smith, A. T. Walden, and P. C. Young (Eds.), Nonlinear and Nonstationary
Signal Processing, pp. 217–246. Cambridge: Isaac Newton Institute and Cambridge University Press.
Saikkonen, P. (2007). Stability of mixtures of vector autoregressions with autoregressive conditional heteroskedasticity.
Statist. Sinica 17, 221–239.
Wong, C. S., W. S. Chan, and P. L. Kam (2009). A Student t-mixture autoregressive model with applications toheavy-tailed financial data. Biometrika 96, 751–760.
Wong, C. S. and W. K. Li (2000). On a mixture autoregressive model. Journal of the Royal Statistical Society B62, 95–115.
Wong, C. S. and W. K. Li (2001). On a mixture autoregressive conditional heteroscedastic model. Journal of the American
Statistical Association 96, 982–995.
Wood, S., O. Rosen, and R. Kohn (2011). Bayesian mixtures of autoregressive models. J. Comput. Graph. Statist. 20,
174–195.
Zhu, F., Q. Li, and D. Wang (2010). A mixture integer-valued ARCH model. J. Statist. Plann. Inference 140, 2025–2036.

3

ブリッジ推定量におけるチューニングパラメーターの選
択について
高崎経済大学・経済学部

宮田 庸一

2021 年 8 月 15 日
′ 記号は, 行列の転置を表すものとする. 観測ベクトル (Yi , x′i )′ は, 以下のスパースな線形モデルに従うとする.

Yi = β0∗ + β1∗ x1i + · · · + βk∗n xkn ,i + ui ,(i = 1, ..., n)= β0∗ + β1∗ x1i + · · · + βk∗n xkn ,i + 0 · xkn +1,i + · · · + 0 · xpn i + ui ,ただし, xi = (x1i , · · · , xpn ,i )′ は pn 次元の説明変数を表すベクトルとし, Yi は被説明変数とする. また β ∗ = (β0∗ , β1∗ , .., βk∗n , 0, · · · , 0)′ ∈ Rpn +1 は真のパラメーターのベクトルとし, βj∗ ̸= 0

(j = 1, ..., kn ), および βj∗ = 0 (j = kn + 1, .., pn ) であるとする. ui は期待値 E(ui ) = 0, 分散
V ar(ui ) = σ02 となる撹乱項とする. ここで, 切片以外の真のパラメーターベクトルにおける非ゼロの要素の添え字の集合を S0,n = {j ∈ {1, ..., pn }|β0,j ̸= 0} とし, その要素の個数を kn := |S0,n | で表すとことにする. ここでは kn と説明変数の個数 pn は, 標本の大きさ n に依存してもよいことにする. このモデルにより生成された標本 (Yi , x′i )′ (i = 1, ..., n) に対して, 以下線形モデルを当てはめることを考える.

Yi = β0 + βj1 xj1 i + · · · + βjk xjk i + ϵ i ,(i = 1, ...n),ただし, j1 , ..., jk ∈ {1, ..., pn }, k = 1, ..., pn とする. このような説明変数の選択問題において,pn の次元がそれほど大きくない場合には, AIC, BIC などの情報量規準を用いて総当たり法を行うのが標準的なアプローチであろう. しかしその一方で, pn の次元が高くなるにと, 評価すべきモデルの数が指数的に増加するため, 総当たり法を行うのが困難になる. このような問題に対する一つの標準的なアプローチは, Tibshirani (1996) による LASSO(Least absolute shrinkage

and selection operator) に代表される罰則付き最小二乗推定量を用いることである. 具体的には,()
y = (Y1 , ..., Yn )′ , X1 = (xji ):n × pn 行列, 1n = (1, ..., 1)′ , X = 1n

X1 とおき, 損失関数

L0,n (β) = ∥y − Xβ∥2 + λn ∥β∥1
を最小にする β の値 β̂L を求めればよい. ただし, ∥β∥1 =

(1)

∑pn

j=1 |βj | とする.

この推定量はLASSO 推定量と呼ばれ, λn → ∞ (n → ∞) とその他の条件のもとで漸近正規性を持つことが知られている. 一方で, LASSO 推定量をモデル選択手法の観点から見たときには, いくつかの問題が生じる. それは収束レートが√n である漸近正規性が成り立つためのチューニングパラメーター

λn のオーダーの仮定の下では, 真のモデルを選ぶ確率は漸近的に 1 にならないことが知られている. この問題を改善するために, 式 (1) における罰則項を, 非凸の形のものに置き換えた以下の損失関数を考え, これを最小にする推定量 β̂ を求める:

Ln (β) = ∥y − Xβ∥2 + λnpn∑

|βj |γ ,(2)

j=1

ただし 0 < γ < 1 とする. この推定量はブリッジ推定量と呼ばれている. 当然のことながら, ブ
リッジ推定においても適切な λn の選び方が重要になる. Umezu et al. (2019) においては, pn はn に依存しない条件, しかし線形モデルより広い一般化線形モデルの下での AIC 型の情報量規準を与えている. 一方で, Wang et al. (2009) では, 線形モデルの下で, モデル選択に関して一致性を持
つ BIC 型の情報量規準を与えている. Huang et al. (2008) では, オラクル性と呼ばれる, 推定量

β̂ が漸近正規性, 漸近有効性を持ち, なおかつ説明変数に関しても一致性を持つような条件を与えている. 例えば γ = 1/2 として, kn = k(定数), pn = log n とおいた場合, λn のオーダーに関する条件は, λn /(n1/4 (log n)3/4 ) → ∞ (n → ∞), λn = o(n1/2 ) となるが, そのような λn の選択肢は,

λn = 3n3/8 でも良いし, λn = 10n3/8 でもよいことになる. すなわち Huang et al. (2008) の結果
は λn の選択に関しては, 何も述べていないことになる. このため, 今回の講演では, Huang et al.

(2008) の条件の下で λn を選ぶための規準を, 一般化周辺尤度に対する近似として与えられることを示し, シミュレーションによりそのパフォーマンスを検証した. ここで用いた一般化周辺尤度とは, Yamanishi (1998) により与えられた, 重みつき対数尤度に (−1) を乗じたものを損失関数とした拡張型確率的コンプレキシティに対応するものである.

参考文献
Huang, J., J. L. Horowitz, and S. Ma (2008). Asymptotic properties of bridge estimators insparse high-dimensional regression models. Ann. Statist. 36 (2), 587–613.
Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. J. Roy. Statist. Soc.
Ser. B 58 (1), 267–288.
Umezu, Y., Y. Shimizu, H. Masuda, and Y. Ninomiya (2019). AIC for the non-concave
penalized likelihood method. Ann. Inst. Statist. Math. 71 (2), 247–274.
Wang, H., B. Li, and C. Leng (2009). Shrinkage tuning parameter selection with a diverging
number of parameters. J. R. Stat. Soc. Ser. B Stat. Methodol. 71 (3), 671–683.
Yamanishi, K. (1998). A decision-theoretic extension of stochastic complexity and its applications to learning. IEEE Trans. Inform. Theory 44 (4), 1424–1439.

３次元分割表における条件付き独立性検定統計量の変換統計量について
北海道教育大・札幌 種市信裕
北海道教育大・釧路 関谷祐里
数学利用研究所 外山淳

1. ３次元分割表における条件付独立性モデル．
３次元の J × K × L 分割表において多項分布モデルを考える. (j, k, l) セルの観測度数∑L∑∑を表す確率変数を Xjkl とする. ただし，Xjkl は Jj=1 K
k=1
l=1 Xjkl = n を満たす非負整数の値を取り，自然数 n は定数とする. また，(j, k, l) セルのセル確率を pjkl とする. このとき，確率変数ベクトル X = (X111 , . . . , XJKL )′ が多項分布 MJKL (n, p)に従う場合について考える. ここで，p = (p111 , . . . , pJKL )′ である. このとき，このJ × K × L 分割表のセル確率に関する条件付き独立性の仮説

H0 : pjkl =pj·l p·kl(j = 1, . . . , J, k = 1, . . . , K, l = 1, . . . , L)p··lを検定するための ϕ-divergence に基づく検定統計量は
(

J ∑
K ∑
L
∑

p̃jkl
Cϕ = 2np̂jkl ϕp̂jkl
j=1 k=1 l=1)である. ここで，p̃jkl = Xjkl /n, p̂jkl = Xj·l X·kl /(nX··l ), であり，ϕ は (0, ∞) において定義された実凸関数であり，ϕ(1) = ϕ′ (1) = 0 と ϕ′′ (1) = 1 を満たすものとする.

2. 独立性検定統計量の分布の漸近展開に基づく近似．
H0 のもとでの検定統計量 Cϕ の分布の漸近展開に基づく次のような近似を考える.
Pr{Cϕ ≤ x|H0 } ≈ W1 + W2 , ここで，W1 は多変量エッジワース展開に基づく項，W2は不連続性を考慮に入れた離散項である. W2 は，多項分布の適合度検定におけるピ
アソンのカイ二乗統計量の分布に対して，Yarnold [4] により与えられた漸近展開に基づく近似式で導入された項に対応する. W1 に関して以下の定理が成り立つ.
定理. (Taneichi et al. [3])関数 ϕ(t) が６回微分可能で，ϕ(6) (t) が t = 1 で連続ならば，W1 は次のように評価される.

1
W1 = Pr{χ2M ≤ x} +3∑ϕv Pr{χ2M +2j ≤ x} + O(n−2 ),n j=0 j

ϕ

ϕ

(1)

′′′

ここで，M = (J − 1)(K − 1)L, v0 = γ0 , v1 = (γ1 + γ2 ) + 2ϕ (1)γ2 + ϕ(4) (1)γ3
′′′
′′′
′′′
+{ϕ (1)}2 γ4 , v2ϕ = 2(−γ2 + γ3 ) − 2ϕ (1)(γ2 + γ4 ) − ϕ(4) (1)γ3 − 2{ϕ (1)}2 γ4 , v3ϕ
′′′
1
1
1
Γ1 − 12
Γ2 + 12(Γ3 + Γ4 ), γ(1 = 14 JKΓ1 + 41 Γ2 − 14 (KΓ= γ4 {1 + ϕ (1)}2 , γ0 = − 12
) 3

+JΓ4 ), γ2 = 18 J 2 K 2 Γ1 + 18 Γ2 − 81 (K 2 Γ3 + J 2 Γ4 ), γ3 = − 12 JK − 81 + 14 J + 41 K Γ1

− 18 Γ2 + 14 (KΓ3 + JΓ4 ) − 81 (Γ3 + Γ4 ), γ4 =

(

1 2 2
3
1
1
1
8 J K + 4 JK + 3 − 2 J − 2 K

)

Γ1

5
+ 24
Γ2 − 18 (K 2 Γ3 + J 2 Γ4 ) − 14 (KΓ3 + JΓ4 ) + 16 (Γ3 + Γ4 ),

Γ1 =

L
∑

p−1
··ℓ , Γ2 =

J ∑
K ∑
L
∑

−1
qjkl
, Γ3 =

j=1 k=1 ℓ=1

ℓ=1

J ∑
L
∑
j=1 ℓ=1

p−1
j·ℓ , Γ4 =

K ∑
L
∑

p−1
·kℓ , qjkl =k=1 ℓ=1

pj·l p·kl
p··l

である.

3. 漸近展開の連続項に基づく変換統計量．
本報告では, 少ない標本数でも Cϕ より速く極限分布に近づく検定統計量の構築を目指す. 導出された離散項 W2 の評価式は非常に複雑であることもあり Cϕ の分布に対して，(1) 式で与えられる W1 の近似式のみを用いて，小標本におけるカイ二乗近似の改良を考える. (1) 式にバートレット修正および改良変換の構築と漸近展開式との関係の理論 (e.g., Fujikoshi [2], Cordeiro and Ferrari [1]) を適用した. (1) 式において,
ϕ′′′ (1) = −1 かつ ϕ(4) (1) = 2 が成り立つ場合には v1ϕ = −v0ϕ かつ v2ϕ = v3ϕ = 0 が成り立つ. このことは, バートレット修正が可能なことを意味する. よってこの場合には,
Cϕ にバートレット修正を施した変換統計量 (Cϕ )B = {1 + 2v0ϕ (nM )−1 }Cϕ を構築できる. (1) 式においてその他の場合には, 改良変換統計量 ([2])[(Cϕ )I=(nα + β)2 log{

1 ((C )2 + γ(C )3 )1 + 1 2 Cϕ + nαϕϕ(nα) (
)}]
3γ
9γ 2
1
3
4
5
1
+
(Cϕ ) +
(Cϕ ) +
(Cϕ )
,
4
20
(nα)2 3

ただし, α = −M (M + 2){2(v2 + v3 )}−1 , β = −(M + 2)v0 {2(v2 + v3 )}−1 , γ =v3ϕ {(M + 4)(v2ϕ + v3ϕ )}−1 , および, バートレット型変換統計量 ([1])ϕ{(Cϕ )CF =

ϕ

ϕ

ϕ

ϕ

}

2v ϕ
2(v0ϕ + v1ϕ )2(v0ϕ + v1ϕ + v2ϕ )1+ 0 +Cϕ +(Cϕ )2 CϕnMnM (M + 2)
nM (M + 2)(M + 4)を構築できる.

参考文献[1] Cordeiro, G. M. and Ferrari, S. L. P.: Biometrika, 78, (1991), 573–582.
[2] Fujikoshi, Y.: Journal of Multivariate Analysis, 72, (2000), 249–263.
[3] Taneichi, N., Sekiya, Y. and Toyama, J.: Journal of Multivariate Analysis, 171,(2019), 193-208.
[4] Yarnold, J. K.: Ann. Math. Statist., 43, (1972), 1566–1580.