令和元年度科学研究費基盤研究 (A)（課題番号：15H01678）「大規模複雑データの理論と方法論の総合的研究」（研究代表者：青嶋誠（筑波大学））によるシンポジウム

統計的推測および確率解析に関する総合的研究
開催責任者：宇野力 (秋田大学教育文化学部)
場所：カレッジプラザ

講堂

住所：秋田市中通 2 丁目 1-51

明徳館ビル２階 (JR 秋田駅近く)

日時：令和元年 12 月 11 日（水）14:00 〜 12 月 13 日（金）12:20

プログラム
１２月１１日（水）14:00〜14:10 開会
座長 種市信裕 (北海道教育大学札幌校)

14:10〜14:50 「Non-stationary Time Series Classification using Topological Data Analysis」
黒木裕鷹 (東京理科大学大学院工学研究科), 塩濱敬之 (東京理科大学工学部)

14:50〜15:30 「Maximum Likelihood Estimation for Hidden Markov Models with Cylindrical
Distributions」
小坪琢人 (東京理科大学大学院工学研究科), 塩濱敬之 (東京理科大学工学部)
座長 塩濱敬之 (東京理科大学工学部)

15:45〜16:25 「Jensen’s alpha in the CAPM market model from the viewpoint of asymptotic
accuracy of estimators of β and of asymmetry of distributions of residuals」
三浦良造 (東北大学大学院経済学研究科, 一橋大学)

16:25〜17:05 「正方分割表における部分対称性に関する種々のモデルとその尺度」
桃崎智隆 (東京理科大学理工学研究科), 福本伸樹 (東京理科大学理工学研究科),
中川智之 (東京理科大学理工学部), 石井晶 (東京理科大学理工学部),
三枝祐輔 (横浜市立大学医学部), 富澤貞男 (東京理科大学理工学部)

17:05〜17:45 「多次元分割表における多要因間独立性と多次元完全独立性の関係」
種市信裕 (北海道教育大学札幌校), 関谷祐里 (北海道教育大学釧路校),
外山淳 (数学利用研究所)

１２月１２日（木）座長 柿沢佳秀 (北海道大学経済学研究院)

10:15〜10:55 「単一強スパイク固有値モデルに対する高次元平均ベクトルの 2 標本検定」
石井晶 (東京理科大学理工学部), 矢田和善 (筑波大学数理物質系数学域),
青嶋誠 (筑波大学数理物質系数学域)

10:55〜11:35 「13 の癌腫の 78 個の Microarray による LSD の確認と Method1 による SM と
BGS 分割の評価 –Big データ時代における LSD の重要な 2 つの特徴–」
新村秀一 (成蹊大学)

11:35〜12:15 「ベイズモデルの母数：事前分布と estimand」
柳本武美 (統計数理研究所)

12:15〜13:30 昼休み
座長 矢田和善 (筑波大学数理物質系数学域)

13:30〜14:10 「Gaussian Semiparametric Estimation for Two Dimensional Intrinsic
Stationary Random Fields」
矢島美寛 (東北大学大学院経済学研究科)

14:10〜14:50 「再生核ヒルベルト空間における射影のモーメントによる二標本検定」
牧草夏実 (千葉大学大学院融合理工学府)

14:50〜15:30 「非負データに対する再帰的な密度推定について」
柿沢佳秀 (北海道大学経済学研究院)

座長 島谷健一郎 (統計数理研究所)

15:45〜16:25 「ノンパラメトリックな極値統計解析法の漸近理論」
吉田拓真 (鹿児島大学大学院理工学研究科)

16:25〜17:05 「ノンパラメトリック重回帰分析における Lasso 型推定の改良」
松島佑樹 (島根大学大学院自然科学研究科), 内藤貫太 (千葉大学大学院理学研究院)

17:05〜17:45 「Semiparametric regression with local divergence」
内藤貫太 (千葉大学大学院理学研究院)

１２月１３日（金）座長 内藤貫太 (千葉大学大学院理学研究院)

10:15〜10:55 「市場創造時点の特定に向けた文献調査：KIP を前提として」加藤淳一 (久留米大学商学部)

10:55〜11:35 「Estimation, Diagnostics, and Extensions of Nonparametric Hawkes Processes
with Kernel Functions」
庄建倉 (統計数理研究所)

11:35〜12:15 「非定常クラスター過程のパラメータ推定法：Palm 尤度法の拡張」
島谷健一郎 (統計数理研究所)

12:15〜12:20 閉会

Non-stationary Time Series Classification using
Topological Data Analysis
Yutaka Kuroki
Graduate School of Engineering, Tokyo University of Science
Takayuki Shiohama
Department of Information and Computer Technology

Abstract
Time series classification (TSC) is an important and challenging problem in data mining. There
are hundreds of algorithms for TSC available with the increase of time series data availability. For
more details, we refer recent work of Bagnall et al. (2017) and review paper of Fawaz et al. (2019).
Some of the machine learning algorithms are based on the bug of patterns, and learning patterns ofsimilarity are key feature extraction for time series data. In this study, we employ feature extractionof time series using Topological Data Analysis (TDA). TDA refers to a collection of methods forfinding topological structure in data. Until recently, topological inference relied on deterministic approaches, and it is well known that these inference is easily affected by outliers and/or noisy datasets.
Moreover, the high computational costs are required for computing persistence homology with complex datasets in time and space. To overcome these circumstances, two approaches are proposedin literature. The first approach is the bootstrap estimation for persistence diagrams and landscapesof Chazal et al. (2014). The second refers the methods of subsampling of Chazal et al. (2015). In
this study, we introduce an unsupervised classification learning for several non-stationary time series
using topological features. The bootstrap and subsampling methods are implemented to compare theperformance of classifications.
Keywords: Bootstrap estimation, time series analysis, topological data analysis, subsampling method.

References
Bagnall, A., J. Lines, A. Bostrom, J. Large, and E. Keogh (2017). The great time series classification
bake off: a review and experimental evaluation of recent algorithmic advances. Data Mining andKnowledge Discovery 31(3), 606–660.
Chazal, F., B. Fasy, F. Lecci, B. Michel, A. Rinaldo, and L. Wasserman (2015). Subsampling methods
for persistent homology. In International Conference on Machine Learning, pp. 2143–2151.
Chazal, F., B. T. Fasy, F. Lecci, A. Rinaldo, A. Singh, and L. Wasserman (2014). On the bootstrap for
persistence diagrams and landscapes. arXiv preprint arXiv:1311.0376.
Fawaz, H. I., G. Forestier, J. Weber, L. Idoumghar, and P.-A. Muller (2019). Deep learning for timeseries classification: a review. Data Mining and Knowledge Discovery 33(4), 917–963.

1

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

1

Introduction

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

Hidden Markov models assume a sequence of random variable to be conditionally independent givena sequence of state variables which forms a Markov chain. To provide the consistency and asymptotic
normality of maximum likelihood estimators (MLE), we need to consider identifiability of the model
parameters. Recently, Miyata et al. (2019a) provide a strong consistency results for the mixtures ofasymmetric circular distribution. Miyata et al. (2019b) investigate the identifiability of circular or cylindrical distribution. In this study, we consider the asymptotic properties for the maximum likelihood
estimation for the HMM with cylindrical distribution.
The rest of paper is organized as follows. Section 2 introduces our HMM models. Several useful distributional properties of cylindrical distributions are explained. Section 3 provides the maximum likelihood
estimation procedures which utilizing EM algorithms. Some technical improvements are also discussed
to avoid the local maximum solution in optimization together with to avoid unboundedness in likelihood
functions. In Section 4, some Monte Carlo simulations are performed to assess the performance of our
proposed procedures. Section 5 provides real data analysis of animal movement trajectories. Section 6
concludes our paper.

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

秋田研究集会発表報告書
発表者
三浦良造 一橋大学名誉教授 東北大学経済学研究科特任（客員）教授
発表題目

Jensen’s alpha in the CAPM market model from the viewpoint of asymptotic accuracyof estimators of β and of asymmetry of distributions of residuals.
(This is a joint work with Mr. Navruzbek Karamatov, Ph.D student at Tohoku University)：１．研究の動機
単純線形回帰モデルの傾きβの推定精度とモデルの偶然誤差項の分布の対称性・非対称性を軸にして、回帰モデルの切片の推定値にどの程度の差異が現れる
か、二通りの推定値のどちらを採用すべきかの問題を扱う。
金融のアカデミアと実務においては傾きの推定に最小二乗推定(LSE)を使っている。これは BLUE（最小分散不偏線形推定）と呼ばれる。しかし、順位統計量に基づ
く推定量は非線形推定量であり、外れ値に対してロバストであることが知られている。偶然誤差項の確率分布が裾が重い場合は LSE よりも順位推定を使う方が推定精度がよいといいう意味で適切だろう。
この立場に立って、二つの異なる回帰残差から求められる切片αの推定値にどの程度の差があるかを調べ、それが残差分布の対称性・非対称性とどのような関係を持つかを調べた。
：２．発表概要
偶然誤差項の確率分布を一般化された Lehmann 対立仮説モデル（セミパラメトリックなモデル）であらわし、残差に基づいて、このモデルが含む、位置母数と非対称度を示すパラメターの推定量を構成し、これらの推定量の漸近正規性の数学的証明を報告した。
さらに、セミパラメトリックモデルに含まれる関数形が未知である対称な確率分布 F を推定する推定量を構成し、その漸近正規性の証明を報告した。
これらの数理的結果を支えとして、実際の株式市場でこれらのパラメータがとる数値を算出した。米国、英国、日本の 3 市場の日次データを用いた。
βの推定を二つの方法で行うので、残差が二通り得られ、それぞれがパラメターα、θ、μの推定に使われる。ここから得られる膨大な数値結果をまとめる軸として、（１）αを偶然誤差項分布の位置と非対称効果（非対称性がもたらす期待値の大きさへの貢献量）の和として表すαの分解、さらに（２）算出されたαの二通りのαの値の差が残差の非対称度、偶然誤差項分布の裾の重さとどのように関係するかをグラフに表した。
αは株式ポートフォリオのパフォーマンス分析に用いられるが、LSE による推定値が順位統計量による推定値に比べて過大推定、あるいは過小推定である様子を示す以下のグラフに沿って報告した。そこでは、過大・過少の基準はβの推定精度が高いほうを正しい推定とみなして、正しいほうがもたらす残差から得たαを正しい推定値とした。

: References.
Jaeckel, L.A. (1972). Estimation of regression coefficients by minimizing the dispersionof the residuals. The Annals of Mathematical Statistics. Pp.1449-1458.
Jureckova, J.(1971). Nonparametric estimate of regression coefficients. The Annals of
Mathematical Statistics. Pp.1328-1338.
Miura, R. and Tsukahara, H.(1993). One-sample estimation for generalized Lehmann’salternative models. Statistica Sinica. Pp.83-101.
Miura, R.(2014). A simple linear regression model with Generalized Lehmann’s
Alternative Model and R-estimators. Shogaku Ronnkyu. Vol.61. No.3. pp.89-108 (inJapanese) (Journal of Business Administration Kwansei Gakuin University)正方分割表における部分対称性に関する種々のモデルとその尺度
東京理科大学大学院理工学研究科
東京理科大学大学院理工学研究科
東京理科大学理工学部
東京理科大学理工学部
横浜市立大学医学部
東京理科大学理工学部

