シンポジウム「統計科学の現代的課題」報告書
本シンポジウムは科学研究費・基盤研究（A)「大規模複雑データの理論と方法論の総合的研究」（研究代表者：青嶋 誠（筑波大学）, 課題番号：15H01678）の助成により開催した。
■日時： 2017年1月27日～1月29日■場所： 金沢大学サテライトプラザ■開催責任者： 星野 伸明 （金沢大学）【 プログラム・目次 】
講演者氏名 所属

日時
セッション
1月27日 １<座長：蛭
川雅之＞

13:00 西埜晴久

広島大学大学院社会科
学研究科

13:40 松井宗也

南山大学経営学部

14:20 生亀清貴

東京理科大学理工学部

演題頁

経済不平等度を予測するためのモデルについて1

組み合わせ論と無限分解可能分布のひとつの接点多変量確率密度関数の二重対称性について3
5

15:00 休憩
セッション
２＜座長：
加藤昇吾
＞

15:15 森山卓

九州大学大学院数理学
府

15:55 蛭川雅之

摂南大学経済学部

16:35 鶴田靖人

金沢大学人間社会環境研究科

二標本ノンパラメトリック検定の連続化と局所漸近検出力
Nonparametric Estimation andTesting on Discontinuity of
Positive Supported Densities: A
Kernel Truncation Approach
円周上のカーネル密度推定量とそのバンド幅選択法の漸近的
性質

7

9

11

17:15 休憩
セッション
３＜座長：
生亀清貴
＞

17:30 加藤昇吾

統計数理研究所

角度データのための統計モデ
ル

13

18:30 終了
セッション
４＜座長：
1月28日
西埜晴久
＞

10:00 国友直人

明治大学政治経済学部

10:40 本山 要

上智大学大学院理工学
研究科数学領域

11:20 塚原英敦

成城大学経済学部

The Simultaneous Multivariate
Hawkes-type Point Processes
and their application to
Financial Markets
ウェーブレット分散と参照形式を利用した商品先物の暴落予
測On Backtesting Risk
Measurement Models

15

17
19

12:00 昼休み
セッション
５＜座長：

13:30

加葉田大志 大阪市立大学大学院医朗
学研究科

交互作用を含む観察データにおける傾向スコア解析手法による結果の差異

21

国友直人
＞

14:10 高橋佳苗
14:50 石原拓磨

大阪大学医学部附属病 診断研究における検査性能の院
比較
大阪大学大学院 医学系 複数の二値評価変数をもつ臨
研究科
床試験における検定法

2325

15:30 休憩
セッション
６＜座長：
塚原英敦
＞

15:45 新村秀一

成蹊大学 経済学部

16:25 山本けい子

函館工業高等専門学校
一般理数系

17:05 山縣 一慶

上智大学大学院理工学
研究科数学領域

17:45 竹下 佳宏

上智大学大学院理工学
研究科数学領域

R.A.Fisher 以後の判別分析の新理論と遺伝子解析の新手法
産業連関表に基づく都道府県クラスターと産業構造推移の可視化
テキストマイニングにもとづくレ
ビューのスコアリングを用いた映画の統計的分類
規模の異なる変量群をもつロジスティック回帰モデルの係数2
段階推定

27

29

31
33

18:25 終了
19:00 懇親会
セッション
1月29日 ７＜座長：松井宗也＞

10:00 伊藤伸介

中央大学経済学部

10:40 紙屋英彦

大阪経済大学経済学部

11:20 吉田知行

北星学園大学経済学部

12:20 終了

国勢調査の匿名化ミクロデータの作成可能性について―地域区分に着目して―Estimation of the shape ofdensity level sets of starshaped distributions
正方分割表の一致率検定のための代数的方法3537
39

経済不平等度を予測するためのモデルについて
広島大学大学院社会科学研究科

西埜晴久∗

日本経済における経済の不平等度の研究は，2000 年代初めに橘木俊詔（1998) および大竹文雄（2005）において行われてきて，大きな関心を集めてきた．そこで，不平等度の時
系列的な変動についても興味深い問題対象であると考えている．
本報告では，Nishino, H. and K. Kakamu (2011) で用いられているように対数正規分布を仮定して分位データを用いてモデル化する．なお，日本の家計調査のデータに対して対数正規分布は適合的である. そして，本報告ではそのモデルに GARCH 型の時系列モデルを組み込むことで，不平等度を予測するための時系列モデルを組み立てることを試みるものである．
不平等度を取り扱う日本のデータとして，本報告では総務省統計局の家計調査を用いる．
対数正規分布は (µ, σ) の二つのパラメータをもつ分布である．一方，対数正規分布のジニ
係数は
(
)
σ
G = 2Φ √− 1,(1)
2とあらわされる．
（たとえば，Kleiber and Kotz (2003) を参照のこと）. つまり ジニ係数はσ だけの単調増加関数であり，σ の大きさが不平等度の大きさを表すことになる．また，2つのパラメータを持つ分布のなかでは対数正規分布が家計調査の 5 分位のデータには AIC
で当てはまりを比較する限り，もっとも適合であることが分かる．
標本サイズを n, 部分標本サイズを k としたときに，分位データは部分順序統計量として {X(n1 ) ≤ X(n2 ) ≤, . . . , ≤ X(nk )} (1 ≤ n1 < n2 < . . . < nk ≤ n) と書ける．家計調査の 5 分位データでは，n = 10000, k = 4, (n1 , n2 , n3 , n4 ) = (2000, 4000, 6000, 8000) とな
る．また，(x1 , x2 , . . . , xk ) を {X(n1 ), X(n2 ), . . . , X(nk )} の観測値とする．また，対数変換値 y = (y1 , y2 , . . . , yk )′ = (log x1 , log x2 , . . . , log xk )′ とする．なお，対数正規分布を対数変換すると正規分布になる．
部分順序統計量の同時密度は，David and Nagaraja (2003) から，f (y|θ) ={ k}∏ (F (yi ) − F (yi−1 ))ni −ni−1 −1F (yi )n1 −1
n!
f (y1 )f (yi )
(n1 − 1)!
(ni − ni−1 − 1)!
i=2(1 − F (xk ))n−nk
×(n − nk )!

(2)と書ける．ただし，F (·) は累積分布関数であり，f (·) は確率密度関数である．
グループ・データを (x1 , . . . , xk ) とし，その対数変換を (y1 , y2 , . . . , yk ) = (log x1 , log x2 , . . . , log xk )とする．対数正規分布を仮定しているため, Y = log X ∼ N (µ, σ 2 ) となる．

∗

〒 739-8525 東広島市鏡山 1-2-1. Email: hnishino@hiroshima-u.ac.jp- 1 -yt = (y1,t , y2,t , . . . , yk,t )′ と表示すると，GARCH(1,1) model は
yt ∼ f (yt |µt , σt ),

(3)

2
2
σt2 = α0 + α1 σt−1
+ β1 σ̃t−1, (GARCH)(4)

σt = exp {α0 + α1 log(σt−1 ) + β1 log(σ̃t−1 )} ,(EGARCH)(5)と書ける．ただし，{σ̃t } は yt = (y1,t , y2,t , . . . , yk,t )′ からの σt の推定値である．
また，推定は以下の疑似最尤法を用いる．つまり，logL =

T∑

log f (yt |µt , σt )

(6)

t=2

を最大化するパラメータを推定値とする．
次に一期先予測は，例えば，EGARCH(1,1) に対しては，{}
σ̂t = exp α̂0 + α̂1 log(σt−1 ) + β̂1 log(σ̃t−1 ) ,(7)とする．ただし，{σ̃t } は yt = (y1,t , y2,t , . . . , yk,t )′ からの σt の推定値であり，α̂0 , α̂1 , β̂1
パラメータの推定値である．さらに，予測の観点からモデルを比較すると平均予測誤差の平方根 (root MSE) を
v
u
T
u 1 ∑
t
{σ̂t − σ̃t }2
(8)
T −1
t=2

用いることにする．なお，root MSE によって家計調査のデータにモデルの当てはまりを判断すれば，共に α1 = 0 となる，GARCH(0,1) 次いで，EGARCH(0,1) が良いことが分かった．しかし，さらに当てはまりの良いモデルを開発することを今後検討している．

References[1] David, H. A. and Nagaraja, H. N. (2003) Order Statistics, 3rd ed. Wiley, New York.
[2] Kleiber C. and Kotz, S. (2003) Statistical Size Distributions in Economics and Actuarial
Sciences, Wiley, New York.
[3] Nishino, H. and Kakamu, K. (2011) “Grouped Data Estimation and Testing of Gini Coefficients using Lognormal Distributions,” Sankhyā B 73, 193–210.
[4] 大竹文雄 (2005) 『日本の不平等』, 日本経済新聞社.
[5] 橘木俊昭 (1998) 『日本の経済格差』, 岩波新書.

- 2 -組み合わせ論と無限分解可能分布のひとつの接点
南山大学

松井 宗也

研究内容
組み合わせ論では、その解析に母関数 (Generating function) がよく用いられ、諸処の問題に有効となることも多い (Wilf [5])。この母関数は、確率・統計において非負整数値の確率変数の確率母関数に対応する。本報告では、対称群の cycle index
polynomials が非負整数値の無限分解可能分布（それは複合ポアソン分布と同値）の確率関数に対応することに注目する。そして無限分解可能分布に関する既知の結果をもとに、組み合わせ論の新しい定理を得る。さらに、その定理の組み合せ論的な証明や、その定理から敷衍される結果も述べる。
より具体的な内容を述べる。非負整数列 (ak )k≥0 、(bk )k≥0 と (ck )k≥1 は以下の関係を満たすものとする。
∞
∑
k=0

k

ak u =

∞
∑
bk uk
k=0

k!

= exp

∞
(∑
c uj )
j

j=1

j

.

(1)

このとき数列 (ak )k≥0 と (bk )k≥0 の性質を (ck )k≥1 のそれから導くという問題を考える。Bender and Canfield [1] は、数列 (ck )k≥1 と c0 = 1 からなる数列が log-concave
であると仮定し、そのうえで数列 (ak ) が log-concave と almost log-convex であることを示した。厳密には不等式：

k+1
ak−1 ak+1 ,kkbk−1 bk+1 ≥ b2k ≥bk−1 bk+1
k+1

ak−1 ak+1 ≤ a2k ≤が成り立つことを示した 1 。それは組み合わせ論的な証明による。研究の主結果は、対応する関係を log-convex のそれに置き換えたものへ拡張したことである。つまり、
数列 (ck )k≥1 と c0 = 1 の log-convexity から数列 (ak ) のそれを導いた。詳しくは以下の同値な不等式を示した。

ak−1 ak+1 ≥ a2k ,

b2k ≤

k
bk−1 bk+1 .
k+1

証明は、式 (1) の確率的解釈と非負整数の値を取る複合ポアソン分布の性質による。
この証明方法は、最初の log-concavity の関係の別証にもなる。
次に、拡張された結果 (log-convex に関する) の組み合わせ論的な証明も与える。
組み合わせ論的な証明では log-concave/convex のみならず、更に強い結果を得られ1ほとんど log-convex であることは第１式の左側の不等式を指す。なお、第１式の右側の不等式と第２式の左側の不等式は同値である。

- 3 -る。以下に述べるのは Bender and Canfield [1] の組み合わせ論による結果とその拡張の両方である。まず非負数列 (ck ) から構成される集合を以下に定義する。

X = {c1 , c2 , . . .},Y = X ∪ {cj ck − cj−1 ck+1 : 0 < j ≤ k},
Z = X ∪ {cj−1 ck+1 − cj ck : 0 < j ≤ k}.
すると、式 (1) を満たす数列 (bk ) の log-concave/convex form を (ck ) のそれの多項式で以下のように表せる。

(n + 1)bm bn − mbm−1 bn+1 ∈ N[Y ],mbm−1 bn+1 − (n + 1)bm bn ∈ N[Z ]/N[X ] for 1 ≤ m ≤ n,
ここで N[Y ] は、集合 Y の要素からなり、かつ非負整数の係数を持つ多項式である。N[Z ] と N[X ] も同様に定義される。最初の式が Bender and Canfield [1] による log-concave に関する結果であり、次の式がそれを log-convex へ拡張したものである。
さらに、式 (1) を満たす数列 (ak ) や (bk ) の convolution や binomial convolution
によりできる数列の性質を与える。式 (1) の積により (ak ) の convolution や (bk ) のbinomial convolution が得られる。この操作により得られる数列は、これまでに確立された convolution に対する log-convex/concave の法則に必ずしも従わない。こ
れをみるために、離散確率変数の確率関数を数列とみなし、それにより具体例を構成した。
キーワード：Cycle index polynomials, compound Poisson, symmetric group, infinitely
divisible, generating function, log-convexity.

参考文献[1] E. A. Bender and E. R. Canfield (1996) Log-concavity and related properties of
the cycle index polynomials, J. Combin. Theory Ser. A 74, 57–70.
[2] B. G. Hansen (1988) On log-concave and log-convex infinitely divisible sequences anddensities. Ann. Probab. 16, 1832–1839.
[3] M. Matsui (2016) Log-convexity and the cycle index polynomials with relation tocompound Poisson distributions, ArXiv:1609.06875.
[4] F. W. Steutel and K. van Harn (2004) Infinite Divisibility of Probability Distributions on the Real Line. Marcel Dekker, New York.
[5] H. S. Wilf (1994) Generatingfunctionology. 2nd. ed. Academic Press, San Diego.

- 4 -多変量確率密度関数の二重対称性について
東京理科大学 理工学部

生亀清貴 ∗Tomizawa (1985) は二元分割表において二重対称モデル，準二重対称モデル，周辺二重対称モデルを提案し, 定理「二重対称モデルが成り立つための必要十分条件は，準二重対称モデルと周辺二重対称モデルの両方が成り立つことである」を与えた. Yamamoto
et al. (2012) はこの定理を多元分割表に拡張した．
本講演では多変量確率密度関数に対して二重対称性, 準二重対称性と周辺二重対称性を定義し, 二重対称性をもつ確率密度関数の分解を与える.
確率変数 (X1 , . . . , XT ) は確率密度関数 f (x1 , . . . , xT ) をもつとする. ただし
f (x1 , . . . , xT ) > 0 for (x1 , . . . , xT ) ∈ DT ,f (x1 , . . . , xT ) = 0 for (x1 , . . . , xT ) ∈/ DT ,
またDT = {(x1 , . . . , xT ) | ai < xi < bi ; i = 1, . . . , T },ただし ai = −∞ かつ bi = +∞, または ai , bi は共に有限とする. (c1 , . . . , cT ) を DT 内の点とする, ただし ai , bi が共に有限の場合は ci = (ai + bi )/2. また任意の i (i = 1, . . . , T )に対して, Xi = xi のとき x∗i = 2ci − xi とおく.
確率密度関数 f (x1 , . . . , xT ) の二重対称性を次のように定義する：(1, . . . , T ) の任意の並べ替え (π1 , . . . , πT ) に対して,f (x1 , . . . , xT ) = f (xπ1 , . . . , xπT )= f (x∗1 , . . . , x∗T ).
また任意の k (k = 1, . . . , T − 1) に対して, k 次周辺二重対称性を次のように定義する：