桃﨑 智隆
福本 伸樹
中川 智之
石井 晶
三枝 祐輔
富澤 貞男

正方 r × r 分割表を考え, (i, j) セル確率を pij (i = 1, · · · , r; j = 1, · · · , r) とする. 本稿では,{pij + pji > 0} を仮定する.

第 1 部 正方分割表における条件付き部分対称性とその尺度
行変数を X, 列変数を Y とする. McCullagh (1978) によって提案された条件付き対称 (CS) モ
デルは,

P (X = i, Y = j | X < Y ) = P (X = j, Y = i | X > Y )(すべての i < j)により定義される.
一般に, CS モデルが成り立たないとき, CS モデルとの隔たりの程度を測ることに関心がある.
Tomizawa and Saitoh (1999) は CS モデルからの隔たりを測る尺度 Φ(λ) を提案した:

1∑ ∑ U(λ)(pij + pL=ji )Ψij
2
r−1

Φ

(λ)

r

(λ > −1),i=1 j=i+1

ただし,
(λ)

Ψij = 1 −

qij =

λ2λ
(λ)
Gij ,
λ
2 −1

(λ)

Gij =

]
1[
1 − (qij )λ+1 − (qji )λ+1 ,
λ

pU
ij

pL
pij L
pji
ji
,
q
=
, pU
, pji =(i < j),ji
ij =LU + pL
δ
δL
pU
+
p
p
U
ij
ji
ij
ji
δU =

∑∑

pst , δL =

∑∑

s<t

pst ,

s>t

CS モデルが成り立たないとき, CS モデルより制約を緩めたモデルとその尺度に関心がある. そこで, CS モデルより制約を緩めたモデルである条件付き部分対称 (CPS) モデルと, CPS モデルからの隔たりを測る尺度 τ (λ) を提案した.
CPS モデルを以下のように定義する:
P (X = i, Y = j | X < Y ) = P (X = j, Y = i | X > Y ) (少なくとも１組の (i, j), i < j).
また, CPS モデルからの隔たりを測る尺度を次のように提案した:

τ (λ) =

r−1
∏

r
[
∏

i=1 j=i+1

τ (λ) は次の性質を持つ.

] U L(λ) (pij +pji )/2

Ψij(λ > −1)

(i) 0 ≤ τ (λ) ≤ 1.
(ii) τ (λ) = 0 ⇔ 少なくとも 1 組の (k, l), k < l に対して,P (X = k, Y = l|X < Y ) = P (X = l, Y = k|X > Y ) が成立.
(iii) τ (λ) = 1 ⇔ すべての i ̸= j に対して, pij = 0(このとき, pji ̸= 0) が成立.
∑
推定尺度 τ̂ (λ) は近似的に N (τ (λ) , n−1 σ 2 [τ (λ) ]) に従う (n = i,j nij , nij は (i, j) セル度数). 講演
では, τ (λ) に関するシュミレーション結果や信頼区間を用いた実データ解析結果について示した.

第 2 部 正方分割表における部分対称と部分非対称からの隔たりを測るベクトル型尺度
本講演では, 部分的な対称セルで非対称構造を持つモデルとして部分完全非対称 (PAS) モデル
を以下で定義した:

pij = 0 (pji ̸= 0) (少なくとも 1 組の i ̸= j)
また, pcij = pij /(pij + pji ), p∗ij = pij /δ, δ =尺度を次のように提案した:

κ(λ) =

r−1
∏

r
[
∏

∑∑

i̸=j pij とする. PAS モデルからの隔たりを測る] ∗ ∗(λ) (pij +pji )1 − ϕij(λ > −1)i=1 j=i+1

ただし,
(λ)

ϕij = 1 −

]
1[
(λ)
(λ)
c λ+1
c λ+1
H
,
H
=
.
1
−
(p
)
−
(p
)
ij
ji
ij
λ
2λ − 1 ij
λ2λ部分対称 (PS) モデルからの隔たりを測る尺度を Υ(λ) とする. PS モデルと PAS モデルからの隔たりを同時に測るベクトル型尺度を次のように提案した:()⊤
Λ(λ) = Υ(λ) , κ(λ)(λ > −1)
Λ(λ) は次の性質を持つ.
(i) Λ(λ) = (0, 0)⊤ ⇔ PS モデルと PAS モデルが同時に成立.
(ii) Λ(λ) = (0, 1)⊤ ⇔ 対称モデル (すべての i ̸= j に対して pij = pji ) が成立.
(iii) Λ(λ) = (1, 0)⊤ ⇔ 完全非対称モデル (すべての i ̸= j に対して pcij = 1 または pcji = 1) が成立.
b (λ) は近似的に N (Λ(λ) , n−1 Σ) に従う. 講演では, Λ(λ) に関するシミュレーション結果
推定尺度 Λ
や, 信頼楕円を用いた実データ解析について示した.

参考文献
McCullagh, P. (1978). A class of parametric models for the analysis of square contingency tables
with ordered categories. Biometrika, 65:413–418.
Tomizawa, S. and Saitoh, K. (1999). Measure of departure from conditional symmetry for squarecontingency tables with ordered categories. Journal of the Japan Statistical Society, 29:65–78.

多次元分割表における多要因間独立性と多次元完全独立性の関係
種市
関谷
外山

信裕
祐里
淳

(北海道教育大・札幌)(北海道教育大・釧路)(数学利用研究所)

1. ２次元分割表の独立性モデル
２次元の J ×K 分割表において多項分布モデルを考える。(j, k) セルの観測度数を表す∑∑
確率変数を Xjk とする。ただし，Xjk (j = 1, . . . , J, k = 1, . . . , K) は， Jj=1 K
k=1 Xjk =n を満たす非負整数の値を取り，自然数 n は定数とする。また，(j, k) セルのセル確率を pjk とする。ただし，pjk (j = 1, . . . , J, k = 1, . . . , K) は，0 < pjk < 1 (j =∑∑1, . . . , J, k = 1, . . . , K), Jj=1 K
k=1 pjk = 1 を満たすものとする。このとき，p = (p11 ,⊤. . . , pJK ) とおき，確率変数ベクトル X = (X11 , . . . , XJK )⊤ が多項分布 MJK (n, p) に従う場合について考える。
周辺確率及び周辺度数を，

pj· =

K
∑

pjk , p·k =

J
∑
j=1

k=1

pjk , Xj· =K∑

Xjk , X·k =

k=1

J
∑

Xjk (j = 1, . . . , J, k = 1, . . . , K)j=1と定義するとき，J × K 分割表における独立性の仮説
(1)

(1)

H0 : pjk = pj· p·k(j = 1, . . . , J, k = 1, . . . , K)を検定するための ϕ-divergence に基づく統計量は，
J ∑K∑
Xj· X·k(1)
Cϕ = 2nϕn2j=1 k=1(

nXjk
Xj· X·k)で定義される。ここで，ϕ(t) は，区間 (0, ∞) における凸関数で，ϕ(1) = ϕ′ (1) = 0 と(1)ϕ′′ (1) = 1 を満たすものとする。帰無仮説 H0 のもとで，n → ∞ としたときの検定(1)
統計量 Cϕ の極限分布は，すべての関数 ϕ に対して，自由度 (J − 1)(K − 1) のカイ二乗分布である。

2. 多次元分割表の独立性モデル
３次元の J × K × L 分割表において多項分布モデルを考える。(j, k, ℓ) セルの観測度数を表す確率変数を Xjkℓ ，セル確率を pjkℓ とする。このとき，確率変数ベクトル
X = (X111 , . . . , XJKL )⊤ が多項分布 MJKL (n, p) に従う場合について考える。ここで，p = (p111 , . . . , pJKL )⊤ である。
周辺確率と周辺度数に対しては，前節で用いた記法と同様に，和を取った添え字を· で表すことによって，記号 pj·· , p·k· , p··ℓ , pjk· , pj·ℓ , p·kℓ 等を定義することとする。このとき，J × K × L 分割表のセル確率に関する以下の仮説を考える。

(2)
(3)

(2)

H0 : pjkℓ = pj·· p·k· p··ℓ (j = 1, . . . , J, k = 1, . . . , K, ℓ = 1, . . . , L)
(3)
H0 : pjkℓ = pj·· p·kℓ (j = 1, . . . , J, k = 1, . . . , K, ℓ = 1, . . . , L)

(2)

(3)

(2) は完全独立性の仮説，(3) は１要因対２要因の独立性の仮説である。帰無仮説 H0 , H0を検定するための ϕ-divergence に基づく統計量はそれぞれ，
J ∑
K ∑
L
∑
Xj·· X·k· X··ℓ
(2)
Cϕ = 2nϕn3
j=1 k=1 ℓ=1
J ∑
K ∑
L
∑
Xj·· X·kℓ(3)
Cϕ = 2nϕn2j=1 k=1 ℓ=1
(2)

(

(

n2 Xjkℓ
Xj·· X·k· X··ℓnXjkℓ
Xj·· X·kℓ

)

)

(3)

で定義される。帰無仮説のもとで，検定統計量 Cϕ , Cϕ の極限分布はそれぞれ，自由度 JKL − J − K − L + 2，自由度 (J − 1)(KL − 1) のカイ二乗分布である。
４次元の J × K × L × M 分割表の多項分布モデルにおいて，以下の仮説を考える。

(4)
(5)
(6)
(7)

(4)

H0 : pjkℓm = pj··· p·k·· p··ℓ· p···m (j = 1, . . . , J, k = 1, . . . , K, ℓ = 1, . . . , L, m = 1, . . . , M )(5)
H0 : pjkℓm = pj··· p·k·· p··ℓm (j = 1, . . . , J, k = 1, . . . , K, ℓ = 1, . . . , L, m = 1, . . . , M )(6)
H0 : pjkℓm = pj·· p·kℓm (j = 1, . . . , J, k = 1, . . . , K, ℓ = 1, . . . , L, m = 1, . . . , M )(7)
H0 : pjkℓm = pjk·· p··ℓm (j = 1, . . . , J, k = 1, . . . , K, ℓ = 1, . . . , L, m = 1, . . . , M )(4) は完全独立性の仮説，(5) は１要因対１要因対２要因の独立性の仮説，(6) は１要因対
(4)
(5)
(6)
(7)
３要因の独立性の仮説，(7) は２要因対２要因の独立性の仮説である。H0 , H0 , H0 , H0の各々の帰無仮説を検定するための検定統計量の，n → ∞ としたときの極限分布はそれ
ぞれ，自由度 JKLM −J −K −L−M +3，JKLM −J −K −LM +2，(J −1)(KLM −1)，
(JK − 1)(LM − 1) をもつカイ二乗分布である。
５次元以上の分割表においても同様に，完全独立性の仮説，および，種々の多要因
間独立性の仮説が考えられる。

3. 分割表の独立性モデルの間の関係
３次元分割表のセル確率 pjkℓ (j = 1, . . . , J, k = 1, . . . , K, ℓ = 1, . . . , L) に対して，
K × L の長方形に並んでいるセル確率をすべて１列に並べ直し，添え字を付け直すことによって得られる新たなセル確率 pji (j = 1, . . . , J, i = 1, . . . , KL) を考えると，３次元の J × K × L 分割表における１要因対２要因の独立性モデル (3) は，２次元のJ × KL 分割表における独立性モデルに帰着する。これを一般化すると，多次元分割表における r 個の要因グループ間の独立性モデルは，r 次元分割表の完全独立性モデ
ルに帰着する。
この関係に着目し，Kobe et al. [1] は，２次元分割表における独立性検定についてのTaneichi and Sekiya [2] の結果を基に，３次元分割表の２要因対１要因の独立性検定の改良をおこなった。同様に，多次元分割表における多要因間独立性検定は，より低次元の分割表における完全独立性検定の改良理論を基に，容易に改良することができる。

参考文献[1] T. Kobe, N. Taneichi and Y. Sekiya, Improved transformed statistics for the test of one
factor independence from the other two in an r × s × t contingency table. J. Japan Statist.
Soc., 45(1), (2015), 77–98.
[2] N. Taneichi and Y. Sekiya, Improved transformed statistics for the test of independence
in r × s contingency tables. J. Multivariate Anal., 98, (2007), 1630–1657.

単一強スパイク固有値モデルに対する高次元平均ベクトルの 2 標本検定
石井 晶
矢田 和善
青嶋 誠

(東京理科大学理工学部)(筑波大学数理物質系数学域)(筑波大学数理物質系数学域)

1. はじめに母集団が 2 つあり，各母集団 πi は，平均に p 次のベクトル μi ，共分散行列に p 次の正定値対称行列 Σi をもつとする．各母集団から，ni 個の p 次元データベクトル xi1 , ..., xini を無作為に抽出し，p × ni データ行列 X i = [xi1 , ..., xini ] を構成する．ここで，p > ni である．Σi の固有値を λi(1) ≥ · · · ≥ λi(p) (> 0) とし，適当な直交行列 H i = [hi(1) , ..., hi(p) ] で
Σi を Σi = H i Λi H Ti , Λi = diag(λi(1) , ..., λi(p) ) と分解する．各母集団において，標本共分散行列を S i = (ni −1)−1 (X i −X i )(X i −X i )T と定義する．ここで，X i = [xi , ..., xi ]，ni
xi = n−1i
j=1 xij である．Aoshima and Yata [3] は，高次元データの固有値モデルを次の 2 つに分類した．1 つ目は，強スパイク固有値モデル (strongly spiked eigenvalue (SSE)モデル) とよばれるもので，以下のように定義される．

lim inf

 λ2

p→∞

i(1)
tr(Σ2i )



>0

(1)

2 つ目は，弱スパイク固有値モデル (non-SSE モデル) とよばれるもので，以下のように定義される．
λ2i(1)→ 0 (p → ∞)(2)
tr(Σ2i )本講演は，典型的な SSE モデル (1) として，次の固有値モデルを考える．
p2
s=2 λi(s)(A-i)= o(1) (p → ∞)λ2i(1)

Ishii et al. [6] は (A-i) を「単一強スパイク固有値モデル (単一 SSE モデル)」とよんで，「p → ∞
単一 SSE モデル (A-i) のもとで高次元共分散構造の検定を与えた．本講演では，だが，n1 , n2 は固定」の枠組みで高次元小標本漸近理論を展開し，単一 SSE モデル (A-i)
における新たな 2 標本検定とその漸近的性質を論じた．

2. 単一 SSE モデルにおける 2 標本検定
次の仮説を考える．

H0 : μ1 = μ2

vs.

H1 : μ1 = μ2(3)

高次元データの 2 標本検定に対して，次の検定統計量はスタンダードである．
2

Tn = x1 − x2  −

2


n−1
i tr(S i )i=1

Bai and Saranadasa (1996, Stat. Sin.)，Chen and Qin (2010, Ann. Stat.)， Aoshima and Yata
[1, 2] は，Tn の漸近正規性を non-SSE モデル (2) のもとで論じ，それぞれ異なる高次元

2 標本検定を与えた．一方，SSE モデル (1) においては，Tn の漸近正規性は崩れることが，Aoshima and Yata [3] により示された．Ishii [4] は，単一 SSE モデル (A-i) のもとでTn を再評価し，高次元 2 標本検定を与えた．Ishii [4] は，第 1 固有空間について，次を
仮定した．
λ1(1)
(C-i)
= 1 + o(1), |hT1(1) h2(1) | = 1 + o(1) (p → ∞)λ2(1)実際のデータに対する仮定 (C-i) の妥当性は，Ishii et al. [5] が与えた第 1 固有空間の同等性検定を用いて確認することができる．Yata and Aoshima [7] のノイズ掃き出し法による λi(1) の推定量を λ̃i(1) とおくと，以下の検定統計量が考えられる．
Tn + 2i=1 λ̃i(1) /ni
Fa = un (n1 + n2 − 2) 2
i=1 (ni − 1)λ̃i(1)
ここで，un = (1/n1 + 1/n2 )−1 とする．Fa の漸近帰無分布は次のようになる．
定理 2.1 (Ishii [4]). 各母集団に (A-i) を仮定する．また，(C-i) を仮定する．そのとき，いくつかの正則条件と H0 のもとで，次が成り立つ．

Fa ⇒ F1,ν

(p → ∞ だが，n1 , n2 は固定)

ここで，ν = n1 + n2 − 2，Fν1 ,ν2 は自由度 (ν1 , ν2 ) の F 分布に従う確率変数を表す．
与えられる α ∈ (0, 1/2) に対して，仮説 (3) について，次の検定方式を考える．

Fa ≥ F1,ν (α) のとき，H0 を棄却する(4)

ここで，Fν1 ,ν2 (α) は Fν1 ,ν2 の上側 α 点を表す．Ishii [4] は，検定方式 (4) の第 1 種の過誤と検出力に関する理論的な結果を与え，検定方式 (4) の漸近的な性能を保証した．しかし，実際のデータ解析においては第 1 固有空間の同等性の仮定 (C-i) は必ずしも成立しない．そこで，本講演では，(C-i) を仮定できない場合に，新たな高次元 2 標本検定を
与た．

参考文献
[1] Aoshima, M. and Yata, K. (2011). Two-stage procedures for high-dimensional data, Sequential
Analysis (Editor’s special invited paper), 30, 356-399.
[2] Aoshima, M. and Yata, K. (2015). Asymptotic normality for inference on multisample, highdimensional mean vectors under mild conditions, Methodology and Computing in Applied
Probability, 17, 419-439.
[3] Aoshima, M. and Yata, K. (2018). Two-sample tests for high-dimension, strongly spiked eigenvalue models, Statistica Sinica, 28, 43-62.
[4] Ishii, A. (2017). A two-sample test for high-dimension, low-sample-size data under the strongly
spiked eigenvalue model, Hiroshima Mathematical Journal, 47, 273-288.
[5] Ishii, A., Yata, K. and Aoshima, M. (2016). Asymptotic properties of the ﬁrst principal component and equality tests of covariance matrices in high-dimension, low-sample-size context,
Journal of Statistical Planning and Inference, 170, 186-199.
[6] Ishii, A., Yata, K. and Aoshima, M. (2019). Tests for high-dimensional covariance structures,submitted.
[7] Yata, K. and Aoshima, M. (2012). Effective PCA for high-dimension, low-sample-size data withnoise reduction via geometric representations, Journal of Multivariate Analysis, 105, 193-215.

13 の癌腫の 78 個の Microarray による LSD の確認と Method1 による SM と BGS 分割の評価― Big データ時代における LSD の重要な 2 つの特徴 ―新村秀一

成蹊大学名誉教授

1999 年から 2004 年に米国の 6 医学研究が医学誌に論文を発表し、研究に用いたデータを公開している。これ
らの医学研究では癌と健常、あるいは異なった 2 種の癌症例を用いて、「癌の遺伝子診断に役に立つ癌の遺伝子の組」を特定することを目的としている。筆者は、2015 年 10 月 25 日に Jeffery らが、これらの 6 種のマイクロアレイを HP に公開していることを知った。そこで 10 月 28 日にそれらをダウンロードし、僅か 54 日間で判別理論の応用問題として簡単に次のように解決した。
1) MNM 基準による Revised IP-OLDF(RIP)で 6 個のデータを判別し、MNM=0 であり、2 群が遺伝子空間で線形分離していることが分かった。線形分離可能なデータ(LSD)である遺伝子空間とその部分空間を Matryoshka と呼ぶ。スイス銀行紙幣データ等の LSD 判別研究で、LSD は「Matryoshka 構造」を持つことが分かった。すなわ
ち LSD はその中に小さな Matryoshka を最小次元の Matryoshka である BGS まで入れ子状に含んでいる。iPS
細胞研究で、山中らは 2 万個以上といわれる遺伝子から、Rats の ES 細胞で活性化している 24 個の遺伝子で万能細胞の塊ができることを見つけた。そして重回帰分析の変数減少法に似た方法で、24 個から 4 個の山中遺伝子を特定した。4 個から 1 個を省くと万能細胞ができないことと BGS から 1 個省くと MNM が 1 以上になることは同じ概念である。高次元データから癌の遺伝子を特定することは NP-hard ではあるが、このような方法によるアプローチが必要である。ただし偏差平方和の最小化でなく MNM 基準で行う必要がある。
2) もし BGS がｋ組あれば、簡単にｋ組の小標本に分割でき高次元の呪いから解放される。このために、数理計画法ソフトの LINGO で BGS を求める Program4 を開発した。しかし、遺伝子数 p が大きいと計算時間がかかるので、Matryoshka Feature Selection Method(Method2)で Small Matryoshka（SM）を見つける Program3 を開発した。
すなわち高次元であろうが小標本であろうがデータが LSD であれば、1)Matryoshka 構造を持ち、2)複数個の排他的な BGS(あるいは SM)と MNM が 1 以上のその他の部分空間に分割できる。そして 6 個の高次元データから求まった小標本を統計分析し、マイクロアレイによる癌の遺伝子診断を世界で初めて提案した。
しかし用いたデータは 2004 年までのものであり、癌腫が限られている。また僅か 6 種類のデータが LSD であっても、他の多くのデータが LSD でないと、医学研究への適用が限られる。そこで Bruno らが開発した Curated
Microarray Database (CuMiDa)が最適と判断した。そしてこの DB から 13 種類の異なる癌の 78 個のデータをダ
ウンロードした。そして、2 群判別で全てが MNM=0 であることが分かった。すなわち、すでに分析した 6 個のデータに加えて、新しい 77 個のデータでも全て LSD である。これで高次元データの呪いから解放され、SM あるいは BGS を統計分析し癌の遺伝子診断が行えるという筆者の主張の正しさと汎用性が確認できた。CuMiDa では、今後 RNA-seq などの遺伝子データの搭載も検討している。今回の研究から、恐らく RNA-seq でも LSD であることが確認できると考えている。もし確認できない場合、マイクロアレイより RNA-seq の方が優れているという研究者もいるが、RNA-seq が癌の遺伝子診断に役に立たないことを示す。以上の結果は、マイクロアレイやその他の遺伝子データによる癌の遺伝子診断を今後大きく発展させることが期待できる。
また医学研究者が、新規の研究でマイクロアレイや RNA-seq を集めた場合、最初に LSD であることを確認すべきである。LSD であれば、これまでの多くの研究方法は全て無用になり、SM と BGS に分割して遺伝子診断を行えばよい。筆者の方法で重要な信号はもれなく列挙できる。しかし数多くある SM と BGS の遺伝子診断の重要度を検討する必要がある。このために、すでに開発した 100 重交差検証法(Method1)で、検証標本の平均誤分類確率
M2 の小さな順に遺伝子診断の重要度のランク付けを行った。