fXi1 ...Xik (xi1 , . . . , xik ) = fXi1 ...Xik (xπi1 , . . . , xπik )= fXj1 ...Xjk (xi1 , . . . , xik )= fXi1 ...Xik (x∗i1 , . . . , x∗ik )(1 ≤ i1 < · · · < ik ≤ T ; 1 ≤ j1 < · · · < jk ≤ T ),ただし (πi1 , . . . , πik ) は (i1 , . . . , ik ) の任意の並べ替え, fXi1 ...Xik は (Xi1 , . . . , Xik ) の周辺
確率密度関数とする.

2010 Mathematics Subject Classification: 62H05
キーワード：二重対称性, 準二重対称性, 周辺二重対称性，確率密度関数, 分解∗〒 278-8510 千葉県野田市山崎 2641
e-mail: iki@is.noda.tus.ac.jp- 5 -一般に, 確率密度関数はf (x1 , . . . , xT ) = α
[
×

∏

T
[∏

][ ∏ ∏]
αi1 (xi1 )αi1 i2 (xi1 , xi2 ) × · · ·

i1 =1

···

∏

1≤i1 <i2 ≤T]
αi1 ...iT −1 (xi1 , . . . , xiT −1 ) · α1...T (x1 , . . . , xT ),(1)

1≤i1 <···<iT −1 ≤Tのように表すことが可能である, ただし (x1 , . . . , xT ) ∈ DT ,{αi (ci ) = αi1 i2 (ci1 , xi2 ) = · · · = α1...T (x1 , . . . , xT −1 , cT ) = 1}.
このとき, 確率密度関数 f (x1 , . . . , xT ) が二重対称であるための必要十分条件は, 式 (1)に次の制約を課すことである：

αi1 ...im (xi1 , . . . , xim ) = αi1 ...im (xπi1 , . . . , xπim )= αj1 ...jm (xi1 , . . . , xim )= αi1 ...im (x∗i1 , . . . , x∗im )(m = 1, . . . , T ; 1 ≤ i1 < · · · < im ≤ T ; 1 ≤ j1 < · · · < jm ≤ T ),ただし (πi1 , . . . , πim ) は (i1 , . . . , im ) の任意の並べ替え. 任意の k (k = 1, . . . , T − 1) に対して, f (x1 , . . . , xT ) の k 次準二重対称性を次のように定義する：式 (1) に対して,

αi1 ...im (xi1 , . . . , xim ) = αi1 ...im (xπi1 , . . . , xπim )= αj1 ...jm (xi1 , . . . , xim )= αi1 ...im (x∗i1 , . . . , x∗im )(m = k + 1, . . . , T ; 1 ≤ i1 < · · · < im ≤ T ; 1 ≤ j1 < · · · < jm ≤ T ),ただし (πi1 , . . . , πim ) は (i1 , . . . , im ) の任意の並べ替え.
このとき, 次の定理を得る.
定理： 任意の固定された k に対して (k = 1, . . . , T − 1), 確率密度関数 f (x1 , . . . , xT ) が二重対称性を満たすための必要十分条件は, f (x1 , . . . , xT ) が k 次準二重対称性と k 次周辺二重対称性の両方を満たすことである.

参考文献[1] Tomizawa, S. (1985). Double symmetry model and its decomposition in a square contingency table. Journal of the Japan Statistical Society, 15, 17-23.
[2] Yamamoto, K., Takahashi, F. and Tomizawa, S. (2012). Double symmetry model and its
orthogonal decomposition for multi-way tables. SUT Journal of Mathematics, 48, 83-102.

- 6 -二標本ノンパラメトリック検定の連続化と局所漸近検出力(Smoothed two-sample nonparametric tests and theirlocal asymptotic powers)
Taku MORIYAMA(Graduate School of Mathematics, Kyushu University)∗1
Yoshihiko MAESONO (Faculty of Mathematics, Kyushu University)∗2

1. Introduction
Let X1 , X2 , · · · , Xm be independently and identically distributed random variables(i.i.d.) from distribution function F (x) and Y1 , Y2 , · · · , Yn be i.i.d. from F (x−θ) whereθ is unknown location parameter. We assume that m + n = N and 0 < limN m/N =limN λN = λ < 1. We consider ‘2-sample problem’ whose null hypothesis H0 : θ = 0
and H1 : θ > 0.
Let us define that ψ(x) = 1 (x ≥ 0), = 0 (x < 0) and Z1 = X1 , · · · , Zm =Xm , · · · , Zm+1 = Y1 , · · · , ZN = Yn . The median and Wilcoxon’s rank sum test statistic
are
n
∑
∑ ∑
M=
ψ(Yj − Z), W2 =ψ(Yj − Xi ).
j=1

1≤i≤m 1≤j≤nwhere Z denotes the sample median of {Z1 , Z2 , · · · , ZN }. For observed values x =(x1 , x2 , · · · , xm )T and y = (y1 , y2 , · · · , yn )T , we put s = s(x, y) and w2 = w2 (x, y) are
the realized value of M and W2 . If the p-values P0 (M ≥ m) or P0 (W2 ≥ w2 ) is smallenough, we reject the null hypothesis H0 .

2. The smoothed nonparametric tests
We propose the following smoothed median test statistic)(m∑
Z − Xi
∗
f
M=
K
h
i=1
where K ∗ is an integral of the kernel k ∗ (t) which satisfies∫ ∞k ∗ (t)dt = 1 and k ∗ (t) = 0 for t ≤ 00and h is a bandwidth which satisfies h → 0 and nh → ∞. Further,we propose thefollowing Wilcoxon’s rank sum test
)
(
m ∑
n
∑
Y
−
X
j
i
f2 =
W
K
h
i=1 j=1
where K is an integral of a symmetric kernel k. We can easily prove their asymptotic
normality and that their Pitman eﬃciency are same as the discretes :f] = eP [M † ],e P [Mf2 ] = eP [W2 ]e P [W

where M † is equivalent to M .
Keywords: nonparametric test, kernel estimator, significance probability, local asymptotic power.
∗1
e-mail: moritaku3542168@gmail.com∗2
e-mail: maesono@math.kyushu-u.ac.jp- 7 -Table 1: Comparison of significance probabilities of M and W2
sample size (m, n)
z0.9
z0.95
z0.975z0.99(10,10) (20,20) (30,30)
1.68
2.21
2.79
6.37
5.66
2.17
3.72
3.05
2.80
33.7
7.54
1.56

(10,20)
7.13
3.44
16.6
6.05

(20,10)
7.07
3.41
14.7
5.29

∗
(Um
, Un∗ )1
3.86
4.21
4.10
4.11

3. Significance probability
Table 1 shows the ratio of frequency of exact p-value of W2 getting smaller than M in
the tale area Ωα
{
}
m(x)
−
E
(M
)
w
(x)
−
E
(W
)
0
2
0
2
√
√
Ωα = x ∈ Rn
≥ v1−α , or
≥ v1−αV0 (M )
V0 (W2 )where v1−α is a (1 − α)th quantile of the standard normal distribution N (0, 1), and
E0 (·) and V0 (·) stand for an expectation and a variance under H0 , respectively. We
count samples that an exact p-value of the test is smaller than the other in the tailarea Ωα , and calculate the ratio of the frequency.
Because the values in Table 1 are larger than 1, we find that W2 tends to havesmaller p-value than M . They can let one to use W2 if one wants the small value. This
comes from that the possible p-value of M is more sparse than W2 like the sign test
and Wilcoxon’s signed rank test as is discussed by Maesono et al.(2016). On the otherf and Wf2 in table 2 are close to 1.
hand, we can find the ratio values of M
f and W
f2
Table 2: Relation of significance probabilities of M
sample size (m, n)
z0.9
z0.95
z0.975z0.99(10,10) (20,20) (30,30)
1.41
1.08
1.31
0.829
1.21
1.39
1.64
1.12
1.09
0.985
1.03
1.32

(10,20)
1.05
1.41
1.23
0.707

(20,10)
1.22
1.24
1.39
0.806

∗
(Um
, Un∗ )1
0.827
0.782
0.755
0.677

References[1] Hájek, J., Sidak, Z. & Sen, P. K. (1999). Theory of rank tests, Academic press.
[2] Maesono, Y., Moriyama, T., & Lu, M. (2016). Smoothed nonparametric tests and theirproperties. arXiv preprint arXiv:1610.02145.

1∗
Um
and Un∗ are realized values of the discrete uniform distribution U ∗ (5, 40)

- 8 -

氏

名：蛭川雅之所

属：摂南大学経済学部

講演題目 ： “Nonparametric Estimation and Testing on Discontinuity of Positive
Supported Densities: A Kernel Truncation Approach” [Benedikt Funke (Technical
University of Dortmund)との共著]
講演内容：
本講演では、未知の確率密度関数が所与の点において不連続で（もしくはジャンプが）あると疑われる場合のジャンプ幅の推定および不連続性の検定を取り扱った。密度関数の不連続性に関する問題は、統計学では 90 年代初頭に既に研究がなされていた。一方、経済学においてこの方面の研究は回帰分断デザイン(regression discontinuity design)と関連してごく最近開始された。本講演では、
回帰分断デザインの共変量に用いられる経済変数の多くが非負値をとる点を踏
まえ、共変量の密度関数の台が非負実数全体である場合を想定して非対称ガンマ・カーネルを利用する推定・検定手法を提案した。
まず、ジャンプ幅（＝所与の切断点における密度関数の左側・右側極限の差）の推定法は以下の手順に従う：①ガンマ・カーネルを切断点で左右に分割する；②左右各部分をそれぞれ積分値が 1 となるよう規準化し、二つの正当な切断カーネル関数を作り上げる；③対応する切断カーネル関数を用いて左側・
右側極限をノンパラメトリック推定し、その差をジャンプ幅推定量とする。ただし、この推定量の分散は通常の O(n-1b-1/2)を維持する一方、片側のみの平滑化を用いる手法のため、バイアスは通常の O(b)より遅い O(b1/2)となってしまう。
このような遅いバイアス収束は検定のサイズに悪影響を及ぼすことが懸念されるため、Terrel and Scott (1980)のバイアス修正法を応用してジャンプ幅推定量のバイアスを O(b) に改善した。
次に、ジャンプ幅推定量が漸近的に正規分布に従うことはリアプノフ中心極
限定理を用いて証明される。この結果から、ジャンプ幅推定量に関する t 値を
ジャンプ幅がゼロ（＝切断点において連続）の帰無仮説の下で検定統計量として使用してよいことがわかる。そこで、二種類の漸近分散推定値に基づく検定- 9 -統計量を提案した。なお、これら漸近理論の導出に不完全ガンマ関数の級数近似が積極的に利用されている点も本研究の新規性の一つである。
また、推定・検定を実行する上で重要な問題である平滑化パラメータの選択法も提案した。具体的には、サブサンプリング法に基づく検出力最適化法を本検定に適合させて用いた。加えて、所与の切断点で密度関数が連続であるとの帰無仮説が棄却された場合に密度関数全体を推定する方法にも言及した。特に、切断点の左右で切断カーネル関数を用いて密度推定を行うと、通常のガンマ・カーネルを用いた場合と同様の統計的特性が得られることを証明した。
さらに、モンテカルロ実験では、提案されたジャンプ幅推定量のバイアスがほぼゼロとなる点、および、検出力最適化法によって計算された平滑化パラメ
ータの値を用いた検定統計量はサイズを損なうことなく検出力を改善でき、
McCrary (2008)の検定統計量に比べて高い検出力を得られる点を確認できた。
Angrist and Lavy (1999)の実データを用いた応用例では、McCrary (2008)の検定
では密度関数が連続であるという帰無仮説が棄却できない分断点の幾つかで、提案する検定は仮説を棄却した。これは本検定の高い検出力の裏付けと見ることができる。
最後に、今後の本研究の発展・拡張に触れる。まず、現状切断点は唯一と想定しているが、複数の切断点を特定し各点におけるジャンプ幅ゼロの同時検定法を開発すべきであろう。さらに、切断点を特定せずに密度関数のジャンプの有無を検定する方法の開発にも取り組みたいと考えている。
参考文献：
Angrist, J. D., and V. Lavy (1999): “Using Maimonides’ Rule to Estimate the Effect ofClass Size on Scholastic Achievement,” Quarterly Journal of Economics, 114, 533 575.
McCrary, J. (2008): “Manipulation of the Running Variable in the Regression
Discontinuity Design: A Density Test,” Journal of Econometrics, 142, 698 - 714.
Terrell, G. R., and D. W. Scott (1980): “On Improving Convergence Rates forNonnegative Kernel Density Estimators,” Annals of Statistics, 8, 1160 - 1163.

- 10 -

円周上のカーネル密度推定量とそのバンド幅選択法の漸近的性質
金沢大学人間社会環境研究科
金沢大学経済学経営学系

1 はじめに周期性を持つデータは、円周上の角度として値を取

鶴田靖人
寒河江雅彦

f (θ) に従うとする．カーネル密度推定量は，n1∑Kκ (θ − Θi ).
fˆκ (θ) :=n i=1

る確率変数として扱うことができるので方向データで与えられる．

と呼ばれる．このような方向データを分析することを∫目的とした統計学が方向統計学である．方向統計学では周期性を持つ確率変数 Θ ∼ f (θ) を扱う．ただし，Marzio et al. (2011) は sin モーメント ηj (Kκ ) :=sin(θ)j Kκ dθ を用いて sin オーダーカーネルを提案した．p 次 sin オーダーカーネル関数 Kκ を，

η0 (Kκ ) = 1,

θ ∈ [−π, π)，f (θ) = f (θ + 2π) とする．

ηj (Kκ ) = 0, 0 < j < p,

ηp (Kκ ) ̸= 0,

本稿は、方向データのモデリングとしてよく用いらと定義する．p 次 sin オーダーカーネル密度推定量のれる巻き込みコーシー (WC) 分布をカーネル関数とし

MISE は以下の式となる：()∑∞
1 + 2 j=1 γj2 (κ)ηp2 (Kκ )R f (p)AMISE =
+
,
(p!)2
2πn
(1)
∫
ただし，R(g) := {g(θ)}2 dθ ，γj (Kκ ) := EK [cos(jθ)]

た WC カーネル密度推定量の漸近的性質を与える．数値実験の結果から分布が集中しているという条件の下では WC カーネル密度推定量は良い性質を持つことが分かっている．

．式 (1) の第 1 項と第 2 項はバイアスの 2 乗と分散にまた，新しいモーメントを定義し，p 次オーダーカー
ネル関数のクラスを提案する．このカーネル関数のクラスは，Hall et al.(1987) が与えた 2 次オーダー
カーネル関数のクラスを p 次オーダーカーネル関数に拡張したものとなっている．p 次オーダーカーネ
ル密度推定量は平均二乗誤差 (MISE) の収束レート
が O(n−2p/(p+1) ) で漸近正規性を持つ．このカーネル
関数は，Jones and Foster (1993) の加法型構成法やTerrell and Scott (1980) の乗法型構成法を用いて，2次オーダーカーネルから 4 次以上の高次オーダーカー

対応している．

WC カーネルを以下のように定義する:
1
1 − ρ2Kρ (θ) =, 0 < ρ < 1,2
2π 1 + ρ − 2ρ cos(θ)
ただし，ρ は集中度を調節するパラメータとする．WC
カーネルの特性関数は以下のように定義できる：

ϕp = ρ|p| .
WC カーネルの γj (ρ)はγj (ρ) = ρj となるので，WCカーネルの AMISE は次式で与えられる：′′

AMISEWC [fˆρ (·)] =ネル密度推定量を構成できる．
方向統計学での平滑化パラメータ (集中度パラメータ) κ の選択法として， least square cross validation(CV) 法（Hall et al. , 1987) や Direct Plug-in (PI)
法 (Mardio et al. , 2011) などがある．κ の推定に

{1 − ρ2 }2 R(f )1+.
16nπ(1 − ρ2 )(2)(2) について 1 − ρ2 = h とおくと，次式が得られる：′′
h2 R(f )
1
ˆ
+
.
(3)
AMISEWC [fh (·)] =
16
nπh
∗
∗
κ と同様にして (3) を最小にする最適な h は次式のようになる：

い．我々は CV 法と PI 法の漸近的な性質を導出した．

)1/3
′′
8
h =
, n > 8(πR(f ))−1 . (4)πR(f ′′ )n(3)，(4) より AMISE の収束レートは O(n−2/3 ) となる．

CV 法による κ の推定量 κ̂ の漸近的な収束レートはWC カーネルの最適な収束レート O(n−2/3 ) は同じ

ついての研究は，数値実験による比較実験が主であり，選択法の理論的な性質はほとんど議論されていな

O(n

−1/10

) であるが，PI 法による κ̂ の漸近的な収束レートは O(n−5/14) となる．この結果は PI 法の方がCV 法よりも優れた選択法であることを示している．

∗(

VM カーネルの収束レート O(n−4/5 ) とは異なる． このことは，Sin 型モーメントは MISE の収束レートに対応したカーネル・モーメントとは言えないことを示唆している．また，Di Marzio et al. (2011) は， sin

2 先行研究型 p 次オーダーカーネルは必ずしもバイアスを修正で対称なカーネル関数を Kκ (θ) とする．ただし，κ > 0

きるわけではないと指摘した．つまり，sin 型 p 次オー

は集中度パラメータ (バンド幅の逆数に対応する平滑化

ダーカーネルの次数と MISE の収束レートは必ずしもパラメータ) とする．標本 Θ1 , . . . , Θn は独立同一分布

対応しているわけではない．

- 11 -

漸近的な性質を示す：

3 高次オーダーカーネル密度推定量pκ̂CV /κ∗ −→ 1,

本稿では Kκ (θ) := Cκ−1 (L)L(κ{1 − cos(θ)}) と定義する．モーメントを µl (L) :=∫∞0

L(r)r(l−1)/2 dr，ただし，r = κ{1 − cos(θ)} とする．の偶数 l に対して，

0 < l < p,)(
∑p/2 bp,2t f (2t) (·)
t=1

PI 法とは，ψ4 をカーネル密度推定法で推定し κ∗ を義する：

量の MISE は以下の式で表せる：

d(L)κ1/2,µ20 (L)κpn(5)
ただし，bp,2t ，d(L) は定数である．式 (5) の第 1 項と第2t!

AMISE =

(7) の収束レートは実数直線上の CV 法で推定したバン

推定する方法である．カーネル密度推定量 ψ̂4 (g) を定µp (L) ̸= 0,

を満たすことである．p 次オーダーカーネル密度推定µ2p (L)R(7)

ド幅推定量 ĥCV の収束レートの対応している．

Kκ (θ) が p 次オーダーカーネル関数であるとは，正µ0 (L) ̸= 0, µl (L) = 0,d2n1/10 (κ̂CV /κ∗ − 1) −→ N(0, σCV).

+ψ̂4 (g) := n−2

n ∑
n
∑

Tg(r) (Θi − Θj ),i=1 j=1−1
ただし，Tg (θ) = Cg (S)S(g{1 − cos(θ)}) は p 次オー

ダーカーネル関数，g は集中度パラメータ．

ψ̂4 (g) のバイアスを最小にする g を g∗ と表す：
g∗ := cn2/(p+5) ,2 項はバイアスの 2 乗と分散に対応している．κ = h−2

ただし，c は定数．このとき，ψ̂4 (g) の平均二乗誤差とおけば，(5) は実数直線上で定義されたカーネル密度(MSE) は，推定量の MISE の形によく似ている．実際に，その収束レートは O(n−2p/(p+1) ) となる．

{
O(n−(2p+1)/(p+5) )inf MSE[ψ̂4 (g∗ )] =g>0
O(n−1 )

p < 4,p ≥ 4,となる．したがって，4 次以上の高次オーダーカーネル

4 高次オーダーカーネルの構成法を用いたとき，ψ̂4 (g) はパラメトリックな収束レートでKκ,[p] (θ) := Cκ−1 (L[p] )L[p] (κ{1 − cos(θ)}) は p 次
オーダーカーネル関数を表す．

ある O(n−1 ) を達成する．
また，PI 推定量は，

κ̂PI := β ψ̂4 (g∗ )2/5 n2/5 ,Jones and Foster (1993) に対応する加法型構成法は(8)

次の式となる：で与えられる．κ̂PI は Tg が２次オーダーカーネル関数

p+12
L[p] (r) :=L[p] (r) + rL′[p] (r),ppただし，r = κ{1 − cos(θ)}，L′[p] (r) = dL′[p] (r)/dr.

のとき，以下のような漸近正規性を持つ．

Terrell and Scott (1980) に対応する乗法型構成法
は，2 次オーダーカーネル密度推定量 fˆ の積で与えら(9) の収束レートは実数直線上の PI 法で推定したバン

れる：

レートを比較すると PI 法の方がスピードが速い．

d2
n5/14 (κ̂PI /κ∗ − 1) −→ N(0, σPI).

−1/3fˆκ[TS] (θ) := fˆκ4/3 (θ)fˆκ/4 (θ),[TS]

ただし，bias2f [fˆκ(θ)] = O(κ−4 ), Varf [fˆκ

[TS]

(9)

ド幅推定量 ĥPI の収束レートに対応する．両者の収束

参考文献(θ)] =[1] Di Marzio, M., Panzera, A. and Taylor, C. C.

O(n−1 κ1/2 ) となる．

(2011).

5 平滑化パラメータの選択法

Journal of Statistical Planning and In-

ference 141, 2156-2173.

2 次オーダーカーネル密度推定量の (5) を最小にする[2] Hall, P., Watson, G. S. and Cabrera, J. (1987).
Biometrika, 74, 751-762.

集中度パラメータ κ∗ は以下の式で与えられる:2/5κ∗ = βψ4 n2/5 ,(6)∫ π (r)ただし β は定数，ψr := −π f f (θ)dθ ．本稿では κ∗の推定法として CV 法と PI 法を挙げている．これらの[3] Jones, M. C. and Foster, P.J. (1993). Journal of
Nonparametric Statistics 3, 81-94.
[4] Scott, D. W. and Terrell G. R. Jounal of
the American Statistical Association, 82, 1131-

推定法の定義と漸近的性質について述べる．

1146.

CV 法とは，n2∑ˆf−i (Θi ),CV(κ) = R(fˆ) −n i=1[5] Sheather S. J. and Jones M. C. Jounal of theを最小にする推定量 κ̂CV を求める手法である．κ̂CV の[6] Terrell, G. R. and Scott, D. W. (1980). The An-

Royal Statistical Society. Series B, 53, 683-690.
nals of Statistics 8, 1160-1163.

- 12 -

角度データのための統計モデル
統計数理研究所加藤 昇吾

本報告では，加藤 (2016) のレビュー論文の内容の一部に基づき，角度データに関連した統計的話題を議論した．具体的には，(1) 角度データの解析における問題点，(2) 角度データのための統計手法の基礎，(3) 円周上のコーシー分布とその拡張，という３つのテーマに分けて報告を行った．

1．角度データの解析における問題点
様々な学問分野において，角度として表される観測値が得られることがある．例えば，気象学における風向の観測はその一例である．風向は，西を −π とし，反時計回りを正の向きとすれば，南を −π/2，東を 0，北を π/2 のように角度で表すことができる．つまり，任意の
風向は −π 以上 π 未満の角度 θ，もしくは円周上の点 (cos θ, sin θ)，として表現できる．他には，犯罪が起こる時刻（0 時から 24 時）を記録したデータも同様に，角度の観測と解釈することができる．その他，医学・地震学・動物行動学などの分野においても角度の観測が存在している．
角度のデータには，統計解析をする上で大きな問題がある．それは，このようなデータの解析する上では，統計学が主に対象としている実数値データのための解析手法をそのまま使うことができないという問題である．この問題は，角度には周期性があり，その位相が実数の位相と異なっていることに起因している．例えば，平均や分散などの要約統計量，実数値
データのための確率分布や回帰モデルを角度のデータにそのまま応用すると，不自然な結果を得ることになってしまう．

2. 角度データのための統計手法の基礎
前節で述べた問題点を克服するため，統計学では角度データの解析法に関する研究がなされてきた．本報告においては，角度データのための統計手法の基礎として，角度データのための要約統計量および確率分布について紹介を行った．
要約統計量については，mean direction と mean resultant length を紹介した．これら２つの要約統計量は次のように定義される．確率変数 Θ1 , Θ2 , . . . , Θn を [−π, π) に値をとる角度の標本とする．このとき，この標本の mean direction と mean resultant length は，
( n
)
n
∑
1 ∑ iΘjiΘj
Θ = arg
e
, R=
e
,
n j=1
j=1
とそれぞれ定義される．ここで，i は虚数を表す．直感的に言えば，mean direction は角度
データの平均に相当し，mean resultant length は角度データのばらつきをあらわす要約統計量である．このように定義をすると，それぞれの要約統計量の値が自然に解釈をすることが可能となる．
確率分布については，フォン・ミーゼス分布と wrapped 分布を紹介した．フォン・ミーゼ
ス分布は，確率密度関数

f (θ) =

1
exp {κ cos(θ − µ)} ,2πI0 (κ)−π ≤ θ < π;

- 13 -−π ≤ µ < π, κ ≥ 0,で定義される円周上の確率分布である．ここで，Ij は位数 j の第一種変形ベッセル関数をあらわす．フォン・ミーゼス分布は，‘円周上の正規分布’ とも呼ばれ，角度データの統計解析において中心的な役割を果たしてきた．実際，フォン・ミーゼス分布には，指数型分布族に属する，最尤推定量が陽な形で表現できる，などの利点がある．一方，正規化定数が特殊関数を含む形となる，よく知られた ‘尺度’ の変換に閉じていない，再生性が成り立たない，などの正規分布が持つ扱いやすい性質が成立しない問題点も知られている．
次に wrapped 分布の紹介を行った．X を実数に値をとる確率変数とする．また，確率変数
Θ を以下のように定義する：
X ≡ Θ (mod 2π).
このとき，Θ の従う分布を，分布 X の wrapped（巻き込み）分布という．wrapped 分布は，実数上の分布から直接的に得ることが可能であり，分布の単峰性を保つなどの利点があるが，一般には確率密度関数が無限級数の形であらわされる等の問題点もある．

3. 円周上のコーシー分布とその拡張
本報告の最後のテーマとして，円周上のコーシー分布とその２変量拡張の紹介を行った．
円周上のコーシー分布は，確率密度関数

f (θ) =

1 − ρ21,
2π 1 + ρ2 − 2ρ cos(θ − µ)−π ≤ θ < π;−π ≤ µ < π, 0 ≤ ρ < 1,で定義される．円周上のコーシー分布は，再生性が成り立つこと，メビウス変換に関して閉
じていること，無限分解可能であること，などの扱いやすい性質を持つことが知られている（例えば，McCullagh, 1996）．一方，最尤推定に関しては，n ≤ 4 のときは最尤推定量の陽な表現が知られているものの，一般には Kent & Tyler (1988) のアルゴリズムなどを用いて数値的に求める必要がある．
最後に，円周上のコーシー分布を２変量分布へと拡張した Kato & Pewsey (2015) の分布を紹介した．この分布の確率密度関数は，特殊関数や無限級数などを含まない陽な形で表現することが可能である．他にも，分布のパラメータの解釈が容易である，周辺分布と条件付分布がともに円周上のコーシー分布となる，効率の良い疑似乱数の発生法が存在する，特性関数や相関係数を簡潔な形で表現できる，などの多くの扱いやすい性質を持っている．パラメータの推定に関しては，モーメント推定に関しては推定量を陽にあらわすことが可能である．最尤推定については一般に数値的に求めることが必要であるが，モーメント推定値を初期値として使用すれば，推定アルゴリズムが早く収束することが経験的に知られている．
参考文献[1] 加藤昇吾 (2016). 「円周上のコーシー分布と関連した統計モデル」. 『日本統計学会誌』, 46,
85–111.
[2] Kato, S. & Pewsey, A. (2015). A Möbius transformation-induced distribution on the torus.
Biometrika, 102, 359–370.
[3] Kent, J.T. & Tyler, D.E. (1988). Maximum likelihood estimation for the wrapped Cauchy
distribution. Journal of Applied Statistics, 15, 247–254.
[4] McCullagh, P. (1996). Möbius transformation and Cauchy parameter estimation. The Annalsof Statistics, 24, 787–808.

- 14 -

The Simultaneous Multivariate Hawkes-type
Point Processes and their application to
Financial Markets
Naoto Kunitomo (Meiji University)
joint work with
Daisuke Kurisu, Yusuke Amano, and Nao Awaya (University of Tokyo)
In economic and financial time series we sometimes observe large jumps. Although
they are relatively rare events, they often have significant influence not only on afinancial market but also several diﬀerent markets and macro economies. By using
the simultaneous Hawkes-type models we introduce, which are multivariate point
processes, it is possible to analyze the causal eﬀects of large events in the sense of
the Granger-non-causality (GNC) and the instantaneous Granger-non-causality(IGNC). We investigate the financial market of Tokyo and other markets, andapply the Granger non-causality tests. We have found several important empiricalfindings among financial markets and macro economies.
We divide the continuous observation period [0, T ] to the discrete observation
periods Iin = (tni−1 , tni ] (i = 1, · · · , n). The initial time is tn0 = 0 and we interpret
Iin as the i−th day, but it may be possible to have more finer observations. Let theobservable price process be Pj (t) (j = 1, · · · , d ; tni−1 < t ≤ tni , i = 1, · · · , n) andin s ∈ Ii we consider the (negative) log-returns of prices Yjn (s) (tni−1 < s ≤ tni ) be
Yjn (s) = − log[Pj (s)/Pj (tni−1 )] (j = 1, · · · , d; i = 1, · · · , n) . Let the first stopping
time when Yjn (s) exceeds the threshold uj in s ∈ Ii be τjn (i, 1). Define Xjn (s) =Yjn (s) for s ∈ tni−1 ≤ s ≤ τjn (i, 1) and Xjn (s) = Xjn (τjn (i, 1)) for s ∈ [τjn (i, 1), tni ].
We define the simple counting processes Njn∗ (s, uk ) by the number of stopping
times that Xjn (s) exceed uj (j = 1, · · · , d) for a particular j but not for otherk ̸= j by the time s. (For the simplicity, we can assume that the jumps of thecounting process Njn∗ (s, uk ) can occur once at each tni of the end of each intervals(tni−1 , tni ].) Then we notice that as n −→ ∞ the interval length goes to zero,that is, maxi=1,···,n |tni − tni−1 | −→ 0 and the simple counting process Njn∗ (s, uk )converges to Nj∗ (s, uk ). The resulting counting process can be interpreted as thelimiting process in the high frequency asymptotics.
The simple point processes Nj∗ (t) (j = 1, · · · , d) satisfies the standard conditions that as ∆t → 0 P (Njn∗ (t+∆t, uj )−Njn∗ (t, u) = 1|Ftn ) = λnj (t, uj )∆t+op (∆t),
P (Njn∗ (t+ ∆t, uj )− Njn∗ (t, u) > 1|Ftn ) = op (∆t), P (Nkn∗ (t +∆t, uj ) −Njn∗ (t, uj ) ≥1|Ftn ) = op (∆t) for k ̸= j, where Ftn is the σ−field generated by the informationn∗(s, ujk ) by the number of stopat t. Also we define the simple proint processes Njk
nping times that Xj (s) exceed uj (j = 1, · · · , d) for a particular j and Xkn (s) exceeduk (k = 1, · · · , d; k ̸= j) for a particular k and other Xln (s) (l ̸= j, k) do not exceedn∗(t, ujk ) withul by the time s in Iin . Then we introduce the point processes Njkn∗n∗n∗
co-jumps of Nj and Nk by P (Nj (t + ∆t, uj ) − Nj (t, uj ) = Nk (t + ∆t, uk ) −n∗n∗Nkn∗ (t, uk ) = 1|Ft ) = λn∗jk (t, ujk )∆t + op (∆t), P (Nj (t + ∆t, uj ) − Nj (t, uj ) >1|Ft ) = op (∆t), P (Nkn∗ (t + ∆t, uk ) − Njn∗ (t, uk ) > |Ft ) = op (∆t) for k ̸= j.

- 15 -

When we have co-jumps of two point processes, we define the point process(1)n∗
Njn (s, uj ) = Njn∗ (s, uj ) + Nj,k(s, ujk ) (j ̸= k; j, k = 1, · · · , d)and the corresponding conditional intensity functions are given by(2)n∗λnj (t, uj ) = λn∗j (t, uj ) + λj,k (t, ujk ) .

Then the resulting processes correspond to as the marginal point processes of thevector point process. We use the self-exciting form of intensity function as
(3)

∗
λn∗
j (t, u) = λj0 +

∫ t

p
∑

−∞ i=1cji (X n (s−))gji (t − s)dNJn∗(s, u) ,i

where the index set is defined as Ji = i (i = 1, · · · , d), Ji = 1, 1 + (i − d) (i = d +1, · · · , 2d − 1), · · · , Jp = 1, · · · , d. Sequentially define Nin∗ (s, ui ) = Nin∗ (s, u) (i =n∗n∗n∗1, · · · , d), Nd+1
(s, u) = N1,2
(s, u), · · · , and Npn∗ (s, u) = N1,···,d(s, u) with the intenn∗sities λj,k (t, ujk ). In this formulation we use the damping function as gji (t − s) =e−γji (t−s) and the impact function as C(X) = (cji (x)). Since we are interested inlarge jumps, it is important to use the probability function in the tails. We use
the generalized Pareto (tail) probability function for x > uj (j = 1, · · · , d) as
[

(4)

j
x
1 + σjξ(s)]−1/ξjP (Xjn (s) > x|Xjn (s) > uj , Fs ) = []−1/ξjj
1 + σjξ(s)uj[=]−1/ξjξj1 + ∗ (x − uj )σj (s),

we set σj∗ (s) = ξj uj + σj (s).
References
Florens, J-P. and D. Fougere (1996) “Noncausality in Continuous Time,”Econometrica,
64, 1195-1212.
Grothe, O., V. Korniichuk and H. Mannera (2014)“Modeling multivariate extreme events using self-exciting point processes,” Journal of Econometrics, 182,
269-289.
Hawkes, A. G. (1971) “Point Spectra of Some Mutually Exciting Point Processes,”
Journal of the Royal Statistical Society. Series B, 33-3, 438-443.
Kunitomo,N., A. Ehara and D. Kurisu (2016) ”A Causality Analysis of Financial
Markets by Multivariate Hawkes-type Models,”(in Japanese), CIRJE, Discussion
Paper J-278, CIRJE, University of Tokyo.

- 16 -

ウェーブレット分散と参照形式を利用した商品先物の暴落予想
上智大学大学院理工学研究科数学領域・院

本山

要

上智大学理工学部情報理工学科 加藤剛

研究目的と背景1

本研究は, ウェーブレット解析の中でも特にウェーブレット分散と呼ばれる概念を活用し, 東京商品取引所(TOCOM) における原油先物価格について, 暴落の前兆現象をとらえることを目的としたものである. 最終的には,暴落特有の特徴をもった参照図を作成し, 新たなデータとの参照形式で暴落の注意喚起や警報に使用することを目標としている.
金融市場や金融商品売買に関する報道が絶えず行われている. 2008 年のリーマンショック, その後の中国の不動産バブル, ギリシャ危機,Google の時価総額が世界首位獲得, 次期アメリカ大統領が決定し金融市場に影響を及ぼす等々. 世界のあらゆる出来事が金融市場に影響を与え複雑化している. また, 日本国内に注目すれば,2001 年から確定拠出年金制度の導入,NISA(少額投資非課税制度), あるいは仮想通貨ビットコインの出現になどより, 個人が資産運用や投資に関わる機会が増えている現実がある. 先物取引の適切なリスク指標をつくることは, 機関投資家だ
けでなく個人投資家にも市場参加を助長し, 市場の活性化につながる. しかし,2009 年 5 月に相場の急激な変動がおきた際の緊急処置として新しくサーキットブレーカー制度が採用されたこともあり, 理論的裏付けがあるリスク指標に関する研究は少ない. なお, 石油を選択した理由は, 価格の下落が頻発しているためである.

参照形式による暴落予想

2

参照形式とは, ある特徴を有する状況の情報を参照図として保存しておき, 将来, 参照図と似ている情報をもった状況が現れた場合, 同じような状況が起こる可能性があると予測する方法である. 本研究は, この方法を東京商品取引所 (TOCOM) が扱う原油先物価格の暴落予想に応用したものである. なお, 値動きの状況を表す情報として様々なパラメータが存在するが, ウェーブレット分散の有効性を確かめるため, ウェーブレット分散の時間的推移のみを値動きの状況を表す情報としている.

3

データ解析

3.1

参照図の作成本研究では, 暴落を原油におけるサーキットブレーカー発動条件 2000 円の 75% である 1500 円を上回る下落がある場合の値とし, 時間の間隔とあわせて次のように定義する.
定義 1 暴落時刻
次の不等式を満たす Xt が存在するとき, 時刻 t を暴落時刻とする.

Xt−(216 +214 ) − Xt > 1500
幅を 216 , すらしを 214 , 着目スケールを 210 とした. 前処理した 2014 年 9 月 1 日 0 時 0 分 0 秒から 2015 年 12月 30 日 23 時 59 分 59 秒までの原油の約定値段のデータ数は 20, 768, 150 であったので, 1263 個のウェーブレッ
ト分散が得られた. 得られたウェーブレット分散の個数と同じ 1263 回暴落時刻を調べた. その結果, 1263 回中

111 回が暴落時刻に該当することがわかった.
ウェーブレット分散から 4 点のパターンを構成するイメージを図 1 に示す. 図 1 で, 垂線が下りている時刻を暴落が起きたか否かを調べた時刻 T とする. 時刻 T に対して, 1〜4 の範囲でそれぞれ計算された 4 つのウェーブ- 17 -

図1

暴落があるか調べた時刻とそれに対応するウェーブレット分散の対応図

レット分散の値のパターンが対応する. これを, ずらし 214 でずらしていくことにより, 1260 個のパターンが得られる.
新たなデータから計算したウェーブレット分散にもとづくパターンと類似したものが参照図に存在するとき, 暴落に対する注意喚起や警報を行う. パターンの類似度は距離で測る.

3.2

参照図の作成以上の議論を踏まえて, 参照図の作成方法を説明する.

1. 使用するデータの期間を決定
2. パターンを構成するウェーブレット分散の個数を決定 (探索的に 4 個から 10 個までとした)
3. 1 つの参照図を構成するパターン数を決定 (4 個から 6 個までとした)
4. パターンを作成する方法を決定(a) quantile based
3.3(b) cluster based

結論

次の 2 点で作成した参照図を評価した.

• 信頼性 参照図がもつ, いずれかのパターンに近いパターンが見つかった時に, 実際に暴落が起こる割合• 性能 参照図がもつパターンが, 全体の暴落の回数のうち, 近いパターンとして検出できた割合
結論として, 信頼性と性能の両方の値が十分に高いといえる参照図はみつからなかった. しかし, 信頼性は高く性能が低い参照図, 逆に, 信頼性は低いが性能が高い参照図, さらに細かく見れば, 信頼性が高いパターンを発見できた. したがって, これらを組み合わせることで, より実用的な参照図を作成できる可能性がある.

主な参考文献[1] Nason, G. P., ”Wavelet Methods in Statistics with R”, Springer, (2008).
[2] 東京商品取引所ホームページ, http://www.tocom.or.jp/jp/- 18 -

リスク計測モデルのバックテスト法について
塚原 英敦（成城大学 経済学部）

1

はじめに — 問題の背景

リスク計測モデリングとは，ファイナンス・データに対する統計モデルを設定し，パラメータ推定を行い，リスク尺度の推定値を得るという一連の作業を指す．そして，(1) 現在
用いている統計モデル・推定法のパフォーマンスをモニターする．(2) 統計モデル・推定法の相対的な比較を行うというのが，一般に言われているバックテストの 2 つの目的である（McNeil et al. [7]）．前者の側面を重視して述べれば，バックテストはモデルと推定法の妥当性をチェックする過程で用いられる手法であり，金融リスク管理を適切に行うためには不可欠である．バリュー・アット・リスク（VaR）の場合，よく用いられているバックテスト
法は，いわゆる VaR 超過数に基づくものであり，直観的にも容易に理解できる方法である（Campbell [2]）．
Consultative Document of Basel Committee on Banking Supervision (October 2013) では，2008・9 年の金融危機における VaR の欠陥を指摘した上で，規制資本量決定のためのリスク尺度として，VaR に代わって期待ショートフォールを採用することを提案している．
しかし，同じ文書において，バックテスト法としては VaR に基づく方法の使用を求めている．その理由としては，多くの人々が VaR は期待ショートフォールや他のリスク尺度よりも
バックテストするのが簡単であると主張していることが挙げられる．その根拠として，期待
ショートフォールに対する既存のバックテストは，帰無分布についてのパラメトリックな仮定に基づいていること，そして検定統計量の帰無分布に対する漸近的な近似が必要であることが指摘されている．この議論によれば，リスク尺度のバックテスト可能性（backtestability）とは，モデルによらない（model-free）やり方で，小標本に基づいてバックテストができることと解釈されるが，VaR の場合はこれが満たされる非常に稀な例であり，他のリスク尺度についてこの条件を課すのは無理な注文であると言わざるを得ない．

2

顕在化可能性

近年，バックテスト可能性の数学的な表現として，顕在化可能性（elicitability）という概念が注目を集めている（Gneiting [6]）．大雑把に言えば，統計的汎関数 T (F ) が分布族 F に対して顕在化可能（elicitable）であるとは，ある評点関数 S が存在して，すべての F ∈ Fに対して t = T (F ) が t 7→ EF [S(t, L)] を最小化するただ 1 つの値であることである．ここで，L は分布 F に従う確率変数である．VaR が顕在化可能であることは容易にわかるが，期
待ショートフォールやそれを含むクラスである歪みリスク尺度（distortion risk measures）
はこの顕在化可能性の条件を満たさないことが最近示された（Ziegel [8]）．さらに，整合性かつ法則不変性をもつリスク尺度で顕在化可能なものは，expectile しかないことも証明された（Bellini and Bignozzi [1]）．ただし，expectile は直観的な理解に欠けるという短所もあり，実用には程遠いであろう．
顕在化可能性はバックテストの目的 (2) のためにはある程度有効である一方，モデルのモニタリングという目的 (1) との関連性は明らかではない．例えば，平均は顕在化可能であるが，平方誤差に基づいて，モデルに依存しない形でバックテストすることは，漸近近似を用
いるかパラメトリックな分布を仮定しない限り不可能である．

- 19 -

3

Davis の較正性観測される損失と補助変数の系列 (Ln , X n )n∈N に対して，モデル P = {P z : z ∈ Z } は(Rd+1 )∞ 上の確率分布族であり，それに対して条件付分布列 F1z (l) := P z (L1 ≤ l), Fnz (l) :=P z (Ln ≤ l | Fn−1 ), n = 2, 3, . . . が定義される．このとき，統計的汎関数 T に対して，較正関数
ϕ : R2 → R と，狭義単調増加な可予測正規化列 b = (bn ) で，limn→∞ bn = ∞, P z -a.s., ∀z ∈Z を満たすものが存在し，
n
1 ∑
lim
ϕ(Lk , T (Fkz )) = 0, P z -a.s., ∀z ∈ Z
n→∞ bnk=1が成り立つとき，T は P で (ϕ, b) 較正可能と呼ばれる（Davis [3]）．例えば，VaRkα に対しては，ϕ(l, t) = 1[t,∞) (l) − α, bn = n ととればよいことが大数の法則から明らかである．
実用上は，観測値 (L1 , X 1 ), . . . , (Lk−1 , X k−1 ) が得られ，T (Fkz ) の予測値 τk が何らかのアルゴリズムで計算されたとき，この予測値の “妥当性”を Jn (L1 , . . . , Ln , τ1 , . . . , τn ) =∑nb−1
nk=1 ϕ(Lk , τk ) が十分 0 に近いかどうかで判断できる．これは，Dawid [4, 5] による逐次予測（prequential）原理「予測評価の妥当性基準は，損失の実現値と予測の実際値のみに依存すべきである（予測値がどう計算されたかは無関係）」を満たしている．
Davis は非常に広い条件の下で VaR が較正可能なのに対して，期待ショートフォールのような平均型の汎関数については，較正可能のためにはかなり厳しい条件が必要であることを示した．しかし，較正可能性もまた漸近的な概念であり，バックテスト可能性を表現する
十分に満足できる答であるとは言えないであろう．
以上の議論から，バックテスト可能性という概念の数学化は非常に困難な，あるいは不可能な問題のように見受けられる．

参考文献[1] F. Bellini and V. Bignozzi. On elicitable risk measures. Quantitative Finance, 15:725–
733, 2015.
[2] S. D. Campbell. A review of backtesting and backtesting procedures. Journal of Risk,
9:1–17, 2007.
[3] M. H. A. Davis. Verification of internal risk measure estimates. Statistics & Risk
Modeling, 33:67–93, 2016.
[4] A. P. Dawid. Present position and potential developments: some personal views.
statistical theory. the prequential approach (with discussion). Journal of the Royal
Statistical Society (A), 147:278–292, 1984.
[5] A. P. Dawid. Prequential analysis. In C. B. Read S. Kotz and D. L. Banks, editors, Encyclopedia of Statistical Sciences, Update Volume 1, pages 464–470. Wiley-Interscience,
1997.
[6] T. Gneiting. Making and evaluating point forecasts. Journal of the American Statistical
Association, 106:746–762, 2011.
[7] A. J. McNeil, R. Frey, and P. Embrechts. Quantitative Risk Management: Concepts,
Techniques, and Tools. Princeton University Press, Princeton, New Jersey, second
edition, 2015.
[8] J. F. Ziegel. Coherence and elicitability. Mathematical Finance, 26:901–918, 2016.

- 20 -

交互作用を含む観察データにおける傾向スコア解析手法による結果の差異
加葉田 大志朗 1,2), 山本 紘司 1), 新谷 歩 2)
1) 大阪大学 大学院医学系研究科 臨床統計疫学寄附講座
2) 大阪市立大学