ベイズモデルの母数：事前分布と 
柳 本 武 美 統計数理研究所

 序
標本分布の族を Ü         とする。ベイズ推定では、事前密度    を仮定して推論を行う。母数  の推定は事後平均      
で行うことが通常である。但し    は関数   の密度  の下での期待値とする。推定量を得るための母数を と呼ぶ。この場合は  が になる。ところでどのよう を選ぶと良いだろうか。確率モデルを仮定する場合には、母数の意味が取りやすいようにに選択する。特に     な事前分布を仮定するためには、意味が取りやすい母数を選ぶ必要がある。しかし、母数を推定する段階では母数の意味について考えるより、リスクあるいはベイズ
リスクを小さくするように選ぶ必要がある。この違いは従来等閑視されてきたように思われる。別の理由として、応用的側面から母数が決まるだろうと考えることにある。この場合も応用的側面と推定量の効率をどう捉えるのかを明瞭にする必要がある。

 損失と 

の選択

理論統計学の基本的枠組みでは、推定量はリスクによって評価される。リスクを小さくする推定   を定めて議論するのがよい。広く用いられ量が良い推定量である。そうだとすると、損失 ¾
る損失は二乗損失     と    損失

 

 



  Ý Ý  Ý 



 !

  とである。ここで Ý   の Ý は予測分布を表す Ü と同じ次元の変数である。この損失を        
 も使われる。    損失一対一に対応する