大学院医学研究科

医療統計学講座

１．はじめに観察研究により得られたデータの解析では、交絡等のバイアスへの対処が必須となる。特に近年医学分野においては、交絡の除去を目的として、傾向スコアを利用した解析手法が普及している。傾向スコアは治療に割り当てられる確率を背景情報などから計算したものである。現在提案されている主流な傾向スコア解析としては、類似した傾向スコアで対照群を選ぶマッチング法、傾向スコアの逆数を用いて重み付けする
IPTW(Inverse Probability of Treatment Weighted)法、IPTW法と多変量回帰モデルを組み合わせたDR(Doubly Robustness)法がある。特にマッチング法は解釈がしやすく、臨床で利用されることも多い。しかしデータによっては各手法も使い分けが必要となるが、その詳細は知られていない。

2. ケーススタディ
本研究では、舌癌患者における予防的リンパ郭清の効果を検証する目的で得られた
実際の観察研究データの解析より着想を得た。この解析では、傾向スコアによるマッチング法と、3因子交互作用項を含む多変量回帰モデルを用いて解析を行った。その結果傾向スコアでは治療が有効であることが示されたが、多変量回帰モデルにおいては治療効果と年齢・腫瘍サイズの間に交互作用が確認された。つまり治療効果は術前の年齢と腫瘍サイズにより変化し、中年層(40〜60歳)で腫瘍サイズが5mm以上の症例において治療は有効であることが示唆された。一方、高齢あるいは若年層で5mm以上の症例においては、予後を悪化させる可能性があることが示唆された。
この解析では感度解析として、IPTW法とDR法を用いた解析も行ったが、どちらもマッチング法のように治療効果があるという結果は得られなかった。