書いたとき、その双対形  
母数の選択に依存しないが、二乗損失は依存する。しかも、二乗損失は母数が多次元の場合などには妥当性が理解しにくい問題がある。

 指数分布族の場合




$ 
指数分泌族  Ü   "#   Ü   
をかんがえる。母数  は共役¼
  の事後平均         Ü で推定する。対母数と呼ばれる。通常は、平均母数  ¼ ½¼

応する推定量を     と書く。この推定量の #  予測子Ý   は 別の予測
子  Ý Ü   Ý       より損失     の下でのベイズリスクが大きくなる。

この困難は



を共役母数  とすることにより解消される。 #  予測子は

Ý Ü



"#





      



% !

   の下で最適になる。この周辺は &   ' となる。損失  
 +, , --() -!! で調べられた。

次の等式はよく知られている。 Ü   Ü    に対して 成り立つ。一方、すべての Ü に対して






Ü Ü
   ÜÜ  



 





Ý Ü



!((() * Ý   Ý  がすべての Ü と Ý    Ü

が成り立つ。このことは、. は損失を小さくし過ぎることを示している。

    分布の場合
平均が  分散が !  ¾ である # 分布の密度関数は

    





"#



   / !で与えられる。この分布は指数分布族に含まれない。歴史的には 標本中央値と . が一致する
ことから、# !00/) 1 !(!! らの関心を惹いた。その結果この推定量の改善は ' 2 !(33 を始め多くの研究結果がある。この分布と推定量についての成果は  4 ら --!に詳しく記述されている。
そこで、ベイズ推定における
 を考察する。最初に考えつく    には問題が多い。 ¾¾と !  の区別が難しいし、二乗誤差    5     では二つの項の和を取ることの意味が取りにくい。
さて    損失は
 "#     5      ! 5
 
 !




と書けるが扱いにくい。そこで、    で近似すると
                !
¾
    5



 
 !



/ 

/ %

となる。この損失の下でのベイズリスクを最小にする
 を求めると     となる。
提案法は、損失 / % の下でベイズリスクを最小にする意味での最適性が満たされている。この事
実から、他の直感的あるいは強い制約条件の下での推定量よりも良い性能を示すと期待される。既存の推定量として! 最尤推定量同時に標本中央値
 最良線形不偏推定量 ' 2) !(33% 
 を )  とした事後平均が挙げられる。この中で上でも述べたように、 ! は  に劣るので実質的には二つとの性能比較になる。そこで     事前分布を仮定した下で性能比較を実行して、良好な結果を得つつある。
の推定量は線形不偏に制約する理由が分かりにくい。結果として性能が劣る。% のベイズ推定量は
 の選択の理由が直感的である。実際、  の場合に  の推定量が存在しなくなる。
実際に数値実験的に比較を行うとほぼ期待通りの結果が得られつつある。興味深いことに、! 
が !  の推定量としてのバイアスが小さいことが観察された。この事実に対応する不偏性は、正規分布では解析的に示される。

 ! & ) 6

 ') & 9

 ', ,) .

.

 ' !((0

7

!(((

   ) ) 380((03(:% '  2) ;!(33     )  )%%: 8  4) < ) 4=  ) >  ?  )  --!         <#   3 # ) ? < !00/         ) ) 3!383 0* ) >  +, , ) > --(    !) ) %-3/%-08 : *  )
>  +, , ) > -!!    !) ) !((---    )  ) /:8:  /  1) 6 .