3. 検証内容
上記の観察研究データ解析におけるマッチング法は、傾向スコアが最も近い症例同士をマッチングするSimple nearest neighbor matchingを用いた。当該手法による解析では基本的にATT(Average Treatment effect for Treated)、主に治療が行われた症例における治療効果を推定することが知られている1)。今回のマッチング法を用いた解析では、中間年齢層で腫瘍サイズが大きい、つまり臨床的に最も治療をうける確率が- 21 -

高い症例における平均的な治療効果が反映されたと考えることができる。
IPTW法とDR法はマッチング法とは異なり解析に用いる症例数が減少しないため、最近では臨床論文でも用いられることが増えてきた。しかし基本的にIPTW法とDR法はATE(Average Treatment Effect)すなわち、集団全体の平均的な効果を検証している。
そのため当事例のように治療効果が対象症例の背景によって変化する状況においては、集団全体の平均的治療効果を推定すると、治療を受けた症例における治療効果が十分に検出されないことが考えられる。またIPTW法とDR法においても重みを調整することで、ATTと同様の解析結果を得られることが提唱されているが2)、実際にそれらを用いた臨床論文は目にすることが少ない。
傾向スコアを用いた手法比較論文においてはATEに着目した報告が行われることが多く、当事例のように、治療効果が交互作用を持つデータにおける解析手法については言及されることが少ない。しかし医学領域において治療効果が患者背景により変化す
ることは多くのケースで考えられ、その影響を考慮して解析手法を選択することは重要であると考える。
本研究では当事例のように治療効果との交互作用を含むデータにおいて、各種傾向
スコア解析手法を適用した際の結果の違いを検証し、目的に応じた適切な解析手法を提案することを目的とし、当事例に類似した擬似データを発生させ、モンテカルロシ
ミュレーションにより各種傾向スコア解析手法を適用した結果の比較を行った内容を報告した。

1) Stuart, E. A., Matching Methods for Causal Inference: A Review and a Look
Forward. Statistical Science : A Review Journal of the Institute of Mathematical
Statistics 25(1): 1–21, 2010.
2) Hirano, K., and Imbens, G. W., Estimation of Causal Effects Using Propensity
Score Weighting: An Application to Data on Right Heart Catheterization. HealthServices and Outcomes Research Methodology 2(3–4): 259–278, 2001.

- 22 -

診断研究における検査性能の比較
高橋 佳苗 （大阪大学医学部附属病院）
山本 紘司 （大阪大学大学院医学系研究科）
１．はじめに医療において、スクリーニング検査は疾患の早期発見・早期治療のために重要な役割を果たしている。スクリーニング検査の性能をあらわす指標としては、感度（sensitivity：SE）、、陽性適中率（positive predictive value：PPV）および陰性的中特異度（specificity：SP）
率（negative predictive value：NPV）がしばしば用いられている。SE は罹患している対象者をスクリーニング検査によって正しく罹患していると判断する確率であり、SP は罹患していない対象者を正しく罹患していないと判断する確率である。一方、PPV は対象者のスクリーニング検査の結果が陽性であるとき、当該対象者が実際に疾患に罹患している確率であり、NPV はスクリーニング検査の結果が陰性であった対象者が実際に疾患に罹患していない確率である。SE、SP、PPV、NPV は臨床的に有用であり、患者の治療方針の決定に影響を及ぼす。
2 つのスクリーニング検査法の性能を比べる目的で上記の指標を比較する際、SE およびSP については McNemar 検定が一般的に用いられている。しかし、PPV および NPV につ
いては主流となる解析手法は確立しておらず、いくつかの手法が提案されているのみである(Leisenring et.al, 2010, Moskowitz and Pepe, 2006, Kosinski, 2013)。さらに、提案されているいずれの手法も、中心極限定理が用いられており、対象者の少ない小規模臨床試験への適用に問題がないかは未だ検討されていない。
さらに、スクリーニング検査の性能については、本来 1 つの指標のみで判断するよりも、複数の指標から総合的に判断されるべきである。2 つのスクリーニング検査間で複数の指標を比較する際は、すべての指標の優越性を確認するよりも、1 つの指標の優越性と他の全ての指標での非劣性を組み合わせたデザインの方が実臨床の場に即していると考えられる。
そこで本講演では、①PPV および NPV の比較に対する解析手法を小規模臨床試験へ適用した場合の性能評価、②いずれか 1 つの指標の優越性が認められれば他の指標は非劣性を確認するのみでよいとする、優越性と非劣性を組み合わせた多重エンドポイントの解析手法(Bloch et al, 2007, Muscha and Turan, 2012) の性能評価について報告した。
２．方法および結果①PPV および NPV の解析手法を小規模臨床試験へ適用した場合の性能評価
多項分布をもとに中心極限定理とデルタ法を利用したもの（恒等変換、対数変換、ロジッ
ト変換）、回帰分析の枠組みを利用したもの（一般化スコア統計量、重み付き一般化スコア- 23 -

統計量）の解析手法がすでに提案されている。Kosinski (2013)で報告されている方法を用
いて、表 1 の設定でシミュレーション用データを生成し、症例数を 10 例から 50 例とした
場合の各手法の Type 1 error rate 及び検出力を求めた。
2 つの診断検査間の PPV、NPV に差がないとしたとき(図 1)、名目有意水準に最も近い
検出力となったのは重み付き一般化スコア統計量を用いた方法であった。重み付き一般化
スコア統計量を用いた方法は、2 つの診断検査間の PPV、NPV に差があるとしたとき(図 2)の検出力も高く、小規模臨床試験においては最も性能の良い方法と考えられる。
②優越性と非劣性を組み合わせた多重エンドポイントの解析手法の性能評価
表 2 に示すシナリオについて、①と同様に Kosinski (2013)の方法を用いてシミュレーシ
ョン用データを作成した。有病割合は 0.4、非劣性マージンは 0.1、名目有意水準は 0.05、症例数は 50 例から 500 例とした。SE および SP の解析には McNemar 検定を、PPV およ
び NPV の解析には重み付き一般化スコア統計量を用いた方法を適用した。
Bloch et al. (2007)の方法と Muscha and Turan (2012)の方法間で結果に大きな違いは見られなかった。また、優越性と非劣性の両方を評価するデザインの検出力が優越性の評価のみと比較して大きく変わらない場合もあることから、いずれか 1 つの指標の優越性が認められれば他の指標は非劣性を確認するのみでよいとする、優越性と非劣性を組み合わせた
多重エンドポイントの枠組みでデザイン・解析することの有用性が示唆された。

参考文献
 Bloch, D. A. et al. A combined superiority and non-inferiority approach to multiple
endpoints in clinical trials. Statistics in medicine 2013, 26:1193-1207.
 Kosinski, A. S. A weighted generalized score statistic for comparison of predictive
values of diagnostic tests. Statistics in medicine 2013, 32:964–977.
 Leisenring W, Alonzo T, Pepe MS. Comparisons of predictive values of binary medical
diagnostic tests for paired designs. Biometrics 2000; 56:345–351.
 Mascha, E. J. and Turan, A. Joint hypothesis testing and gatekeeping procedures.

Anesthesia and Analgesia 2012, 114:1304-1317.
 Moskowitz, C. S. and Pepe M. S. Comparing the predictive values of diagnostic tests:
sample size and analysis for paired study designs. Clinical Trials 2006; 3:272–279.

- 24 -

複数の二値評価変数をもつ臨床試験における検定法

石原
山本

拓磨
紘司

(大阪大学大学院・医学系研究科)(大阪大学大学院・医学系研究科)

臨床試験の検証的段階では，薬剤や治療処置による治療効果を評価するために複数の評価変数を主要評価変数として用いる場合がある．本発表では特に複数の評
価変数のうち，少なくとも 1 つの変数について優越性の根拠が得られ，その他すべての変数に対しては非劣性の根拠が得られるときにのみ治療効果があると判断す
るような試験を考える．ここでは簡単のため 2 群比較を想定し，2 つの群 (i = 1, 2)の p 個の評価変数に対して，次のような複合帰無仮説 H0 を考える．
{} {}H0 : max µj ≤ 0 ∪ min (µj + ϵj ) ≤ 0
1≤j≤p

1≤j≤p

(1)

ここで µj は，第 i 群における j 番目の評価変数の処置効果を µij (i = 1, 2; j =1, ..., p) としたときに µj = µ1j − µ2j で表される処置効果の群間差であり，ϵj (j =1, ..., p) は ϵj ≥ 0 の非劣性マージンである．
このとき Bloch et al.（2001），Perlman and Wu（2004），Nakazuru et al. (2014)
は上記の仮説（1）に対する検定手順を提案した．
これらの研究では，処置効果の差は連続変数として議論されている．そこで本発表では，特に評価変数が二値変数である場合において，第 i 群の j 番目の評価変数に対する標本比率を π̂ij (i = 1, 2; j = 1, ..., p)，第 j 番目の評価変数の標本比率の群間差を ∆j = π̂1j − π̂2j としたとき，次の場合に帰無仮説（1）を棄却する検定手
順を次のように提案した（有意水準 α）：−1

χ2 ≡ ∆t Σ̂t ∆ > χ2α and∆j + ϵj√> Zα for j = 1, ..., p.
π̂1j (1−π̂1j )π̂2j (1−π̂2j )+n1
n2
ただし，∆t = (∆1 , ..., ∆p )t であり，χ2α 及び Zα はそれぞれ χ2 分布及び標準正規分
布の上側 α パーセント点である．
提案した検定手順の性能についてはモンテカルロシミュレーションを用いて検出力と Type I error rate を評価した．シミュレーションや実際のデータ解析結果を報告した．
参考文献- 25 -

Nakazuru, Y., Sozu, T., Hamada, C. and Yoshimura, I. (2014). A new procedureof one-sided test in clinical trials with multiple endpoints. Japanese Journalof Biometrics, 35, 17-35.
Perlman, M.D. (1969). One-sided testing problems in multivariate analysis. The
Annals of Mathematical Statistics, 40, 549-567.
Perlman, M.D. and Wu, L. (2004). A note on one-sided tests with multiple
endpoints. Biometrics, 60, 276-280.
Bloch, D.A., Lai, T.L. and Tubert-Bitter, P. (2001). One-sided tests in clinical
trials with multiple endpoints. Biometrics, 57, 1039-1047.
Bloch, D.A., Lai, T.L., Su, Z. and Tubert-Bitter, P. (2007). A combined superiority and non-inferiority approach to multiple endpoints in clinical trials.
Statistics in Medicine, 26, 1193-1207.

- 26 -

R.A.Fisher 以後の判別分析の新理論と遺伝子解析の新手法
成蹊大学

経済学部

特任教授

新村秀一

Fisher は「Fisher の仮説」を考えることで、計算機の助けを借りないで判別理論を確⽴した。しかし、Fisher の仮説を検定する良い統計量がないこと、誤分類確率や判別係数の標準誤差が定式化できなかったこと、などの重要な事実に多くの研究者や利⽤者は注意を払わなかった。Fisher も指摘する通り、「Fisher の仮説を満たさないデータに Fisher の LDF を適⽤してはいけない」が、それが分からないため多くの分野に適⽤され、⼀⾒成果を上げてきた。しかし⼤学卒業後の 1971 年に「⼼電図の⾃動解析システムの診断論理を LDF と QDF
で 4 年間アプローチして、医師の開発した枝分かれ論理に⻭が⽴たなかった」。この経験から、医学診断や、各種の格付けや、試験の合否判定を得点合計で判別すると、線形分離可能（MNM=0）であっても誤分類確率が最⾼ 20%になることが実証研究で分かった。即ち既存の判別理論には、世界中の誰も研究していない 4 つの深刻な瑕疵があることを判別分析の多くの実証研究で⾒つけ、数理計画法による 5 個の最適判別関数を開発し解決してきた。
また、判別分析は推測統計⼿法でないので、「⼩標本のための 100 重交差検証法（新⼿法 1）」を開発し、 誤分類確率と判別係数の 95%信頼区間を求めることができた。また、検証標本の平均誤分類確率（M2）が最⼩モデルを Best モデルとして選ぶと、Vapnik が提起した汎化能⼒に優れたモデルを選ぶことになる。各⼿法の全ての変数の組み合わせモデルで Best
モデルを選び、8 種の異なった LDF（改定 IP-OLDF、改定 LP-OLDF、改定 IPLP-OLDF、
H-SVM、SVM4、SVM1、ロジスティック回帰、Fisher の LDF）の Best モデルの中で最⼩の M2 を⽐較して、多くの場合で改定 IP−OLDF の Best モデルの M2 が最⼩であった。
Vapnik の定義した汎化能⼒は、単に固定された p 変数の LDF での話であることが理解されていない。また、新⼿法 1 は、LOO 法よりはるかに優れている。以上で 4 種の深刻な問題を全て解決したが、2015 年の富⼭での科研費シンポジュームで「世界的に著名な
Microarray データが公開されているのを知り、それらを判別すると、世界中の統計や医学研究者が⾼次元空間の分析と称し、特に Feature Selection の研究を 15 年以上⾏ってきたが(問題 5)が、改訂 IP-OLDF で解決した（Matroska Feature Selection Method、新⼿法 2）。
そして Microarray データは、MNM=0 になる⼩さな部分空間の排他的な和集合になっていることが分かった」。即ち問題 5 の認識と同時に、僅か 41 ⽇間で新⼿法 2 を開発した。これによって、新⼿法 1 と新⼿法 2 で Fisher 以後の新しい判別理論を確⽴した（New Theory
of Discriminant Analysis After R. Fisher, Springer(2017)）。
参加者と以下の点について真摯に議論したい．

- 27 -

1) Fisher は Fisher の仮説に基づいて判別分析理論を確⽴したが，限界を知っていた．⽥辺の指摘[48]の他，検証にアヤメのデータを⽤いている点，仮説を満たさない場合に QDFが提案されていることが重要である．
2) Fisher 以後，RDA や LASSO 等の分散共分散に基づく理論が開発され，これらが Fisherの後継者と考えられているが間違いである．Fisher 理論を適⽤してはいけない医学分野
で，Cox は Cox 回帰やロジスティック回帰を提案したが，彼こそが Fisher の正当な後継者である．そして Vapnik は MP で判別分析にアプローチした。統計や OR の分野を避け、多くの実証研究をパターン認識の分野で汎化能⼒や LSD 判別やカーネル SVM を広めた第三世代の後継者である．
3) 新村は，既存の判別理論には 4 つの問題があることを実証研究で⾒つけ，それらを解決した．その応⽤研究として，統計や医学の多くの研究者が従来の統計⼿法で 10 年から
15 年以上に渡り⾼次元空間の分析を試みたが成果は得られなかった．しかし 3 種類のOLDF は⾼次元空間を少数の遺伝⼦空間に⾃然に縮約し，Matroska 構造をもつこと，
Datasets は少数の排他的な SM の和集合であることを⽰した．すなわち，これらの各 SMは統計的に⼩標本であり，通常の統計⼿法で分析可能である．
4) LASSO 等のアプローチは無意味である．通常のスイス銀⾏や⽇本⾞データで LSD であることが分からず、Feature Selection をできないのに⾼次元でできると考えるのは論理的ではない．

- 28 -

産業連関表に基づく
都道府県クラスターと産業構造推移の可視化
山本けい子（函館高専）
・寒河江雅彦（金沢大学）

概要：産業連関表データの新たな活用法として，地域間の比較や経済構造の推移に着目した分析を行う．平成 17 年
版都道府県単位の産業連関表（取引基本表の内生部門）データに基づいて，都道府県間の経済構造のクラスター分析を行った．また，平成 17 年版と平成 23 年版の同一都道府県における産業連関表データの主成分分析を通して，産業構造の推移の可視化を試みた．これらの結果から，産業連関表を用いた地域間の比較や経年変化の抽出方法について議論する．

1 はじめに

3 分析方法

産業連関表は，国内経済において一定期間（通常 1 年

3.1 データ間）に行われた財・サービスの産業間取引を一つの行列

平成 17 年版産業連関表を各都道府県の Web サイトかに示した統計表で，全国版，地域版，都道府県版など，

らダウンロードし，取引基本表の大分類の内生部門デー

5 年ごとに作成されており，主に，地域内の生産波及効

タを抽出した．ただし，都道府県によって，部門数が若果の計算や経済構造の分析に用いられている．我々は，干異なっており，部門を分割・結合して，全国データに恣意的な選択なく作成過程も統一的な産業連関表データ

記載の 34 部門に統一した．また，平成 23 年版についてを，地域間で比較することにより産業構造の類似性を抽は，大分類が 37 部門へと拡大されているため，分析に出できるのではないかと考え，地域クラスタリングを試使用した都道府県については，平成 17 年版と同じ 34 部

みた．本稿では地域として都道府県単位を設定し，地域門となるようデータを結合した．これらの操作により，間の類似度として，行列ノルムを用いてクラスター分析

都道府県ごとに 34 部門 × 34 部門の行列データを作成した結果について議論する．また，同一地域における産した．

業構造の推移を抽出するため，平成 17 年版と平成 23 年

3.2 方法

版でのデータをそれぞれ主成分分析し，主成分の比較や3.2.1 都道府県クラスターの可視化3.1 で都道府県ごとに作成した 34 × 34 の行列データ

可視化を通して手法の有効性を検証する．

に基づいて，都道府県間の類似性を算出し，都道府県の2 産業連関表

経済構造のクラスター分析を行う．クラスター分析は，

産業連関表は，経済構造を総体的に明らかにし，経済各データを１つのクラスター（まとまり）としてスター波及効果分析 や経済指標の基準改定を行うための基礎

トし，クラスター間の距離（類似度）に基づいて，2 つ資料を提供することを目的として，総務省を中心に各省のクラスターを逐次的に結合していく階層的クラスタリ

庁で 5 年ごとに作成される統計表である．ある産業にングを採用し，樹形図によって可視化した．距離の指標

新たな需要が発生した場合にどのように生産が波及するには行列の要素ごとのノルムか計算することができることから，さまざまな経済分析に用いられている．一定期間において，財・サービスがv
u 34 34 (
)2
u∑ ∑
(A)
(B)
(A)
(B)
X −X
d (X , X ) = t
ij

ij

各産業部門間でどのように生産され，販売されたかにつ

i=1 j=1

いて，行列の形で一覧表にとりまとめられており，取引を用いた．ただし，都道府県 A の行列データを X(A) ，

基本表，投入係数表，逆行列係数表などから成り立って都道府県 B の行列データを X(B) とする．ノルムの他にいる．

距離の指標として，主成分や非負値行列因子分解などの本稿では，都道府県単位で作成されている産業連関表のうち平成 17 年版の大分類から取引基本表の中間投入・利用も考えられる．なお，実装には，R 言語を使用した．

3.2.2 経済構造推移の可視化

中間需要を表す行列（内生部門データ）および平成 23 年

同一都道府県における平成 17 年版データと平成 23

版の大分類から内生部門データを用いた分析を行う．

年版データを比較することにより，経済構造の推移を可視化する．各年の特徴を的確に捉えながら推移の概観を- 29 -

可視化するため，34 × 34 の行列データを主成分分析

データの情報をどの程度表しているかの目安となる数値によって次元縮約し，得られた主成分を可視化することである．

によって推移の読み取りを試みる．主成分は，データの愛知 のH17年̲主成分

4 結果
4.1 都道府県クラスターの可視化

0.0
-0.2
-1.0

鉄鋼
0.0

0.2

0.4

0.6

0.8

輸送機械

-0.4

第2主成分 (0.102 )

0.0
-0.2
-0.4

金属
機械

-0.8

示することによって推移を可視化した．

輸送機械

-0.6

１主成分と第２主成分として抽出し，これらを各年で図第2主成分 ( 0.03 )分，・・・として得られることから，産業構造の特徴を第愛知 のH23年̲主成分
商業
対事業所サ
食料
電力ガス
電機
機械
金属

-0.6

商業
食料
対事業所サ
運輸
電機
建設

-0.8

もつ情報を多く保持するものから第１主成分，第２主成

1.0

鉄鋼
0.0

0.2

0.4

0.6

0.8

第 1主成分 ( 0.91 )第 1主成分 ( 0.8 )東京 のH17年̲主成分

東京 のH23年̲主成分

1.0

取引基本表のデータは，各都道府県内で取引される金クラスター分析を行う．標準化しないデータと標準化し

0.6
0.4

公務
金融
商業 運輸

情報通信
対事業所サ

金融
-0.4

-0.2

第 1主成分 ( 0.69 )都道府県クラスタリング(L2)

愛知

図2

0.0

-0.7

-0.5

-0.3

-0.1

第 1主成分 ( 0.7 )

愛知県と東京都の産業構造の推移

東京

4e+06

図 2 より，愛知県では，第１主成分の寄与率が高く，「輸送機器」と「鉄鋼」の２大産業を有し，平成 17 年と

0e+00

神奈川
埼玉
静岡
北海道
茨城
兵庫
広島
福岡
愛媛
和歌山
大分
福島
長野
宮城
岩手
熊本
岐阜
滋賀
山形
富山
石川
山梨
奈良
香川
沖縄
長崎
青森
宮崎
秋田
福井
徳島
佐賀
島根
鳥取
高知
鹿児島
新潟
京都府
岡山
山口
三重
栃木
群馬
千葉
大阪府

2e+06

Height

6e+06

8e+06

第2主成分 (0.109 )

0.0

不動産

-0.6

たデータを用いた分析結果を図 1 に示す．

不動産

0.0

は，34 × 34 の行列データを標準化したものを使用して

商業
運輸

-0.4

一方，金額よりもその産業の取引構造に関心がある場合

公共サ
対事業所サ
対個人サ

-0.8

行列データそのものを使用したクラスター分析を行う．

第2主成分 ( 0.17 )く異なる．この金額の違いに着目する場合は，34 × 34 の

0.8

0.4

情報通信

0.2

額であることから，都道府県によって金額の単位が大き都道府県クラスタリング(L2̲N)を合わせて 80% の寄与率であり，平成 17 年では「情報

北海道

40

れたが，平成 23 年では，「サービス」系が「情報通信」

30

東京と同程度まで伸びており，さらに「不動産」も一つの大

鳥取
島根
山梨
岐阜
静岡
山形
京都府

宮崎
福井
青森
徳島

香川
熊本
鹿児島
佐賀

富山
三重
兵庫

岩手
福島
愛媛
大分
高知
愛知
広島
秋田
和歌山
奈良
大阪府
新潟
滋賀
埼玉
群馬
長野
千葉
神奈川
茨城
栃木
岡山
山口
福岡

石川

長崎

沖縄

きな特徴として現れていることがわかる．
宮城

Heightない．一方，東京都の場合は，第１主成分と第２主成分
通信」と「金融」，「サービス」系に特徴的な傾向がみら

hclust (*, "complete")

20
10
0

平成 23 年では，産業構造にほとんど変化がみられてい5 考察
本稿では，産業連関表の新たな活用方法として，地域

図 1 都道府県のクラスター分析 (上/下：標準化なし/あり)hclust (*, "complete")

クラスタリングと産業構造の推移の可視化法について提案・検証した．通常，経済データの分析には専門的な知識が必須となるが，ここで扱った手法は，クラスター分

図 1 の標準化しないデータを用いた都道府県クラス
ターの樹形図では，愛知が特出しており，次いで東京，
神奈川／埼玉・静岡のクラスターとその他の都道府県でクラスターが形成されていることがわかる．また，標準化したデータでは，北海道と東京がそれぞれクラスターを形成していることがわかる．

析や主成分分析と可視化を組み合わせたデータマイニン
グ手法であるため，経済に関する事前知識は不要である．
ただし，結果の解釈については，専門家の協力が想定される．また，産業が 34 部門という多く複雑な状況を次
元縮小法によって特徴抽出し可視化することは，容易な解釈と直感的なイメージが得られるという利点もある．

4.2 経済構造推移の可視化

今回使用した産業連関表（取引基本表の内生部門）は

平成 17 年と平成 23 年のデータをそれぞれ主成分分析

各都道府県内のみで取引される金額のため，今後は，都分析し，産業（需要側）の 34 部門についての第１主成

道府県内外との移出・移入も含めた経済構造の分析へ展

分と第２主成分を図示して比較する．図 2 に愛知県と東開が期待できる．

京都の例を示す．軸ラベルのカッコ内は寄与率で，元の- 30 -

テキストマイニングにもとづくレビューのスコアリングを用いた映画の統計的分類
上智大学大学院理工学研究科数学領域・院

1

山縣

一慶

研究目的と背景