!(!!

On Gaussian semiparametric estimation for
two-dimensional intrinsic stationary random fields
Yoshihiro Yajima
Tohoku University
We propose a Gaussian semiparametric estimator for semiparametric models of twodimensional intrinsic stationary random fields (ISRFs) observed on a regular gridand derive its asymptotic properties. Originally this estimator is an approximate
likelihood estimator in a frequency domain for long memory models of stationary
and nonstationary time series(Robinson(1995); Velasco(1999)). We apply it to twodimensional ISRFs. These ISRFs include a fractional Brownian field, which is a
Gaussian random field and is used to model many physical processes in space. The
estimator is consistent and has the limiting normal distribution as the sample size
goes to infinity. We conduct a computational simulation to compare the performanceof it with those of diﬀerent estimators.
Let {X(s) : s ∈ Rd } be a random field. If {X(s)} satisfies that for any fixedh(∈ Rd ), the increment Zh (s) = X(s + h) − X(s) is a stationary random field,{X(s)} is called an ISRF. Then {X(s)} is characterized by
E(X(s + h) − X(s)) = 0,
Var(X(s + h) − X(s)) = 2γ(h),where 2γ(h) is the variogram function. Hereafter we also assume that X(0) = 0.
Under some assumptions, the variogram function is expressed by∫
1 − cos((λ, h))2γ(h) =
g(λ)dλ,d(2π)dRwhere g(λ) satisfies∫∥ λ ∥2g(λ)dλ < ∞.
2
Rd 1+ ∥ λ ∥respectively.
Hereafter we assume that d = 2 and {X(s)} is a two-dimensional Gaussian ISRF.
Then we also denote λ by (λ1 , λ2 ) and the spectral density function g(λ) by g(λ1 , λ2 )respectively.
We introduce the following assumptions.
1Assumption 1 g(λ1 , λ2 ) is expressed byg(λ1 , λ2 ) =∥ λ ∥−2H−2 go (λ1 , λ2 ),0 < H < 1,

where go (λ1 , λ2 ) is a nonnegative with go (0, 0) > 0, symmetric , go (λ1 , λ2 ) =go (−λ1 , −λ2 ), twice continuously diﬀerentiable function for −∞ < λ1 , λ2 < ∞ and isbounded with bounded first and second order partial derivatives.
Then we have the following asymptotic properties of the estimator.
Theorem 1 Let Ĥn be the Gaussian semiparametric estimator where n is the sample
size. Under Assumption 1 and additional conditions, Ĥn converges to H0 in probability as n → ∞. Moreover for p ≥ 2, m1/2 (Ĥn −H0 ) converges to N (0, 1) in distribution
as n → ∞. where m is the number of the discrete Fourier transforms used to construct
the estimator and p is the order of data tapers.
We conducted some computational simulations and compared the performance of ourestimator with that of the estimator proposed by Zhu and Stein(2002) which are basedon empirical variograms. Then our estimator is less biased than the estimator of Zhu
and Stein if go (λ1 , λ2 ) is not a constant.
References
Robinson,P.M. (1995). Gaussian semiparametric estimation of long range dependence.
Ann. Statist. 23 1630-1661.
Velasco,C. (1999). Gaussian semiparametric estimation of non-stationary time series.
J. Time Ser. Anal. 20 87-127.
Zhu,Z. and Stein,M.L. (2002). Parametric estimation for fractional Brownian surfaces.
Statist. Sinica12 863-883.

2

再生核ヒルベルト空間における射影のモーメントによる
二標本検定
千葉大・融合理工学府 牧草 夏実

1 はじめにi.i.d.

i.i.d.

P ，Q をヒルベルト空間 H 上の確率分布とするとき，二標本 X1 , . . . , Xn ∼ P ，Y1 , . . . , Ym ∼Q に基づく検定
帰無仮説 H0 : P = Q vs. 対立仮説 H1 : P ̸= Q

を考える．ユークリッド空間での二標本検定はすでに様々な検定方法が議論されているが，ヒルベルト空間に値をとる確率変数に対する二標本検定を考えることで，高次元データに対する二標本検定の議論を与える．高次元データに対するアプローチとして，カーネル法という方法がある．カー
ネル法を用いた二標本検定として，Maximum Mean Discrepancy(MMD) に基づく二標本検定が[1] によりすでに議論されているが，この MMD と同様の考え方により，2 次のカーネルの MMDによる二標本検定を考える．

2 再生核ヒルベルト空間でのモーメントの定義
確率変数 X ∼ P ，Y ∼ Q を正定値カーネル k によって，この k に対応する再生核ヒルベルト空間 H(k) 上に，それぞれ k(·, X)，k(·, Y ) により変換を行う．このとき，この k(·, X)，k(·, Y ) の平
均まわりの 2 次モーメント Σk (P )，Σk (Q) は，それぞれヒルベルト空間 H(k)⊗2 = H(k) ⊗ H(k)での期待値 Σk (P ) = EX∼P [(k(·, X) − µ(P ))⊗2 ]，Σk (Q) = EY ∼Q [(k(·, Y ) − µ(Q))⊗2 ] によって定められている．ここで，µ(P )，µ(Q) は k(·, X) の 1 次モーメント µ(P ) = EX∼P [k(·, X)]，µ(Q) = EY ∼Q [k(·, Y )] であり，⊗ はテンソル積を表しており，任意の f ∈ H(k) に対し，f ⊗2 =f ⊗ f = ⟨f, ·⟩H(k) f である．

3 検定統計量の構築
この k(·, X) と k(·, Y ) の f ∈ H(k) への射影のモーメント差

sup||f ||H(k) =1| ⟨f, µ(P ) − µ(Q)⟩H(k) | = ||µ(P ) − µ(Q)||H(k)により，2 つの分布の違いを測るのが，Maximum Mean Discrepancy (MMD) と呼ばれるものである．同様の考え方により，(k(·, X) − µ(P ))⊗2 と (k(·, Y ) − µ(Q))⊗2 の A ∈ H(k)⊗2 への射影のモーメント差

sup||A||H(k)⊗2 =1| ⟨A, Σk (P ) − Σk (Q)⟩H(k)⊗2 | = ||Σk (P ) − Σk (Q)||H(k)⊗2

により 2 つの分布の違いを測る．これは，MMD のようなある種の分布の違いを測っている．この違い ∥Σk (P ) − Σk (Q)∥2H(k)⊗2 は

b k (P ) − Σb k (Q)
Tb2 = Σ2
H(k)⊗2によって推定することができる．ただし，
n
n
1∑
1∑
⊗2
b
Σk (P ) =(k(·, Xi ) − µb(P )) , µb(P ) =
k(·, Xi ),
n i=1
n i=1
m
m
∑
1 ∑
b k (Q) = 1Σ(k(·, Yi ) − µb(Q))⊗2 , µb(Q) =
k(·, Yi )
m i=1
m i=1

である．
本発表では，この検定統計量 Tb2 の漸近挙動について報告を行った．特に，帰無仮説 H0 : P = Qのもとで，退化 V 統計量の結果に帰着させることで，Tb2 の漸近分布が，独立な自由度１の χ2 分
布の重み付きつき無限和の形で得られること ([2] 参照)，その重みをデータに基づいて推定する方法について報告を行った．また，対立仮説 H1 : P ̸= Q のもとでの Tb2 の漸近分布が平均 0 の正規分布になっていることについて報告を行った．

参考文献[1] Gretton, A., Borgwardt, K. M., Rasch, M. J., Schölkopf, B. and Smola, A. (2007). A
kernel method for the two sample problem. Advances in Neural Information Procesing
Systems, 19 of MIT Press, Cambridge.
[2] Serfling, R. J. (1980). Approximation theorems of mathematical statistics. Wiley, NewYork.

非負データに対する再帰的な密度推定について
柿沢佳秀 (北大経済)
１．はじめに
位置尺度型のカーネル推定量 (Rosenblatt(1956)，Parzen(1962)) は，畳み込み近似が背後にあって，supp(f ) = R のときに望ましい性質を持つことが知られているが，非負または区間データの密度推定においては，境界付近で一致性がない．そのために境界バイアス問題の対処も議論されている (リノーマライゼーション法，リフレクション法，
一般化ジャックナイフ法は Jones(1993) を参照のこと; なお，このカテゴリーに，局所
多項式法による密度推定法も含まれる)．Chen(1999,2000) 以降，非対称カーネル密度推定法が境界バイアス問題の対処の１つとして脚光を浴び，ベータ，ガンマ，逆ガウス(IG)，相反逆ガウス (RIG)，バーンバウム・サンダース (BS)，対数正規，逆ガンマ等
個別に研究が進められ，現在のところ，次の (1)–(3) のように分類される (積カーネルを使えば，多次元データも扱えるので，簡単のため１次元データとする):(1) Amoroso カーネル密度推定量族 (Igarashi and Kakizawa(2018))．
(2) IG と RIG との -混合 (MIG) カーネル密度推定量族 (Igarashi and Kakizawa(2014))，
及び，加重対数正規カーネル密度推定量族 (Igarashi (2016))．また，これらはN(0, 1) ベースであったから，対称分布ベースの q-MIG カーネル密度推定量族(Kakizawa(2018)) へも拡張されている (q → 0 としたサブクラスが，対数対称
カーネル密度推定量族である)．
(2 ) (2) の対称分布ベースの q-MIG カーネル密度推定量族の類似だが，歪分布ベースの q-BS カーネル密度推定量族 (柿沢; 2017 年度統計関連学会連合大会) へも拡張されている; なお，対称分布を外した為，IG と RIG 相当の拡張は定義できない．
(3) インデックスが ν = ±1/2 の一般化逆ガウス (GIG) カーネル密度推定量 (なお,ν = −1/2(IG), 1/2(RIG) を (2) に分類した; Igarashi and Kakizawa(2014) は GIGの代わりの別名 “変形ベッセル”を採用)．
さて，非対称カーネル密度推定量とは，無作為標本 X1 , . . . , Xn ∼ iid f が与えられ，付随する平滑化パラメータを bn  0 (n → ∞) とおくとき，点 x ∈ supp(f ) で

fbn (x) =n1
k(Xi ; bn , x)n i=1のような形式で定義されるが，これは “非再帰的”である．すなわち，仮に追加データがあった場合に，計算は最初からやり直さねばいけない．他方，再帰的ならば，計算結果を更新するので，追加データを効率的に扱えるメリットがある．この種の話題は，
Rosenblatt–Parzen カーネル法が流行した頃に Yamato(1971) を初めとして，しばしば議論されていた (磯貝 (1993) を参照) が，当時，Chen(1999) 以降の密度推定量は対象となっておらず，本報告では，再帰的な非対称カーネル法を考察した．
２．再帰的な非対称カーネル密度推定量
正数列 γn (ステップサイズと呼ぶ) を選び，f0 (x) = 0 として, i = 1, 2, . . . , n に対し
fi (x) = (1 − γi )fi−1 (x) + γi k(Xi ; bi , x) のように再帰的に更新する．任意の g についてf (x) − g = Π
n

n

n

i=1


Π−1
i γi {k(Xi ; bi , x) − g} + Πn {f0 (x) − g}を示すのは容易である．ここに，Πn = ni=1 (1 − γi )．
以下，Galambos and Seneta(1973) の定義に従い，正数列 hn について (hn ) ∈ GS(d)であるとは，limn→∞ {n(1 − hn−1 /hn )} = d ∈ R を満たすこととする (Bojanic and
Seneta(1973) も参照)．
主結果 [supp(f ) = [0, ∞) を念頭にしたが，[0, 1] でも同様な結果になる] (γn ) ∈ GS(−γ)，(bn ) ∈ GS(−c)，limn→∞ (nγn )−1 = ξ を仮定する．ただし，γ, c, ξ について(A) x > 0 ならば，0 < γ ≤ 1，0 < c < 2γ ，0 ≤ ξ < max{1/c, 2/(γ − c/2)}(B) x = 0 ならば，0 < γ ≤ 1，0 < c < γ ，0 ≤ ξ < max{1/c, 2/(γ − c)}
とする．
B(x)
(i) (1) x > 0 に対して 0 < c ≤ (2/5)γ ならば Bias[fn (x)] ≈ bn1−cξ，2γ > c > (2/5)γならば Bias[fn (x)] = o((γn b−1/2)1/2 )．
nbn δ(2) 0 < c ≤ (1/3)γ ならば Bias[fn (0)] ≈ 1−cξ，γ > c > (1/3)γ ならば Bias[fn (0)] =1/2o((γn b−1
)．
n )
−1/2

V (x)
n bn
(ii) (1) x > 0 に対して 2γ > c ≥ (2/5)γ ならば V [fn (x)] ≈ γ2−(γ−c/2)ξ
，0 < c < (2/5)γならば V [fn (x)] = o(b2n )．
γn b−1n v(2) γ > c ≥ (1/3)γ ならば V [fn (0)] ≈ 2−(γ−c)ξ
，0 < c < (1/3)γ ならば V [fn (0)] =o(b2n )．
⎧⎪
if 0 < c < (2/5)γ,⎪⎨ GS(−2c) ,(iii)(1) x > 0 に対して AMSE[fn (x)] ∈ GS(−(4/5)γ) ,if c = (2/5)γ,⎪⎪⎩ GS(−(γ − c/2)) , if 2γ > c > (2/5)γすなわち，AMSE[fn (x)] を最小化するならば，c = (2/5)γ でなければいけない; このとき，AMSE[fn (x)] ∈ GS(−(4/5)γ), γ ∈ (0, 1] なので，γ = 1 でなければいけない．
また，γ = 1，c = (2/5)γ = 2/5 のとき

V (x)b2n B 2 (x)γn b−1/2n+{1 − (2/5)ξ}2
2 − (4/5)ξ
2/5
5
γn2
≥ 4/5
V 4/5 (x){B 2 (x)}1/5 .
3
4
4{1 − (2/5)ξ}AMSE[fn (x)] =なお，γn = 1/(ξn) (ξ ∈ (0, 5/2)) とすれば，

n−2
n−2
γn2
=
≥
{1 − (2/5)ξ}3
ξ 2 {1 − (2/5)ξ}3(3/5)3(等号は ξ = 1 に限る).

(2) AMSE[fn (0)] を最小化する場合, 同様にして γ = 1，c = (1/3)γ = 1/3 となり，
γn = 1/(ξn) (ξ ∈ (0, 3)) とすれば，ξ = 1 に至る．
(iv) AMISE については (iii-1) の積分版に相当の結果になる (B 2 (x), V (x) を，それぞれ2
I B = 0∞ B 2 (x) dx, I V = 0∞ V (x) dx とせよ)．
(v) 漸近正規性，強一致性も成立．
注意．AMSE あるいは AMISE の最小化で考えると，再帰的推定量は非再帰的推定量に劣る (上記のような最小値を，非再帰版と比較すれば分かる)．他方，過小平滑化をすれば，分散 (分散積分) の主要項の最小化での比較になり，このときは再帰的推定量が良い．

ノンパラメトリックな極値解析法の漸近理論
鹿児島大学大学院理工学研究科 吉田拓真

1

はじめに

降水量, 人口推移, 金融リスクに関するデータでは, 平均よりも最大値や最小値, もしくはそれに近い分位点に興味がある. 別の言い方をすれば, 分布の中心やピークではなく, 裾に興味がある. 分布の裾は大きくわけて３つのタイプに分けられる. ひとつは, 正規分布などに代表される裾確率が標準的に小さい分布 (light tail), もうひとつは一様分布, ベータ分布に代表される裾が有限である分布 (short tail), そして t-分布などいわゆる正規分布よりも重い裾を持つ分布 (heavy tail) である.
本講演では, heavy tail を持つ分布の裾の情報を推測するための方法と理論について述べた.

2モデル

一般に, Y ∼F が heavy tail するとは, y が十分大きいところでは

P (Y > y) = 1 − F (y) = y −1/γ L(y)(1)という構造を持つことをいう. ただし, γ > 0 は Extreme value index (EVI) と呼ばれるパラメー
タで, L(y) は L(y) → c > 0, y → ∞ を満たすノンパラメトリック関数である. この分布は Paretotype-tailed model と呼ばれており, t-分布, Burr 分布など多くの裾の重い分布を含む. 実際に, (1)は分布の裾が次数 1/γ の多項式オーダーで 0 に収束していくことのみを仮定しているにすぎない.
このモデルを考える上で最も重要なのは, γ を推定することといえる. EVI は分布の裾の重さを表すパラメータで, 例えば t-分布の自由度を ν とすると, ν = 1/γ と対応する. すると, γ = 1 でコー
シー分布の裾と同等であり, γ = 0.5 で分散を持たないほど重い分布であることがわかる.
さて, 仮に γ の推定量 γ̂ を得たとすると, γ̂ − γ の収束レートが気になる. しかし, それを数学的に示そうとすると結局 L(y) → c の収束オーダーが必要となってくる. よって, γ̂ の数学的性
質解明のため, L はしばしば, パラメータ A > 0, ρ < 0 とある関数 u を用いてL(y) = 1 + u(y), |u(y)| ≤ Ay ρ/γを満たすと仮定される. 特に,ρ は 2 次オーダーと呼ばれ, γ̂ − γ の収束レートを決定する重要なパラメータである. この仮定をしてなお多くの分布が (1) に含まれることが知られている (see, Hall
and Welsh 1984). 実際に γ̂ を推定する方法としては最尤法があげられるが, ((1) から導かれる密∑度関数を f (y|γ) と書く) そのまま ni=1 log f (yi |γ) を最大化する γ を求めても良い推定量とはならない. なぜならば, (1) はあくまで裾の情報であって, 平均周りに分布するデータなどは (1) に従っている保証はないからである. ではどうするかというと, ある閾値 w を超えるデータのみで推定を行う. このとき, 考える分布は()
Y
P> y|Y > w = 1 − Fw (y) = y −1/γ L̃w (y)wと, 再び Pareto-type モデルになるため, このモデルに基づく最尤法などで γ を推定する. ただし, Fw は Y > w の下での Y /w の分布関数, L̃ は L̃w (y) → c̃ > 0(y → ∞) を満たす関数であ∑る. 最尤法は ni=1 fw (yi /w|γ)I(yi > w) を最大化させる. ここで, fw は Fw に対応する密度関数であり, I は指示関数である. このように, 閾値を超えるデータのみを用いて推定を行う方法はPeak-over-threshold(POT) 法と呼ばれ, EVI の一般的な推定方法となっている. 閾値 w の選び方に関する研究論文も数多く発表されている.
説明変数 x を用いた回帰分析では, 基本的には上記のモデルすべてが x の関数となり, Paretotype-tailed regression model
P (Y > y|X = x) = 1 − F (y|x) = y −1/γ(x) {c(x) + u(y|x)}, |u(u|x)| ≤ A(x)|y|ρ(x)/γ(x)(2)というモデルに基づいて γ(x) の推定を行う. ここで, 任意の x ∈ X に対して, γ(x) > 0, A(x) > 0,
ρ(x) < 0 である. この回帰分析の枠組みでは, γ(x) のパラメトリック推定・ノンパラメトリック推定が主な目的となる. 一方で, EVI が x に依存しない: γ = γ(x) という仮定や仮説検定などは多くの研究者が議論している. EVI の推定に POT 法を用いるときは, 閾値も x に依存する (w(x) と書く). そして, γ̂(x) − γ(x) の n → ∞ の下での収束レートを調べることが重要となる. ここで, nはデータ数である. 関連する先行研究は多いが, そこでは基本的に閾値関数 w(x) は固定した下で行っている. しかし, 閾値関数 w(x) は極値に位置するデータをピックアップするためにあることから, データ数が変化するとそれに対応して変化するのが自然である. もし, n → ∞ で, w(x) を固定したままだと, 相対的に極値ではない, 言い換えると, (2) に従っていないデータを推定に含むことになる. つまり, n → ∞ の下では, w(x) = wn (x) → ∞ と仮定するのが自然である. 実際に,シミュレーションにおいて何らかの閾値選択 (クロスバリデーションや QQ-plot を用いた方法など) を行うとそのような傾向が見られる. 上記で述べた EVI 推定の先行研究などは Dey and Yan(2016) に詳しくまとめられている.

3

主結果

本研究では回帰分析の文脈で関数 γ(x) をノンパラメトリック法で推定し, supx∈X |γ̂(x) − γ(x)| の収束レートを導く. 真の γ(x) を m 次微分が連続な関数とし, γ̂(x) を最尤法をベースにしたノンパ
ラメトリック推定量とする. このとき, 適切な条件の下で
()m/(2m+1) ()−1/γ(x)log{supx∈X |wn (x)|n}ρ(x)/γ(x)sup |γ̂(x) − γ(x)| = OP
+ O sup |wn (x)|supx∈X |wn (x)−1/γ(x) |nx∈Xx∈Xが成り立つ. 第 1 項はノンパラメトリック法によるバイアス-分散で, 第 2 項は (2) のモデルバイア
スである. 特に, 閾値関数として, X = x を与えた下での Y の分布関数 F (·|x) に対応する分位点
wn (x) = F −1 (τn |x), τn → 1(n → ∞) とすると,
({
} m )
(
)
log{(1 − τn )n} 2m+1
sup |γ̂(x) − γ(x)| = OP
+ O (1 − τn )−ρ
(3)(1 − τn )nx∈Xとなることがわかった. ただし, ρ = supx∈X ρ(x) である. (3) はノンパラメトリック分位点回帰をベースにした EVI 推定量の収束レート (Yoshida 2019) と一致していること, 回帰ではない場合の自然な展開になっていることを報告した. また, 分位点回帰への応用について報告した.

References
1.Dey,D.K. and Yan,J. (2016). Extreme value modeling and risk analysis: Methods and applications. Chapman & Hall/CRC.
2.Hall, P. and Welsh, A.H. (1984). Best attainable rates of convergence for estimates of parameters of regular variation. Ann. Statist. 12 1079–1084.
3.Yoshida,T. (2019). Simultaneous confidence bands for extremal quantile regression with
splines. Extremes. Accepted. DOI:10.1007/s10687-019-00360-4

ノンパラメトリック重回帰分析におけるLasso 型推定の改良
島根大

自然科学研究科

千葉大

1

理学研究院

松島 佑樹
内藤 貫太

はじめに

ノンパラメトリック重回帰分析において，局所線形回帰の手法と Lasso 型推定を組み合わせた方法について考える．この方法は Bertin and Lecué (2008) ですでに研究されているが，先行研究では推定したい回帰関数の値にも罰則を付けた上で推定量を計算しているため推定精度が悪くなる問題が挙げられる．この問題を改善するため，回帰関数に関する罰則項を取り除いた推定手法について考え，先行研究で用いられた推定量と我々の提案手法による推定量，これら2 つを比較した結果を報告する．

2

設定

目的変数 Y ∈ R，d 次元の説明変数 X ∈ Rd のデータ (Y1 , X1 ) , . . . , (Yn , Xn ) にノンパラメト
リック重回帰モデル
Yi = f (Xi ) + ei , i = 1, . . . , ni.i.dを想定し，回帰関数 f の推定を考える．ここで，ei ∼ N (0, σ 2 ) は誤差であり，n ≪ d を仮定しておく．f の推定を行うため，カーネル関数を用いた局所重み付き線形回帰の手法を用いる．
推定を行う点を x ∈ Rd とし，モデル式を線形式

Z = Aθ∗ + ε

(1)

)T
(
T
= (A0 · · · Ad )，ε = (ε1 · · · εn )T でと変形する．ここで，Z = (Z1 · · · Zn ) ，A = ÃT1 · · · ÃTn
あり，i = 1, . . . , n に対して，
(
)
(
)T )
(
Xi − x1
Xi − xZi = αi Yi , αi = √K, Ãi = αi 1
, εi = αi ei +αi f (Xi )−ÃTi θ∗dhhnhである．また K はカーネル関数で原点対称な d 次元密度関数，h > 0 はバンド幅である．θ∗ =(θ0∗ θ1∗ · · · θd∗ )T = (f (x) h∂1 f (x) · · · h∂d f (x))T がこのモデルでの回帰パラメータベクトルとなり，∂i は第 i 変数による偏微分作用素を表す．(1) における回帰係数ベクトル θ∗ の推定を行う上で，特に f (x) の推定を考えているため，θ0∗ = f (x) の推定が重要となることに注意しておく．

3

Lasso 型推定

Bertin and Lecué (2008) では，f (x) において本質的となる変数の選択のために Lasso を使った推定minθ∈Rd+1{}∥Z − Aθ∥22 + 2λ∥θ∥1 , θ = (θ0 θ1 · · · θd )T(2)が用いられ，変数選択の一致性が示されている．しかし，この手法では，θ0 にも罰則を付けているため，f (x) の値にも制約を課した上で推定量が求められていることから，θ0∗ = f (x) の推定精度が悪くなる可能性が考えられる．この問題を改善するため，θ0 には罰則を付けない推定方法

minθ∈Rd+1{∥Z − Aθ∥22 + 2λ∥θ−0 ∥1}(3)を提案する．ここで θ−0 ∈ Rd は θ から θ0 を取り除いた部分ベクトル θ−0 = (θ1 · · · θd ) であり，λ は正則化パラメータである．
(2)，(3) の解の同値条件をそれぞれ導き比較することで，アクティブな変数とノンアクティ
ブな変数に関しての条件が変わらないことがわかった．この事実から，提案手法 (3) から得られる推定量においても (2) と同様に変数選択の一致性が導かれることを報告した．
また，理論的比較をするため，しかるべき条件を満たした説明変数 X に対し，それぞれ (2)と (3) のユニークな解として求められる 2 つの理想的推定量を紹介した．この 2 つの理想的推定量は θ0∗ を含んでいたため，実装することはできず，その意味でこれらを「理想的」と呼んだ．
n → ∞，h → 0，nhd+4 → ∞ における 2 つの理想的推定量の漸近挙動を見ることで，(2) と (3)の理論的比較を行った．その結果，(3) から得られる推定量のバイアスは O(hλ) = O(h2 ) であったが，(2) から得られる推定量のバイアスは O(λ) = O(h) となり，このことから，提案手法が推定量のバイアスを改善していることを報告した．また分散は Ruppert and Wand (1994) の結果と同じになることを述べた．
T

研究集会においてはこのような理論的結果の詳細および数値実験の結果について報告を行った．

参考文献[1] Bertin, K. and Lecué, G. (2008). Selection of variables and dimension reduction in highdimensional non-parametric regression. Electronic Journal of Statistics. 2, 1224–1241.
[2] Ruppert, D. and Wand, M. P. (1994). Multivariate locally weighted least squares regression.
T he Annals of Statistics, vol. 22, Number 3, 1346–1370.

Semiparametric regression with local divergence
Kanta NAITO
Department of Mathematics and informatics, Chiba University, Japan1

Setting

Let (Y1 , X 1 ), . . . , (Yn , X n ) ∼i.i.d. f (y, x) = p(y|x)q(x), where (Yi , X i ) ∈ R × Rd , f is the joint
density of (Y, X), p is the conditional density of Y given X = x, and q is the density of X.
Let t ∈ Rd be a target point at which we want to estimate the value of regression functionµ(t) = E[Y |X = t]. Our parametric model for µ is having the form()
e ,m(x, θ) = G−1 θ T x
e T = [1 xT ] ∈ Rd+1 is the vector of explanatory variables, θ = [θ0 θ1 · · · θd ]T ∈ Θ ⊂ Rd+1
where xis the parameter vector, and G is the link function.

2

The Bregman divergence

Now we fix a strictly convex function U . Then the discrepancy between µ(·) and its parametricmodel m(·, θ) = mθ (·) can be measured by the Bregman divergence defined as
DU ∗ (u(mθ ), u(µ))∫[U ∗ (u(m(x, θ))) − U ∗ (u(µ(x))) − µ(x){u(m(x, θ)) − u(µ(x))}] q(x)dx=d∫R=[U ∗ (u(m(x, θ))) − y · u(m(x, θ))] f (y, x)dydx
d
R×R
∫
+
[−U ∗ (u(µ(x))) + µ(x)u(µ(x))] q(x)dx,

(1)
(2)

Rd

where u = U ′ : the derivative of U , U ∗ is the convex conjugate of U : U ∗ (s) = supz∈R {zs−U (z)},and we have used a fact that (U ∗ )′ = u−1 = (U ′ )−1 in (1).
The usual parametric regression can be carried out by using a certain estimator θ̂ of thetrue value of θ. Necessary tools for this estimation scheme areρ(y, x, θ) = U ∗ (u(m(x, θ))) − y · u(m(x, θ)),u′ (m(x, θ))∂e,ρ(y, x, θ) = {m(x, θ) − y} ′ψ(y, x, θ) =x∂θG (m(x, θ))by which the estimator θ̂ and the true value θ ∗ of θ can be defined as∫θ ∗ = arg minρ(y, x, θ)dF (y, x),θ∈Θ R×Rd∫
θ̂ = arg minρ(y, x, θ)dFn (y, x),
θ∈Θ

R×Rd

1

(3)
(4)

(5)
(6)

where Fn is the empirical distribution function based on (Y1 , X 1 ), . . . , (Yn , X n ), F is the cumulative distribution function with its density f and 0d+1 is the zero vector in Rd+1 .
Note that θ ∗ in (5) is the minimizer of (2), and the minimizer of the empirical version of (2)is nothing other than θ̂ in (6). The regression function estimator can be obtained by plugging
θ̂ into θ in m(·, θ):µ̂G (x) = m(x, θ̂).

3

(7)

The local Bregman divergence

We organize a scheme of estimation of θ depending on t locally. Necessary functions are listed
as follows:()
x−tρ(t, y, x, θ) = K
ρ(y, x, θ),(8)h∂ψ(t, y, x, θ) =ρ(t, y, x, θ).
(9)∂θ
Here K(z) is a smooth unimodal integrable function symmetric around z = 0 satisfying K(0) =1, and h > 0 is the scalar bandwidth which controlls the degree of localization.
We note that (8) and (9) are localized version of (3) and (4) respectively, with the use of
the kernel K. Using these functions, we define the true parameter θ ∗ (t) and its estimator θ̂(t)as follows:∫θ ∗ (t) = arg minρ(t, y, x, θ)dF (y, x),θ∈Θ R×Rd∫
θ̂(t) = arg minρ(t, y, x, θ)dFn (y, x).
θ∈Θ

R×Rd

This local estimator θ̂(t) of θ ∗ (t) also suggests us to make a regression estimator defined asµ̂L (x) = m(x, θ̂(x)),(10)

which we call the local estimator of µ(x), because the involved estimator of parameter is determined locally. Since θ̂(x) can vary depending on x, µ̂L would be expected to be more flexible
than µ̂G . On the other hand, we call the estimator µ̂G (x) in (7) the global estimator of µ(x).

4Problem

Our interest is to investigate whether the local estimator µ̂L in (10) performs better than the
global estimator µ̂G in (7). By its construction and the shape of K, it is seen intuitively that||θ̂(t) − θ̂|| converges to zero in some sense, as n → ∞ and h → ∞, hence µ̂L ≈ µ̂G . Even insuch a situation, we think it is still not obvious that the global estimator µ̂G certainly defeats
the local estimator µ̂L .
In the symposium, several theoretical results which claim superiority of µ̂L over µ̂G werereported.

2

市場創造時点の特定に向けた文献調査：KIP を前提として 1
加藤淳一（久留米大学 商学部）
1．目的
本論文の目的は、先行研究のレビューにより、市場創造時点の特定方法の方向性を明らかにすることである。この目的に向けて、次のように議論を進める。2 章で、議論に先立ち、市場創造の意味を共有する。3 章で、2 章で確認した市場創造について、これまでの研究がどのように扱ってきているのかを吟味する。4 章で、本論文の結論をまとめる。
2．市場創造とは
本論文は、市場創造についての事前の理解を必要とする。そこで、以下の議論に先立って、市場創造について共有しておく。まず、「市場」創造を扱うのであるから、本論文の分析単位は市場となる。消費者個々人ではなく、それらを集計した市場に焦点を当てる。
次に、市場「創造」を扱う。ここで創造とは、評価軸（統計的には主成分軸）の移り変わりである。この「創造」について、既存のマーケティングとの対比を通じて述べる。マーケ
ティングでは、新しい市場（製品）の創造をその一生に例えて、製品ライフサイクルという時系列のグラフによりとらえる。この製品ライフサイクルにおいて、市場創造は売上の水準の上昇（縦軸の高さ）により捉えられている。
このようにマーケティングでは、これまで市場創造を売り上げの急速な上昇としている。
これに対して、本論文の定義は、評価（主成分）軸の継時的な移り変わりである。観測値の水準（縦軸で測ったときの高さ）ではなく、縦軸あるいは横軸という「軸」である。以上で、市場創造とは何かを説明した
3．文献レビュー
時系列分析の分野において、構造変化(structure change)と呼ばれる研究が蓄積されている。
それらは、とりわけ単位根過程との関連で例えば、黒住(2007)、市橋・金子(2009)、Byrne andPerman (2006)において議論されている。これら先行研究は、大まかにまとめれば、ある時点での傾きや切片の変化として構造変化をとらえようとしている。よって、これらの研究で想定していた構造変化は、本論文の問題意識（市場創造）とは異なっている。
そこで、より本論文の問題意識（市場創造）と近い研究を探索する。すると、近年になり
状態空間モデルが研究されている。状態空間モデルは、例えば、近年マーケティング分野を1本論文の成果の一部は、平成 31 年（令和元年）度久留米大学の科研費申請・採択に係る個人研究費支援「観光目的地への消費者ニーズの経時的な変化点の特定に向けた文献調査」として研究費の支援を受けました。ここに記して感謝を表します。
1

含めて次のような様々な応用研究がなされている。山口・土屋・樋口(2004)、佐藤・北川（2002）、
伊東（2016）、本橋・磯崎・長尾・樋口（2012）、青柳・佐藤（2015）、奥野・中村（2016）、
佐藤・樋口（2008）、本橋・樋口（2013）などである。
ここで状態方程式が、（評価）軸の変化に最も関連しそうである。だが、状態方程式は、状態の数値の変化としてモデル構築しなければならない。このように理解すると、本論文の市場創造とは異なる。そこで、時系列データへのもう 1 つのアプローチとして回帰分析、とりわけ SCT(Stepwise Chow Test)に注目する。
さて、時系列データの全期間の端から Chow test を繰り返す。すると、説明変数で捉えられていない時点が分かれば、その時点が評価軸の移り変わった時点と考えられる。このようにして考えれば、SCT により市場創造時点の特定できそうである。SCT は市場創造時点の特定方法の 1 つの選択肢として使用できるのではないかと考えられる。
以上で、時系列データを巡る既存研究を吟味し、現在のところ問題は残されているものの、評価軸の変化を捉える 1 つの方法として SCT の使用を提案した。最後に、本論文の目的に戻り、結論を整理する。
4．結論
本論文の目的は、先行研究のレビューにより、市場創造時点の特定方法の方向性を明らかにすることであった。本論文の議論を通じて、これまでの時系列分析の構造変化の研究や、
状態空間モデルの状態方程式は、どちらも軸を固定してレベルの変化を議論しており、本論文が想定した市場創造とは異なっていた。
そのような中でも、SCT の変化点は、別の軸を必要としている時点と解釈すれば、市場創
造時点の特定にも結びつくと考えられた。まだ不十分なところも残されているものの、本論文のレビューを通じて市場創造時点の特定方法に一定の方向性を示せた。今後、更にこの解釈の適切さを巡って考察を深めたい。
参考文献
青柳憲治・佐藤忠彦(2015)、「3 階層多変量状態空間モデリングによる動的市場反応形成メカニズムの解明」、『日本オペレーションズ・リサーチ学会和文論文誌』、Vol. 58、pp. 70-100。
市橋勝・金子慎治(2009)、「長期データ系列における構造変化の検定法比較分析」、『Discussion Paper IDEC, Hiroshima
University』、pp.1-41。
伊東宏樹(2016)、「状態空間モデルの実行方法と実行環境の比較」、『日本生態学会誌』、Vol. 66、p. 361-374。
Joseph P. Byrne & Roger Perman,(2006). "Unit Roots and Structural Breaks: A Survey of the Literature," Working Papers, pp.1-23.
黒住英司(2007)、「経済時系列分析と単位根検定：これまでの発展と今後の展望」、『Discussion Paper Series』、No.228、pp.
1-26。
本橋永至・樋口知之(2013)、「市場構造の変化を考慮したブランド選択モデルによる購買履歴データの解析」、『マーケティング・サイエンス』、Vol.21、N0.1、pp.37-59。
本橋永至・磯崎直樹・長尾大道・樋口知之(2012)、「状態空間モデルによるインターネット広告のクリック率予測」、『オ
ペレーションズ・リサーチ』、Vol. 57、No. 10、pp. 574-583。
奥野拓也・中村和幸(2016)、「個人別セールスプロモーション効果の推定」、『情報処理学会論文誌数理モデル化と応用（TOM）』、9 巻、3 号、pp.61 - 74。
佐藤忠彦・樋口知之(2008)、「動的個人モデルによる消費者来店行動の解析」、『日本統計学会誌』、第 38 巻、第 1 号、pp.
1-19。
佐藤忠彦・北川源四郎(2002)、「マルコフ切り換え回帰モデルによる POS データの解析」、『行動計量学』、第 29 巻、第 2号、pp.233-245。
山口類・土屋映子・樋口知之(2004)、「状態空間モデルを用いた飲食店売上の要因分解」、『オペレーションズ・リサーチ』、
Vol. 49、No. 5、pp. 316-324。

2

Estimation, Diagnostics, and Extensions of Nonparametric
Hawkes Processes with Kernel Functions
Jiancang Zhuang, Institute of Statistical Mathematics, Tachikawa, Tokyo 190-8562, JapanThe Hawkes self-exciting model has become one of the most popular point-process models inmany research areas in the natural and social sciences because of its capacity for investigating
the clustering effect and positive interactions among individual events/particles. It describes thestochastic excitations among a series of events that occur in a continuous time domain or in a
spatiotemporal domain. A spatiotemporal marked Hawkes process, supposing its realization N ={ti : i ∈ Z} with Z being the set of all integers, has a conditional intensity in the formZλ(t, x, m) = µ(t, x, m) +g(t − s, x − u) f (m | m0 ) dN (ds × du × dm0 )(1)(−∞,t)×Rd ×M

where x and m denote the location in the space of Rd and the mark in the space of M, respectively,
and f (m | m0 ) gives the conditional p.d.f. for the magnitudes of direct offspring from an event ofmagnitude m0 .
For such a Hawkes process, the probability that an event, say j, is a background event, is
ϕj = Pr{Event j is a background event} =µ(tj , xj , mi )λ(tj , xj , mi )(2)

and the probability that event j is triggered by another event i is
ρij = Pr{Event j is triggered by i} =g(tj − ti , xj − xi )f (mj | mi ).
λ(tj , xj , mj )(3)

Making use of the above probabilities, we can construct an Expectation-Maximization algorithm
for parametric Hawkes models and an iterative stochastic reconstruction algorithm for nonparametric or semi-parametric Hakes models. When both the background rate and triggering term arenonparametric, for example, to estimate a Hawkes with conditional intensityXλ(t) = νµ(t) + A
g(t − ti )(4)i:ti <tas the conditional intensity function, and the estimates for µ and g becomeXXµ̂(t) ∝
ϕi Zh1 (t − ti )and ĝ(t) ∝ρij Zh2 (tj − tj − t).
i

(5)

i,j

RT
R∞
with restrictions 0 µ̂(t) dt/T = 1 and 0 ĝ(t) dt = 1. In the above equations, the parameters νand A are the relaxation coefficients and are estimated by the MLE or EM method.
Algorithm 1.
1. Stochastic declustering (Expectation). Calculate the background probability and triggering probabilities.
2. Reconstruction (Maximization I). Estimate the nonparametric function in the modelusing nonparametric methods such as kernel functions.
3. Parametrization (Maximization II). Use the MLE method or EM algorithm to estimate
the parameters in the parametric functions and the relaxation coefficients for the nonparametric functions.
The above method is applied to the data analysis in the studies of earthquakes and crimes. In
the crime data, periodicity of the background rate is considered.

非定常クラスター過程のパラメータ推定法：Palm 尤度法の拡張
島谷健一郎（統計数理研究所）
空間点過程のパラメータ推定は、定常モデルから非定常へシフトしている。非定常モ
デルは定常モデルから定常性の仮定を外せば作られるが、どの定常性を外すかにより、いくつかの非定常化が考えられる。植物生態学では、定常ポアソンで分布する親の周りに散布された娘集団（親は既に死亡して位置は不明）という定常クラスター点過程(homogeneous Neyman-Scott process) がよく使われる。ここに非定常性を入れる場合、
1. 定常ポアソンの親集団を非定常ポアソンにする。
2. 定常な親から散布された娘に非定常な生残率を与える。
3. 定常な親集団から、散布距離（散布カーネルのパラメータ、正規分布なら標準偏差）を非定常にする。
4. 1-3 の非定常性の 2 つ以上を同時に入れる。
などの選択肢がある。つまり、実データからのパラメータ推定では、同時に、どの非定常性を加味したモデルを考えるべきかという、モデル選択の問題に直面する。
ここでは、定常な場合に田中・尾形によって提唱された Palm 尤度法を、上記 1 と 2の非定常に拡張し、それでモデル選択とパラメータ推定を行う手法を提唱する。Palm
最尤法に加え、Palm 尤度を用いるベイズ推定も提唱する。そして、簡単なシミュレー
ションで方法論の有効性と限界を検証した結果を報告する。人工点配置データとして、親集団の配置も娘集団に働く生残率も非定常とする full model に加え、いずれかを定常とするモデルでも点配置データを生成した。それらの人工データに、full model に加え、親集団または娘生残を定常（一定）とするモデルの Palm 尤度によるベイズ及び最尤推定も行った。
Palm 尤度を用いたベイズ推定では、MCMC は 1000 回程度できれいに単峰型に収束を示す場合と、振動する場合があり、親密度と娘生残の両方を非定常とするモデルでは、後者が多かった（図 1）。これは、現在の点配置データだけから、両者を分離して推定することが、一般には困難であることを示す。
一方、いずれかを定常として生成した点配置では、MCMC は明白な単峰型を示し、最尤推定値も求まった。さらに、点配置を生成したモデルの最大対数 Palm 尤度は、fullmodel のものと 1.0 を下回る差しか持たず、赤池情報量規準 AIC で真の生成モデルが選ばれる結果となった（表 1）。
以上の結果をまとめると、一般に、親密度と娘の生残の同時推定は困難で、パラメー
タ推定値は振動したり、両者の中間になったりする。でも、うまく推定できる場合もある。また、単純化されたモデルで生成された点配置なら、AIC で真のモデルが選ばれ、パラメータ推定も悪くない。
なお、点配置データでは、サンプルサイズは、点の数でなく（これはデータの一部）、
プロット（window)の大きさが対応する（漸近正規性も window の拡大で示されている）。
点の数が多ければ推定精度が上がるわけではない。Neyman-Scott 過程では、親数が肝心で、実際、非定常ポアソン過程の推定では、200 点くらい欲しい。しかし、親は既にいないので、今のデータからは不明で、推定して初めて目星がつく。人工データでは 親
数 > 200 で生成したからうまく推定できた。実データではやるまでわからない。
従って、いきなり複雑なモデルを膨大な MCMC で推定するより、本手法のような、簡略された推定法で何が推定でき何はそもそも無理そうか、試すのが賢明と言えよう。
/2

1

図1
1
2
3
0 4
-1 05
-2 6
-3 7
8
-4
9
10
2.011
1.5
12
1.0
0.513
0.014
-0.5150
-1.0
16
-1.5
-2.017-2.51819
5.0 20
4.5 21
4.0 223.5
23
3.0
2.5 24
2.0 251.5
26
1.0
0.5 27
0.0 280
2930
4.5 31
4.0 32
3.5 33
3.0 342.5
35
2.0
1.5 36
1.0 37
0.5 38
0.0 39
400
41
2
1

0
0.00000.05590.1190200
0.1190400
0.12400.12400.13300.13300.13300.17170.17330.17330.1733200
0.2126400
0.21260.21260.21260.27530.27530.27530.29820.29820.27900.27900.29200.29200.3113
2000.3057
400
0.33690.34530.35780.35780.39830.45640.44860.4486
0.52250.5225200
0.5225400
0.5225-1.1555 3.2740-1.1555 3.2740-1.2003 3.2757-1.2356 3.2758
600
1000-1.2356800 3.2758-1.1970 3.2911-1.1970 3.2911-1.1714 3.3281-1.1714 3.3281-1.1714 3.3281-1.2779 3.3401-1.2647 3.3647-1.2647 3.3647-1.2647 3.3647
600
8003.3732
1000-1.2336-1.2336 3.3732-1.2336 3.3732-1.2336 3.3732-1.2419 3.3931-1.2419 3.3931-1.2419 3.3931-1.2945 3.3310-1.2945 3.3310-1.2707 3.3403-1.2707 3.3403-1.3054 3.3636-1.3054 3.3636-1.3768 3.3430
600
1000-1.3944800 3.3366-1.4312 3.3637-1.4709 3.3932-1.4794 3.3706-1.4794 3.3706-1.4486 3.3584-1.4287 3.3115-1.4219 3.3547-1.4219 3.3547-1.5366 3.3878-1.5366 3.3878
600
1000-1.5366 8003.3878-1.5366 3.3878