今日の日本は, 武士道に由来する武道や茶道などの伝統的様式のみならず, 映画や漫画, アニメなど現代文化が世界的に注目されており, それらは「クール・ジャパン」として総称される. しかし, 統計手法を用いた現代文化の分析は少なく, 分析の需要は今後高まっていくと考えられる. 本研究では, 数ある” クールジャパン” の 1 つである映画に注目する. 中でも, 定義の明確でない映画のジャンル分けに焦点を当て, テキストマイニングを用いた新たな映画のジャンル分けを行う.
テキストマイニングは, 単語や文節で区切り, 各々の出現件数や傾向を解析し, 有用な情報を取り出すことが目的である. 日本語で書かれたテキストを分析する際, 語の最小単位である形態素に分解する必要があり, 一連の工程を形態素解析という. 本研究では形態素解析ソフトの 1 つである TTM [5] を使用した. TTM は同義語を登録（以下, 辞書作り）することで, 同義語の出現件数を分析することも可能である. 下の表は, 辞書作りの一例である. 辞書 1 では, 我々が”SF” という語から連想するであろう語を登録した. 対して辞書 2 は, 同義語辞書（シソーラス）の Weblio [7] を用いて”SF” の同義語を登録した. また, 映画「スター・ウォーズ

フォースの覚醒」に関するレ

ビュー 493 件に対して作成した辞書 2 つと TTM を用い, レビュー中で”SF” の語が使われている件数をレビュー
点数ごとにまとめ, 比較した. 対象レビューは Yhoo!映画 [6] から Ruby 言語で作成したプログラムを用いて抽出
している.
辞書番号

単語

同義語

1

SF

宇宙, ロケット, 船, 星, 彗星, 惑星, 月, NASA2SF

エスエフ, サイエンス, フィクション, 空想科学小説, 近未来小説

辞書１を用いた抽出件数

辞書 2 を用いた抽出件数

レビュー点数

抽出語

出現件数

レビュー点数

抽出語

出現件数

3

SF

23

3

SF

1

4

SF

24

4

SF

4

5

SF

18

5

SF

5

辞書内容の違いが語の出現件数に影響を及ぼすことがわかる.
従来の辞書作りでは, 同義語辞書（シソーラス）を用いることが多い. 短時間で同義語を見つけることができる
が, 分類対象の内容に合った同義語を抽出できないという欠点を持つ. 同義語抽出の先行研究としては, 藤村 [2] が挙げられる. 藤村の研究では, ノート PC に関するレビューをもとに, 肯定的な意味を持つ語と否定的な意味を持
つ語を分別するためのスコアリングを提案している. また, Anne [1] は ”personalization” と”customization” という曖昧な定義を持つ 2 語が分類可能か否か研究している. Anne は 両語について書かれたテキスト 883 件と1544 件を対象に, テキスト内でそれぞれの語と同時に出現する語（以下, 共起語）の出現件数を調べた. そして, テ
キスト内である語と同時に出現する共起語の出現件数の相関から 2 語が分類可能であることを述べている. 本研究では, 映画に関するテキストデータから共起語の出現件数を測定し, クラスター分析を用いて, 互いに関連性の強い語で構成された語群を得た. これらの語群を辞書登録し, 映画分類に用いる辞書 (共起語辞書) を作成した. 映画に関するレビューから共起語辞書で語の抽出を行い, それぞれの特徴をスコアリングすることで新たな映画分類を- 31 -

得るのが本研究の目的である.

2

共起語の出現件数を用いたスコアリング手法

本研究では, IMDB [3] の視聴ランキングに掲載され, レビュー数が充実している映画 32 本を分類対象とした. 対象映画のレビューは,Yhoo!映画から, Ruby 言語を用いて作成したプログラムで抽出した. また, 映画レビューが,映画内容についての文章とレビュー著者の感想について書かれた文章で構成されていると考え, それぞれについて言及される際に使われる単語群をジャンル属性, 感情属性とした. 提案するスコアリングは, 国土交通省 [4] で提案された評価アンケートに対する重みづけの手法に, 対象データの語数による重みを加えたものである.
以下, 共起語辞書を用いたときのスコアリングについて説明する. 映画番号 i について書かれた, 評価点数 j のレ
ビューにおける, ジャンル属性に含まれる単語の総出現件数を eij とし, 感情属性に含まれる単語の総出現件数を

cij とする. また, ジャンル属性, 感情属性に含まれる特徴的な語の番号を k とする. この時, 映画番号 i について書かれた評価点数 j のレビューにおける, ジャンル属性に属する属性番号 k の語の出現件数を cijk , 感情属性に属す
る属性番号 k の語の出現件数を eijk とする. 以上を用いて, スコアリングの定義を行った. 映画番号 i について書かれた評価点数 j のレビューにおけるジャンル属性と感情属性の重みを, それぞれ,

wcij = 100 ∗cij
eij + cijweij = 100 ∗eij
eij + cijとする. また, 映画番号 i における点数 j のレビュー総単語数を sij とする. このとき, 映画番号 i における, ジャ
ンル属性に所属する k 番目の特性が持つスコア Cik と, 感情属性に所属する k 番目の特性が持つスコア Eij を

Cik =

5
∑

(

j=1

cijk

wcij · ∑22

k=1 cijk

sij

)

· ∑5

Eik =

j=1 sij

5
∑
j=1

(

eijk

weij · ∑8

k=1 cijk

sij

)

· ∑5

j=1 sijとする. これらのスコアリングを用いて 32 本の映画全てに対し各属性番号ごとのスコアをつけた.

3

スコアリング結果の一例

下の表はジャンル属性のスコアリングの例である. 比較することでタイトル特有の強い特性の存在が確認できる.

参考文献[1] Anne, S. and Johanna, B., ”Applying text-mining to personalization and customization”, Expert Systems
with Applications, 39, 10049 - 10058 (2012)．
[2] 藤村 滋, 豊田 正史, 喜連川 優, Web からの評判および評価表現抽出に関する-考察”, 社会法人情報処理学会研究報告, (2004)．

[3] IMDB, http://www.imdb.com/chart/top/[4] 国土交通省 公共事業評価システム研究会, ” 公共事業評価の基本的考え方 (公共事業評価システム研究会報告)について”, 国土交通省, (2002).

[5] TinyTextMiner β version , http://mtmr.jp/ttm/[6] Yahoo!映画, http://movies.yahoo.co.jp/[7] Weblio 類語辞典, http://thesaurus.weblio.jp/- 32 -

規模の異なる変量群をもつロジスティック回帰モデルの係数 2 段階推定
上智大学大学院理工学研究科数学領域・院

1

竹下

佳宏

• SNP ごとにすべての回帰係数を推定し直すが，研究の背景

本来は環境要因は SNP から独立しているため，

21 世紀に DNA の全塩基配列情報を把握することが可能となり，ゲノムワイド関連解析（以下，GWAS）が行われるようになった [1]．GWAS では

環境要因の回帰係数が j 番目の SNP に依存していることが不自然である• SNP の数だけ環境要因の回帰係数を推定するのは計算効率が悪い

一塩基多型（以下，SNP）と呼ばれる数十万単位ある
特定の塩基配列情報に着目して，生物の性質をひも親しまれているが，機械学習の技術が応用されるな• １つのモデルにつき１つの SNP しか考慮しな
いため，複数の SNP の交互作用が説明できない

ど，最近では新しい技術も活用されている [2]．本研

これらの問題を解決する回帰係数の推定を本研究で究では，遺伝子情報の解析に用いられるロジスティッ

行った．

とく．遺伝情報の解析手法では，古典的な統計解析が

ク回帰モデルの新たな係数推定方法を提案する．

2

従来の研究手法
ある i 番目の個体に対し，病気である (yi = 1) かそうではない (yi = 0) かという 2 値の応答変数を，
遺伝子の情報や，その個体の環境情報を考慮して説明するモデルとして，ロジスティック回帰モデルがある．一般には病気の有無以外の 2 値の応答変数にもこのモデルは用いられるが，本研究では病気の有無についてのみ言及する．
me
me 個の環境要因 {Ei,k }k=1と ms 個の遺伝要因でms
ある SNP 情報 {SN Pi,j }j=1のうち j 番目の SNP のみを用いて，個体 i が病気にかかる確率 pi,j を次のロジスティック回帰モデルで表現する [2]．
(
)
pi,j
log
1 − pi,jme∑= α0,j +αk,j Ei,k + βsnp,j SN Pi,j3回帰係数の 2 段階推定

3.1

推定方法の概要

ロジスティック回帰モデルは，応答変数が 0 か１
かの 2 値の他に，応答変数が生起または反応の確率である場合にも利用できる [3]．現実には生起の確率そのものが応答変数として観測されることはほとん
どないため，標本全体に対する生起数の割合で応答変数の値を推定する．
この考え方を利用して，2 段階で回帰係数の推定を行う．

1. 標本の個体 i ∈ {1, 2, . . . , n} について，i と同じ環境下で罹患している標本数を求め，標本全体に対する割合で罹患確率 pbe,i を推定

2. 次のような環境要因だけのロジスティック回帰
モデルを考え，係数 α = (α0 , . . . , αme ) を推定
(1)

k=1

(
log

SNP 情報は数十万単位あるため，一度に SNP をす
べて解析するのではなく，１つの SNP 情報に対してその都度モデルを構築し、その中から病気を十分に説明できるようなモデルを考えることで，病気に関pbe,i1 − pbe,i)= α0 +me∑