00.00000.09040.14581200
1400
0.14580.14640.14640.14560.14560.14560.16070.17270.17270.17271200
1400
0.19110.19110.19110.19110.18130.18130.18130.20250.20250.21230.21230.20470.20470.24851200
1400
0.30290.33610.37720.39190.39190.38760.4556
0.55350.55350.5314
0.53141200
1400
0.5314
0.53141
0 -1.4162 4.0754
2 0.0000 -1.4162 4.0754
3 0.0548 -1.3907 4.1525
4 0.0860 -1.3730 4.1971
1
5 0.1250 -1.3103 4.2345
6 0.1643 -1.3484 4.2276
7 0.1643 -1.3484 4.2276
0
8 0.1643 -1.3484 4.2276
0
200
400
600
800
9 0.1643 -1.3484 4.2276
10 0.2271 -1.4002 4.2399
0.0
11 0.2312 -1.4045 4.2081
0 12 0.2884
200
400-1.3542
600 4.2568
800-0.5
13 0.2884 -1.3542 4.2568-1.0
14 0.2884 -1.3542 4.2568
15 0.2884 -1.3542 4.2568-1.5
16 0.3241 -1.3357 4.1898-2.0
17 0.3241 -1.3357 4.1898
18 0.4296 -1.3257 4.1622-2.5
19 0.4543 -1.2908 4.2317
20 0.4543 -1.2908 4.2317
2.5
21 0.4543 -1.2908 4.2317
2.0
22 0.4814 -1.2841 4.3292
23 0.5566 -1.2999 4.2438
1.5
24 0.5962 -1.3976 4.1852
1.0
25 0.6548 -1.4805 4.1813
26 0.6548 -1.4805 4.1813
0.5
27 0.7956 -1.3742 4.1818
0.0
28 2000.7956 400-1.3742
4.1818
0
600
800
-0.5
29 0.7956 -1.3742 4.1818
30 0.7956 -1.3742 4.1818
31 0.7899 -1.3739 4.1359
5.0
32 0.7971 -1.3366 4.1885
4.0
33 0.7971 -1.3366 4.1885
34 0.7971 -1.3366 4.1885
3.0
35 0.7971 -1.3366 4.1885
2.0
36 0.8252 -1.3630 4.1880
1.0
37 0.8562 -1.4180 4.1202
38 0.8562 -1.4180 4.1202
0.0
0.8562400 -1.4180
0 39 200
600 4.1202
800
40 0.8542 -1.4597 4.1071
2

親密度傾き

1600

1800

2000

2200

2400

2600

親密度切片

1600

1800

2000

2200

2400

生残傾き

大きく振動
1600

2600

1800

2000

2200

2400

2600

散布

最尤推定値に意味は薄い
1600

1800

2000

2200

0 -89853.72
親密度傾き#### -89853.72#### -89819.30#### -89757.54#### -89663.11#### -89594.56#### -89594.56#### -89594.56
1000
1200
1400
1600
1800
2000
2200
2400
2600
#### -89594.56#### -89522.97#### -89475.46
親密度切片1000 -89359.21
1200
1400
1600
1800
2000
2200
2400
2600
####
#### -89359.21#### -89359.21#### -89359.21#### -89344.42#### -89344.42#### -89250.86#### -89221.07#### -89221.07
生残傾き#### -89221.07#### -89194.34#### -89066.58#### -89061.63#### -89049.42#### -89049.42#### -88926.09####-88926.09
1000
1200
1400
1600
1800
2000
2200
2400
2600
#### -88926.09#### -88926.09#### -88871.12
散布#### -88842.19#### -88842.19#### -88842.19#### -88842.19#### -88787.01#### -88732.86#### -88732.86#### 1200-88732.86
1000
1400 1600 1800 2000 2200 2400 2600#### -88724.47

1点の周辺を振動

最尤推定値が求まる

2400

2600

表1
親と生残が同じ減少傾向 θ 1 = φ 1
データ生成：親も生残も非定常、x に沿って減少
親密度 親密度 種子散 生残率
親密度 親密度 種子散 生残率
傾き切片
布
傾き 最大対数 傾き
切片
布
傾き 最大対数尤
尤度
度
データ 3
θ1
Θ0
φ
θ
Θ
φ
σ
σ
1
1
0
1
1
-1.5
3
2
2
-2.2
3
2
真値
Full model
1.30
-1.51
2.91
1.97 -88148.8
2.55
-2.99
3.29
2.52 -142545.8
θ1 = φ1
1.66
-1.85
2.94
1.66 -88149.3
2.54
-2.97
3.28
2.54 -142545.8
親：非定常
5.45
-5.29
3.08
0 -142640.2
生残：定常
親：定常 生0 -0.33
2.91
3.10 -88155.7
0
-0.57
3.26
4.58 -142565.4
残：非定常
データ 2
親は非定常、生残は定常
1
-1.5
3
2
2
-2.2
3
0
真値
Full model
-1.53
0.64
3.08
3.88 -69091.1
1.94
-2.80
3.95
-0.08
-84661.5
θ1 = φ1
1.39
-1.99
3.12
1.39 -69135.6
0.72
-1.73
3.76
0.72
-84678.2
親：非定常
1.81
-2.69
3.94
0 -84661.9
生残：定常
親：定常 生
0
-1.11
3.76
1.23
-84700.3
残：非定常
親は x に沿って減少、生残は x に沿って増加
親は定常、生残は非定常
2
-2.2
3
-2
0
-0.8
3
2
真値
Full model
0.54
-1.61
3.20 -0.93 -46827.9
-0.31
-0.81
3.63
1.97
-68471.7
θ1 = φ1-0.30-1.13
3.22 -0.30 -46841.4
0.95
-1.76
3.52
0.95
-68489.1
親：非定常
-0.59
-0.91
3.19
0 -46856.0
2.00
-2.34
3.26
0 -68540.6
生残：定常
親：定常 生0 -1.31
3.23 -0.55 -46833.2
0
-1.06
3.63
1.73
-68472.5
残：非定常

参考文献
Tanaka, U., Ogata, Y. & Stoyan, D. 2008. Parameter estimation and model selection for
Neyman-Scott point processes. Biometrical Journal, 50, 43–57.
/22