αk Ei,kb = (b
3. 回帰係数の推定値 ααk )0≤k≤me を用いて，個体 i の罹患確率の推定値 pee,i を次の式で計算回帰係数に推定値を代入して pi,j について解き，そb ′i )
exp(αEb ′i )
1 + exp(αEの値があるしきい値を上回れば病気 (y = 1)，下回

ここで，E i = (1, Ei,1 , . . . , Ei,me )連を持つ SNP を見つける．実際はモデル式 (1) のれば健康 (y = 0) と予測する．
従来の研究手法においては，以下の点で問題があった．

- 33 -(2)k=1

pee,i =

4. pbe,i と pee,i の差，すなわち，環境要因だけでは説明できない罹患確率の部分を SNP 情報を用

いて説明．回帰係数の問題としては，次のモデ

たがって最小自乗解を得る．以下，環境要因で説明

ル式で回帰係数 β = (β0 , . . . , βms ) を推定
(
)
(
)
pbe,i
pee,i
log
− log
1 − pbe,i1 − pee,ims∑= β0 +
βj SN Pi,j
(3)

し残した部分を次のように定める．

j=1

また，SNP 情報を n × (ms + 1) 行列 S として表す

2 段階に分けることにより，環境要因の回帰係数は
SNP の影響を受けずに推定できる．そのため, (2) の回帰係数は (1) 式とは異なり, 添字 j が消去される．
3.2
3.2.1

′

r = (r1 , . . . , rn ){}{}pbe,i (δ)
pbe,iri = log− log
1 − pbe,i (δ)
1 − pbe,iとき，最小自乗解は次のように表される．

()−1 ′β̂ = S ′ S
Sr(5)

ただし，行列 S の 1 列目は切片項に対応するため，

環境要因の回帰係数の推定成分が 1 だけの列とする．

応答変数の補正

ところが，一般に SNP の情報はサンプル数 n よ

現実のデータでは，pbe,i = 0, 1 もあり得る．そのりも SNP 数 ms の方が多いため，S にはランク落

場合，(2) と (3) が定義できない．そこで，十分小ちが生じる．そこで，(5) ではなく，β は S のムー

さい δ > 0 を用いて，pbe,i を次のように補正する．

ア・ペンローズ型一般逆行列 S + を用いた次の形で δ
pbe,i (δ) =
pbe,i


1−δ

推定する．

: pbe,i = 0: 0 < pbe,i < 1: pbe,i = 1

b = S+rβ+(4)

推定量 β̂ + は次の性質を持つ [4]．

pbe,i を pbe,i (δ) で置き換えて，3.2.2 節以降で説明する回帰係数の推定を行う．
3.2.2

b = arg min ∥r − Sβ∥
1. β
+
β

2. β̂+ は不偏推定量

最尤法による推定3. 任意の不偏推定量 β̃ に対して，
( )
( )
b
e
V β
+ ≤V β環境要因の回帰係数 α は，{bpe,i (δ)}ni=1 を応答変数の値として用いた最尤法で推定する．実際には最尤解は陽な形で得ることができないので，ニュートン法によって数値的に求める．

b δ は，補正量 δ > 0{bpe,i (δ)} にもとづく最尤解 αの影響を受けており，本来の {bpe,i }ni=1 にもとづく解b とはずれがある．けれども，次の命題から，αbδ がαb の近似解となっていることがわかる．
αb を補正がない場合 (δ = 0) の {b命題 αpe,i }ni=1 による最尤推定量とするとき，{}b δ − α)b 2 =0lim E (αここで，V (・) は分散共分散行列で，不等号は 2
次形式に関する大小関係を表す．

4

数値実験

上記の推定方法を実際のデータに適用し，従来の推定方法との比較を行った．数値実験の結果もシンポジウム当日に発表した．

参考文献

δ→+0

b は漸近的よく知られているように，最尤推定量 αに正規分布に従う．命題の結果から，十分小さな δ にb δ も近似的に正規分布に従う．よって，αbδ
ついて，αについて正規性にもとづく有意性の検定ができる．

[1] 鎌谷直之．遺伝統計学入門, 岩波書店, 2007.
[2] 上辻茂男．”ゲノムワイド関連研究に学ぶ遺伝統計学”, 計算機統計学，25 巻，第 1 号, 17-39 (2012).
[3] 中村永友．多次元データ解析法 (R で学ぶデータ
サイエンス 2), 共立出版, 2010．

3.3

遺伝要因の回帰係数の推定

3.1 節の (3) 式にもとづいて，SNP の回帰係数 β を推定する．回帰モデルとしては，重回帰モデルにし- 34 -[4] 柳井晴夫，竹内啓．射影行列・一般逆行列・特異値分解, 東京大学出版, 1983．

国勢調査の匿名化ミクロデータの作成可能性について―地域区分に着目して―*中央大学経済学部 伊藤 伸介**
(独)統計センター 星野 なおみ
阿久津 文香
総務省統計局
1．はじめにわが国では、平成 25 年 12 月と平成 26 年 3 月にそれぞれ、平成 12 年と 17 年の国勢調査の匿名データの提供が開始されている。一方、わが国の国勢調査の提供済匿名データでは、都道府
県及び人口 50 万以上の市区が最小の地域区分となっているが、より詳細な地域区分を用いた分析が可能なミクロデータに対するニーズが存在すると思われることから、国勢調査の匿名デー
タに関しては、より詳細な地域区分の提供可能性が検討されてよいと考える。そこで、本稿では、秘匿性と有用性の両面から、詳細な地域区分の作成可能性を検討する。

2. 「地域の人口規模の閾値」に基づく秘匿性の検証
匿名化ミクロデータ(個票データに匿名化処理が施されたデータ)の作成・提供においては、秘匿性に関する閾値を設定した上で、その閾値を超えない形で、様々な匿名化技法が適用されることが考えられる。例えば、Hawala(2001)は、アメリカ人口センサスの Public Use Microdata
Sample の作成において用いられる 10 万人という地域区分の閾値に関して、その秘匿性に関する事後検証を行うために、母集団一意(population unique)の比率を用いて地域の人口規模と秘匿性の指標との関連性を明らかにしている。本研究では、平成 22 年国勢調査の A 県の調査票情報(個票データ)を用いて、地域の閾値を変更した場合のキー変数におけるリコーディング(区分統合)の可能性を探った。
最初に、秘匿性に関する第 1 の研究として、伊藤・星野(2014)に基づき、リコーディングを行ったデータに対して、①住宅の建て方、②住居の種類、③性別、④配偶者の有無、⑤国籍、⑥労働力状態、⑦従業上の地位、⑧年齢、⑨産業と⑩職業の 10 個のキー変数をもとに母集団一意(population unique)の比率を計測した。具体的には、10 変数の中で住宅の建て方等の 7 変数については、提供済匿名データ(現在提供されている国勢調査の匿名データ)の区分が用いられているが、残りの 3 変数である年齢、産業と職業については、様々なリコーディングおよびトップコーディングを施した上で、全 81 パターンにおける母集団一意の比率の計測を行った。なお、リ
コーディングにおいては、いわゆる「0.5%基準(単変量において母集団の 0.5%を下回る区分を統合すること)」を考慮した区分も設定されている。
本研究では、地域の人口規模と属性の分類区分との関連性を実証的に明らかにするために、A県を対象にし、地域の閾値については、(1)人口 20 万人以上地域(県庁所在市に該当する地域も含む)、(2)人口 10 万人以上地域、(3)人口 5 万人以上地域、(4)人口 3 万人以上地域、(5)人口 2 万人以上地域、(6)人口 1 万人以上地域、(7)人口 5000 人以上地域、および(8)人口 1000 人以上地域の 8 パターンを設定し、それに該当する 20 地域を選定した。
人口 20 万人以上地域を対象にした場合の母集団一意の比率の算出結果によれば、年齢、産業と職業が原区分である場合の母集団一意の比率は約 18.43%であるが、年齢のみを 5 歳区分に統合した場合、比率が 8.83％に減少することがわかった。その一方で、産業や職業についてのみ
匿名データと同様の区分でリコーディングを行っても、母集団一意の比率はそれぞれ 17.47％と17.64％であり、大きな変化は見られないことが確認された。
一方、本研究では、地域の人口規模と母集団一意の比率との間にトレードオフの関係がみられることが実証的に確認された。このことは、年齢、職業、産業と分類区分が設定された場合に、母集団一意の比率に関する閾値を適切に定めることができれば、秘匿の観点から、提供可能な地
域の人口規模の閾値を導出することが可能なことを意味している。
つぎに、秘匿性に関する第 2 の研究として、年齢、産業、職業と地域の人口規模に秘匿処理を施した場合の母集団一意の比率に対する影響に関する重回帰分析を行った。具体的には、母集団一意の比率を被説明変数とし、説明変数に関しては、年齢のリコーディングおよびトップコーデ
ィングに関する 9 パターン、産業のリコーディングに関する 3 パターン、職業のリコーディン* 本稿は、伊藤・星野・阿久津(2016b)に基づいている。なお、本稿の内容は個人的な見解を示すものであり、統計センターの見解を表すものではないことに留意されたい。
** (独)統計センター非常勤研究員- 35 -

グに関する 3 パターンのそれぞれに関するダミー変数(パターンに該当する場合には 1、そうでない場合には 0)、および対数変換された地域の人口規模をモデルに設定した。
本分析結果によれば、地域の統合が、年齢、産業、職業と比較しても母集団一意の比率に大きな影響を及ぼすことが確認された。このことから、秘匿性の観点から見た場合、匿名データの作成において、地域区分の設定が相対的に重要な要素を占めることが定量的にも明らかになっている。

3. 情報量損失に基づいた有用性の定量的な評価
つぎに、本研究では、様々な地域を対象に有用性の検証を行った。ミクロデータにおける有用性の定量的な評価方法については、クラーメルの V といった関連性の指標の算出や原データからの絶対距離の平均値(average absolute distance)の計測等を行うことが考えられる(伊藤・星野(2014))。また、情報量損失に関する指標の 1 つであるエントロピーを用いて、秘匿の観点から許容可能な分類区分の組み合わせに関する情報量損失の計測を行うことも可能である(伊藤ほか(2016a))。一方、本研究では、原区分からリコーディングを行った場合の距離を計測することによって、情報量損失の計測を行う。そのために、原区分と統合区分におけるクロス表の差に関する指標を作成し、原区分と統合区分におけるクロス表の差の検証を行った。具体的には、区分統合を行った場合、リコーディング後の度数をリコーディングの対象となった区分数で除することによって、度数の按分を行う。つぎに、リコーディング前のクロス表とリコーディング後に按分済みの度数が入力された表を用いて、情報量損失を算出した。
最初に、分析対象である 20 地域を対象に、81 の変数のパターンにおいて算出される情報量損
失値に関する基本統計量を算出した。人口規模が小さい地域については、情報量損失値の平均値が相対的に小さくなっているだけでなく、標準偏差も小さいことが確認できる。このことから、地域の人口規模が小さい場合には、区分統合しても情報量損失の変化が小さいことが推察される。
つぎに、情報量損失と地域規模との関連性を確認したところ、地域の人口規模が大きくなるほど、情報量損失が大きくなる傾向にあることが実証的に明らかになった。その一方で、地域の人口規模が相対的に小さい場合、情報量損失における大きな違いが見られないことを確認することができる。

4. むすびにかえて
本稿では、国勢調査の匿名化ミクロデータを用いて、地域の人口規模に基づく閾値を設定した
場合の母集団一意の比率と情報量損失の検証を行った。本研究の結果を踏まえると、個人単位で抽出した匿名化ミクロデータにおいては、地域の人口規模と母集団一意の比率との間にトレードオフの関係があることを実証的に明らかにすることができた。さらに、情報量損失に関する分析結果からは、地域の人口規模が大きいほど、情報量損失が大きくなること、秘匿処理が情報量損失に及ぼす影響は、属性によって異なることがわかった。
一方、本分析結果によれば、年齢、産業、職業を含むキー変数を用いて算出された母集団一意の比率をもとに、秘匿性の閾値を適切に設定することができれば、提供済匿名データにおいて定められている地域の人口規模 50 万以上という区分より細かな地域区分の設定についても、議論することが可能であると思われる。
こうした議論に基づいて、秘匿性と有用性の両面からの実証研究を踏まえた上で、複数ファイ
ルにおける地域詳細化データや年齢等の属性に関する詳細化データの作成可能性を模索する必要があるように思われる。

参考文献
Hawala, S.(2001) “Enhancing the "100,000 rule" On the Variation of the Per Cent of Uniques in A
Microdata Sample and the Geographic Area Size Identified on the File”, Proceedings of the Annual
Meeting of the American Statistical Association.
伊藤伸介・星野なおみ(2014)「国勢調査ミクロデータを用いたスワッピングの有効性の検証」『統計学』107号，2014 年 9 月 30 日，1～16 頁
伊藤伸介・星野なおみ・阿久津文香(2016a)「国勢調査における匿名化ミクロデータの有用性と秘匿性の定
量的な評価」『製表技術参考資料』No.32, 1～33 頁
伊藤伸介・星野なおみ・阿久津文香(2016b)｢国勢調査ミクロデータに対する匿名化措置の可能性に関する研究｣『製表技術参考資料』No.34，1～59 頁- 36 -

Estimation of the shape of density level sets ofstar-shaped distributions
Hidehiko Kamiya∗
Faculty of Economics, Osaka University of Economics
December 2016
Elliptically contoured distributions generalize the multivariate normal distributions
in such a way that the density generators need not be exponential. In this way, theclass of elliptically contoured distributions includes, for example, distributions whose tails
are heavier than those of the multivariate normal distributions. However, as the namesuggests, elliptically contoured distributions remain to be restricted in that the similar
density contours need to be elliptical. Hence, skewed distributions are not included inthis class.
Kamiya, Takemura and Kuriki [5] proposed star-shaped distributions in which thedensity level sets are allowed to be arbitrary similar star-shaped sets (see also [6], [4]).
Essentially the same idea can be found in v-spherical distributions by Fernández, Osiewalski and Steel [2] and center-similar distributions by Yang and Kotz [7]. Asymmetry isallowed in star-shaped distributions. Hence, besides distributions which are symmetricwith respect to the center such as elliptically contoured distributions and lq -spherical distributions, the class of star-shaped distributions includes asymmetric distributions suchas multivariate skewed exponential power distributions.
Kamiya, Takemura and Kuriki [5] studied distributional properties of star-shaped distributions, e.g., the independence of the “length” and the “direction,” and the robustnessof the distribution of the “direction.” However, they did not investigate inferential problems about star-shaped distributions. From the perspective of [5], the most important
problem in the inference for star-shaped distributions is the estimation of the shape of
their density contours
In this paper, we propose a nonparametric estimator of the shape of the density contours. The point is that the density of the usual direction under a star-shaped distribution
is in one-to-one correspondence with a function which determines the shape of the density contours. Thus, by nonparametrically estimating the density of the direction, we
can obtain a nonparametric estimator of the shape. We prove its strong consistency withrespect to the Hausdorﬀ distance.
Our main result is the following:∗This work was partially supported by JSPS KAKENHI Grant Number 25400201.

- 37 -

sample from a starTheorem 0.1. Let x1 , . . . , xn ∈ X = Rp \ {0}, p ≥ 2, be an∑i.i.d.
np−1
ˆshaped distribution h(r(x))dx. Let fn (u) = (C(η)/(nη )) i=1 L((1 − uT ui )/η 2 ) be akernel estimator (Hall, Watson and Cabrera [3], Bai, Rao and Zhao [1]) of the density
f (u) of u = x/∥x∥ ∈ Sp−1 , x ∼ h(r(x))dx, based on ui = xi /∥xi ∥, i = 1, . . . , n.
Assume the equivariant function r : X →∫ R>0 under the action of the positive real−pnumbers is continuous andso that Sp−1 r(u)du = 1, and that L : R≥0 → R≥0∫ ∞normalized∫∞(p−3)/2
is bounded and satisﬁes 0 L(v)v
dv > 0 and 0 supw: |√w−√v|<1 L(w) · v (p−3)/2 dv <∞. Moreover, suppose η = ηn > 0 is taken in such a way that ηn → 0 and nηnp−1 / log n →∞ as n → ∞.
Then, Ẑn = {fˆn (u)1/p u : u ∈ Sp−1 } is a strongly consistent estimator of the shape
Z = {x ∈ X : r(x) = 1} of the density contours of the star-shaped distribution in thesense that the Hausdorﬀ distance δH (Ẑn , Z) between Ẑn and Z satisﬁesδH (Ẑn , Z) → 0
In addition, Ẑn =δH (Ẑn , Z) → 0 a.s.

∪0≤c≤1 cẐn(n → ∞)a.s.

is a strongly consistent estimator of Z =∪0≤c≤1 cZ:

References[1] Z. D. Bai, C. R. Rao and L. C. Zhao, Kernel estimators of density function of
directional data, Journal of Multivariate Analysis 27 (1988), 24–39.
[2] C. Fernández, J. Osiewalski and M. F. J. Steel, Modeling and inference with vspherical distributions, Journal of the American Statistical Association 90 (1995),
1331–1340.
[3] P. Hall, G. S. Watson and J. Cabrera, Kernel density estimation with spherical
data, Biometrika 74 (1987), 751–762.
[4] H. Kamiya and A. Takemura, Hierarchical orbital decompositions and extended
decomposable distributions, Journal of Multivariate Analysis 99 (2008), 339–357.
[5] H. Kamiya, A. Takemura and S. Kuriki, Star-shaped distributions and their generalizations, Journal of Statistical Planning and Inference 138 (2008), 3429–3447.
[6] A. Takemura and S. Kuriki, Theory of cross sectionally contoured distributions andits applications, Discussion Paper Series 96-F-15, July 1996, Faculty of Economics,
University of Tokyo.
[7] Z. Yang and S. Kotz, Center-similar distributions with applications in multivariate
analysis, Statistics & Probability Letters 64 (2003), 335–345.

- 38 -

正方分割表の一致率検定のための代数的方法
吉田 知行 (北星学園大学 経済学部)

t-yoshida@hokusei.ac.jp
概要: 2 元以上の正方分割表の一致率の検定について，正確な p-値を求める方法を紹介する．この方法は，フィッシャーの並べ替え検定，フィッシャーの正確確率法の流れをくむ方
法で，代数的には，対称群のデータセットへの作用と指標理論の応用分野である．また，有限群上のランダムウォーク (RW) とも関係している．
キーワード： 正確確率法，並べ替え検定，有限群上のランダムウォーク，分割表の列挙問
題，有限群の表現．

一致率検定の正確な p-値を求める方法とその背景にある代数理論を紹介した．主結果は次である．正方分割表 X = (xij ) の周辺和を a = (xi+ ), b = (x+j ) とする．周辺分布が a, b であるような正方分割表の集合を tab(a, b) とする．P (k) を，一致数が k 以上になる確率とする:P (k) = Prob(x ∈ tab(a, b) | Tr(x) ≥ k).
数列 {q(k)}, {p(k)} (k = 0, 1, · · · , n) を次で定義する：∏

2 F0 (−ai , −bi ; z) =

λ

∑ (n)

k≥0

p(k) =

k

( )
j−k j
(−1)
q(j),
k

k!q(k)z k ,

n
∑

n
∑

j=k

k=0

k

p(k)z =

n∑q(k)(z − 1)kk=0

0 ここで 2 F0 は超幾何多項式である．このとき P (k) は次で与えられる．
n
∑
k=0

k

P (k)z = 1 + zn∑q(k)(z − 1)k−1

k=1

この結果は，対称群の直積 Sn × Sn のデータセット DS(a, b への可移な作用と，対称群上のランダムウォークから来ている:
Sn 上の RW −→ Sn × Sn /Sn∆ 上の RW −→ DS(a, b) 上の RW −→ tab(a, b) 上の RW- 39 -

本来なら，分割表をランダムサンプリングして，その中から Tr(x) ≥ k となるものの割合を求める．しかし一致数の場合は超幾何多項式で表せたのである．一致数の p-値は，比較言語学，とくにポリアによる二項検定法やオズワルトのシフトに現れる．後者は安本美典による日本語の起源探求に使われた.

参 考 文 献[吉田 07] 吉田知行，分割表の一致率検定への有限群論と組合せ論の応用，代数学シンポ
ジウム，於神戸大学．ネットに pdf ファイルあり．

[吉田 09] 吉田知行，Finite Gelfand pairs and Markov chain Monte-Carlo method, 「表現論と組合せ論」RIMS 講究録第 1689 (201),164-170,- 40 -