令和４年度科学研究費 基盤研究(A) (課題番号：20H00576)「⼤規模複雑データの理論と⽅法論の⾰新的展開」(研究代表者：⻘嶋誠(筑波⼤学))によるシンポジウム

多様な分野における統計科学の理論とその応⽤

シンポジウム報告書

令和４年１０⽉２７⽇（⽊） 〜 １０⽉２９⽇（⼟）
研究分担者：⽥畑耕治 (東京理科⼤学)
開催責任者：⽥畑耕治, 安藤宗司, ⽯井晶, 中川智之 (東京理科⼤学)
場所：
東京理科大学 野田キャンパス 7 号館 6 階 講堂とのハイブリッド（Zoom）
内容・目的：
多種多様な分野において現れる大規模複雑データの解析において，統計科学への期待と関心は高まる一方である．本シンポジウムでは，カテゴリカルデータ解析，医療統計学，高次元データ解析，ロバスト統計，ベイズ統計などの新しい理論や方法論に関する研究，実データを用いた事例研究，データサイエンス時代の統計教育法など幅広い分野からの講演を募集します．多様な分野の研究者が知識や意見を交換することで，数理データサイエン
スにおける最新動向や問題点を共有し，問題解決の糸口を創出する場になることを目的と
します．

共催：
東京理科⼤学,東京理科⼤学データサイエンスセンター
東京理科⼤学 統計科学研究部⾨

科研費シンポジウム「多様な分野における統計科学の理論とその応用」
日時：2022 年 10 月 27 日（木）〜 10 月 29 日（土）場所：東京理科大学野田キャンパス

7 号館 6 階講堂とのハイブリッド（Zoom）

12:30–13:00 開場
13:00–13:10 開会
13:10–15:20 セッション 1

10 月 27 日（木）◇ 「NDB サンプリングデータセットを用いた精神神経疾患の合併症の網羅的分析」
石井一夫（公立諏訪東京理科大学, 久留米大学）
吉永泰周（福岡歯科大学）
坂上竜資（福岡歯科大学）
小路純央（久留米大学）森川渚（久留米大学）野原夢（久留米大学）
野原正一郎 （久留米大学）
福本義弘（久留米大学）◇ 「連検定による重心動揺データの解析」
竹内直子（大阪府立大学）
綿森葉子（大阪公立大学）
久利彩子（大阪河崎リハビリテーション大学）
有末伊織（関西福祉科学大学）◇ 「GMANOVA モデルにおける新たな経時変動の推定方法とその解釈」
永井勇（中京大学）

15:20–15:30 休憩
15:30–17:40 セッション 2◇ 「統計科学におけるエントロピー概念の誤り訂正」
得丸久文◇ 「多様体学習による銀河進化の探求」
竹内努（名古屋大学, 統計数理研究所）
Suchetha COORAY （名古屋大学）◇ 「Fréchet 距離を用いた分類問題について」
伊森晋平（広島大学）
若木宏文（広島大学）

9:30–10:00 開場
10:00–12:10 セッション 1

10 月 28 日（金）◇ 「正方分割表における松下距離に基づく周辺同等性からの隔たり尺度」
篠田覚(横浜市立大学)
吉本拓矢（中外製薬株式会社）
田畑耕治（東京理科大学）◇ 「方向データのためのカーネル密度推定量のバイアス修正」
鶴田靖人（長野県立大学）◇ 「Outcome-adaptive lasso と adaptive lasso を用いた AIPW 推定量の性能評価」
本江渡（東京理科大学大学院, ノバルティスファーマ株式会社）
安藤宗司（東京理科大学）
土田潤（同志社大学）
寒水孝司（東京理科大学）

12:10–13:40
13:40–15:50-お昼休憩セッション 2◇ 「特異ウィシャート行列の固有値の正確分布論」
清水康希（東京理科大学）
橋口博樹（東京理科大学）◇ 「ロジスティック回帰モデルにおける安定な非凸スパース正則化法」
塘由惟(東京大学, 国立精神・神経医療研究センター)
小川光紀(東京大学)
片井みゆき (政策研究大学院大学医学, 東京女子医科大学)
大庭幸治(東京大学)
松山裕(東京大学)◇ 「ブリッジ推定量を用いた BIC の妥当性とその周辺」
宮田庸一（高崎経済大学）

15:50–16:00
16:00–18:10-休憩セッション 3◇ 「ロジスティック分布における母数推定について」
作村建紀（法政大学）
柳本武美（統計数理研究所）◇ 「対数オッズ比の推定量と検定を改善するためのベイズ法」
小椋透（三重大学）,
柳本武美(統計数理研究所)◇ 「共役解析の再構成と拡張の試み」
柳本武美（統計数理研究所）

9:00–9:30 開場
9:30–11:40 セッション 1

10 月 29 日（土）◇ 「大規模時空間データに対するベイズモデル」
若山 智哉 （東京大学）◇ 「観察研究の効果推定値を標的集団に一般化／移送する方法の検討」
堀江悠生（東京理科大学大学院）
篠崎智大（東京理科大学）◇ 「トーリックモデルからの直接抽出の代数的アルゴリズム」
間野修平(統計数理研究所)
高山信毅（神戸大学）

11:40–11:50

閉会

キャンパス情報＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
HP:https://sites.google.com/view/sympo221027/
アクセス:https://www.tus.ac.jp/info/access/nodcamp.html
キャンパスマップ:https://www.tus.ac.jp/info/campus/noda.html＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝NDB サンプリングデータセットを用いた精神神経疾患の合併症の網羅的分析
石井 一夫 1,4,吉永 泰周 2,坂上 竜資 2,小路 純央 3,森川 渚 4,野原 夢 4,野原 正一郎 4,福本 義弘 41,4

公立諏訪東京理科大学工学部,2 福岡歯科大学口腔歯学部, 3 久留米大学高次脳疾患研究所/久留米大学医学部神経精神医学講座, 4 久留米大学医学部心臓・血管内科
日本では少子高齢化が進行し、要介護者の増加による医療・介護システムの崩壊が危惧されている。我々は、生活習慣病と精神疾患および歯科疾患が生活習慣病や介護の重症化に密接に関連していることに着目し、医療ビッグデータを用いた生活習慣病患者のメンタルケアおよびデンタルケアの研究を展開している。今回、NDB サンプリングデータセットを用いて「認知症」、「うつ病などの気分障害」、「ストレスなどの不安障害」、「睡眠障害」、「依存症」などの精神疾患と生活習慣病との関連を網羅的に調査したので報告した。このうち「認知症」と「睡眠障害」について以下にまとめた。
1. 認知症と生活習慣病との関連性調査（1）概要
近年、日本では少子高齢化が問題となっており、高齢化率と要介護者の急増から、医療・介護の崩壊が起こると危惧されている。要介護の要因には認知症、脳血管疾患、高齢による衰弱等があるが、中でも認知症は全体の 18.7%を占めている。認知症は主に「アルツハイマー型認知症」、「脳血管性認知症」、「レビー小体型認知症」、「前頭側頭型認知症」の 4 種に分類され、アルツハイマー型認知症が 67.6%、脳血管性認知症が 19.5%と大半を占めている。
調査にあたり利用可能なデータベースの 1 つに、厚生労働省が「高齢者の医療の確保に関する法律」に基づき、2009（平成 21）年よりレセプト情報並びに特定健診・特定保健指導情報を収集した「レセプト情報・特定健診等情報データベース(National Databaseof Health Insurance Claims and Specific Health Checkups of Japan: 以下、NDB)」があり、NDB は国民皆保険制度下にある日本においては国民の医療の実態を全数に近い割合で評価できることから、非常に貴重なデータであり、幅広い分野と多くの産業に活用されることが期待されている。本研究では認知症と生活習慣病の関連性調査を目的に、NDB から事前に抽出された NDB サンプリングデータセットに基づいて調査を行った。
（2）方法
対象のデータより、ICD-10 に基づき、認知症と生活習慣病の有病者を特定し、疾患の有無と疾患の関連について効果判定を行った。疾患は、傷病名コード、診療行為コード、医薬品コードに基づいて分類した。認知症は「アルツハイマー型認知症」、「血管性認知症」
を、生活習慣病は、「標準的な健診・保健指導に関するプログラム（確定版）」に基づき、

「がん」
、
「糖尿病」
、
「高血圧」
、
「脂質異常症(高脂血症)」、「脳血管疾患」、「虚血性心疾患」、「高尿酸血症」、「肝機能障害」、「動脈閉塞」「血圧性腎臓障害」を用いた。疾患の併存状況は、オッズ比および調整済オッズ比、カイ二条検定に基づいて評価した。
（3）結果および考察
アルツハイマー型認知症と血管性認知症共に多くの生活習慣病に分類される疾患においてオッズ比が 1 を超えていることから、アルツハイマー型認知症と血管性認知症での傾向の多少の違いはあるが、どちらも生活習慣病との関連性が示唆された。がんは入院においてほぼ関連性が認められなかった。また、他の月についても同様の傾向が見られた。演者
らは、NDB オープンデータを用いて都道府県レベルで、精神神経疾患と生活習慣病の関連を示してきた。今回、NDB サンプリングデータセットを用いて、個人レベルで、認知症と生活習慣病の合併を示すことができた。
2.睡眠障害と生活習慣病との関連性調査(1) 概要睡眠障害における生活習慣の乱れは、近年において大きな社会問題となっている。日本は、特に世界的にも睡眠時間が少ないとされている。また、考えられる
様々な要因や背景によって、さらに質のよい睡眠が損なわれている可能性がある。睡眠障害の一つである不眠症や睡眠時無呼吸症候群は、生活習慣病の要因になる因子であると示唆されている。我々は、レセプト情報・特定健診等情報デー
タベース(NDB)サンプリングデータセットを用いて、精神神経疾患と生活習慣病の関連を調べる網羅的な研究を開始した。今回、サンプリングデータセットの医科
入院外、医科入院、DPC、調剤の４つのレセプトを用いて、睡眠障害の実態を調査し、さらには睡眠障害と生活習慣病の関連性を調査した。
(2)方法及び結果
NDB サンプリングデータセット（2011〜2019 年 1 月、4 月、７月、10 月及び、2020 年 1月）を用いたデータ分析を行った。医科入院外、医科入院、DPC の各レセプトの傷病名コードに対し ICD-10 コードに基づき、各疾患（生活習慣病、及び精神疾患）を集計した。睡眠障害（G47）の有無について、各疾患の罹患率を特定し、オッズ比を算出して評価した。その結果、睡眠障害がある人は様々な生活習慣病（糖尿病、高血圧、高脂血症、高尿酸血症、虚血性心疾患、肝機能障害など）を併発していることや、アルコール依存症やパーキンソン病、双極性障害、統合失調症、その他の不安障害といった精神疾患を併発していることが有意に多いことが示された。
謝辞
本研究は JSPS 科研費 22K10587 の助成を受けたものである。

連検定による重⼼動揺データの解析
⽵内直⼦ 1,3

綿森葉⼦ 2

久利彩⼦ 3

有末伊織 41 ⼤阪府⽴⼤学 2 ⼤阪公⽴⼤学 3 ⼤阪河﨑リハビリテーション⼤学 4 関⻄福祉科学⼤学

はじめに
連検定は、データ列のランダム性を調べる検定法であり 1)、様々な解析ソフトで関数が⽤意されているため容易に活⽤することができる。連検定をヒト⽴位重⼼動揺データへ適⽤させることで⽴位制御の規則性を調べた。その具体的な適⽤⽅法と結果、および有⽤性を紹介した。⽴位の安定度によって加速度の向き変化に違いがあることを捉え、動揺制御の質的評価の可能性が⽰唆された。

背景とデータ
健康寿命の延伸に転倒予防は重要であり、⾃⾝の⽴位バランス能⼒を把握することが有効な対策になる。転倒リスクが特に⾼くなる動作変換時のバランス制御を質的に評価するため、⼒の代替として
加速度に着⽬した。制御⽅法を時系列として捉えるために連検定を⾏った。
両脚⽴ちの姿勢から⼀歩踏み出し、⽚脚⽴位保持の後、両脚⽴位保持まで連続して⾜圧中⼼（Centerof Pressure：COP）の計測を⾏った。踏み出し動作は 4 パターン（前⽅ステップ、前⽅昇段、側⽅ス
テップ、側⽅昇段）で、それぞれ 12 試⾏ずつ計測した。踏み出した後、⽚脚⽴位保持を約 10 秒間、
その後両脚⽴位となり約 10 秒間、合計 20 秒間のデータを取得した。サンプリング周波数は 100Hz。
被験者は健常成⼈ 3 名。重⼼動揺計に荷重がかかり始めてから 7 秒までの各 1 秒間を解析区間 T1〜
T7 とし、計測最後の１秒間を解析区間 Tb とした。

解析
COP 加速度の向きを算出し、その⾓度データの偏差Δa により以下の 2 種類の 2 値化を⾏った。
① 「同」：−90°＜Δa＜90°、「逆」：−180°≦Δa<−90°または 90°＜Δa≦180°② 「右」：−180°＜Δa＜0°、「左」：0°＜Δa＜180°
2 種類の 2 値化で得られたデータ列それぞれに対し、解析区間ごとに、帰無仮説「ランダム」、対⽴
仮説「混合不⾜」の⽚側検定、および、帰無仮説「ランダム」、対⽴仮説「混合過剰」の⽚側検定の連検定を⾏った（p＜0.05）。検定には統計解析ソフト R の関数 tseries::runs.test を⽤いた。

結果
2 値化①の場合の結果を図 1 に⽰す。横軸は各解析区間、縦軸は 12 試⾏中「混合不⾜」「混合過剰」「ランダム」となった頻度を表す。「混合過剰」は重⼼の位置が留まっている時の特徴であり、「混合不⾜」は不安定な状況での特徴と考えられた。2 値化②の場合の結果を図 2 に⽰す。「混合過剰」が現れず、「混合不⾜」は⽚脚⽴位・両脚⽴位共に⾒られ、「右」「左」は固まる傾向にあった。
⽴位の安定度によって連検定の結果が異なる傾向が⾒られた。さらに、従来の評価指標である軌跡
⻑と⽐較したところ、2 値化①のデータに対する連検定の結果との関連がうかがえた。加速度の向き変化のランダム性を調べることにより、動揺制御の質的評価の可能性が⽰唆された。

図 1 加速度の向き変化が同⽅向か逆⽅向かで 2 値化した場合の連検定の結果

図 2 加速度の向き変化が右側か左側かで 2 値化した場合の連検定の結果
⽂献
1) JD Gibbons, S Chakraborti : Nonparametric statistical inference. Chapman and Hall/CRC, 2020

GMANOVA モデルにおける新たな経時変動の推定方法とその解釈 (報告書)中京大学 教養教育研究院 永井 勇

本講演では, n 個の各個体に対して, 全ての個体で測定時点を揃えて p 回測定して得られる

経時測定データの分析を考えた．このようなデータはバランス型経時測定データと呼ばれ, 各個体で測定時点が揃っていないものはアンバランス型経時測定データと呼ばれる. これらのデータの分析の目的は, データの裏に潜む経時変動を上手く捉えることある. 本講演では, バ
ランス型経時測定データの分析について考えた.

バランス型経時測定データの分析の際には, Pothoff and Roy (1964) で提案された次の一

般化多変量分散分析 (Generalized Multivariate Analysis of Variance; GMANOVA) モデルがよく使われるため, 本講演でもこのモデルを用いた;

Y = 1n µ! X ! + AΞX ! + E,

(1)

ここで, 1n は n 次元の全てが 1 からなるベクトル, 0r は r 次元の全てが 0 からなるベクトル,

Y は各行が各個体で測定して得られる経時測定データからなる n × p 行列, A は各個体の特徴を表す測定時点に無関係な k 個の変数からなる rank(A) = k の n × k 行列とし, A! 1n = 0k

(各説明変数で中心化されている) を満たしているとし, X は後述のように各行が測定時点の関数からなる p × q 行列であり, これらの Y , A, X は既知である. また, µ は q 次元未知ベク
トル, Ξ は k × q 未知行列であり, E は E[E] = 0n 0!p , Cov[vec(E)] = Σ ⊗ In の n × p 誤差行列とし, Σ は正則な p × p 未知行列とする. このモデルにおいて, E[Y ] = 1n µ! X ! + AΞX !

の部分が, 経時測定データの分析で目的としている経時変動に対応していることを報告した.
また, 測定時点を t1 < t2 < · · · < tp として, X の i 列目を (t0i , t1i , . . . , tiq−1

) とすることは, 経時変動を測定時点の (q − 1) 次多項式で推定することに対応することを報告した.

このモデル (1) において, 未知の µ と Ξ の推定としてよく使われる推定量は, 次のリスクを最小にする µ と Ξ を求めることで得られることを報告した;
!"
#
"
#! $
R(µ, Ξ|Σ) = tr Y − 1n µ! X ! − AΞX ! Σ−1 Y − 1n µ! X ! − AΞX !
.

(2)

実際に µ̂Σ = arg minµ R(µ, Ξ|Σ) と Ξ̂Σ = arg minΞ R(µ, Ξ|Σ) を求めると, (X ! Σ−1 X)µ̂Σ =X ! Σ−1 Y ! 1n /n, Ξ̂Σ X ! Σ−1 X = (A! A)−1 A! Y Σ−1 X の解となる. これらの方程式を解き µ̂Σや Ξ̂Σ を得るために, モデル (1) において rank(X) = q を仮定することが多い. この仮定は経時変動の推定に用いる関数を制限していると考えられることを報告した.

もし rank(X) = q ならば (X ! Σ−1 X)−1 が存在するので, µ̂Σ = X ! Σ−1 Y ! 1n (X ! Σ−1 X)−1 /n,Ξ̂Σ = (A! A)−1 A! Y Σ−1 X(X ! Σ−1 X)−1 となる. 実際に経時変動を推定する際は Σ が未知のため, その不偏推定量 S = Y ! {In − 1n 1!n /n − A(A! A)−1 A! }Y /(n − k − 1) が代わりに使われる. このとき, E[µ̂S ] = µ, E[Ξ̂S ] = Ξ であることが知られていることを報告した.

一方 rank(X) < q の場合, (X ! Σ−1 X)−1 が存在しないため, これらの推定量は得られない.

そこで本講演では, 永井 (2021, 2022) と同様, q1 + · · · + qr = q となる正の整数 qi を用いて,µ = (µ1 , . . . , µr )! (µi ; qi 次元ベクトル), X = (X1 , . . . , Xr ) (Xi ; rank(Xi ) = qi の p × qi 行

列), Ξ = (Ξ1 , . . . , Ξr ) (Ξi ; k × qi 行列) として, モデル (1) を次のように書き換えた;

Y =

r
!

1n µ!i Xi! +

i=1

r
!
i=1

AΞi Xi! + E.

(3)

ここで Xi の各列が 0p でなければ, r = q (つまり 1 列ずつ分ける) とすると, rank(Xi ) =qi = 1 (i = 1, . . . , r) となる. つまり, X の各列が 0p でない場合も含んでいることを報告した.

このように書き換えたモデルにおいても, リスク (2) と同様のリスクを最小にするようなµ1 , . . . , µr や Ξ1 , . . . , Ξr を求めて並べれば, rank(X) < q でも µ および Ξ の推定量が得られるというのが, 本講演のアイデアであった. 実際, リスク (2) の µ と Ξ と X に分割した形のものを代入すれば, モデル (3) に対応したリスク R! ができ, そのリスクを最小にする µ1 , . . . , µr
および Ξ1 , . . . , Ξr を求めればよいこととなった.

ここでリスク R! を最小にするそれぞれの推定量を µ̂i·Σ , Ξ̂i·Σ (i = 1, . . . , r) として, これらを求めることを考えた. これらを求めて (µ̂1·S , . . . , µ̂r·S ) や (Ξ̂1·S , . . . , Ξ̂r·S ) とすると, µ とΞ の推定量が得られ経時変動の推定ができる. このリスク R! を A! 1n = 0k であることに注意して展開すると, R! =(µ1 , . . . , µr に関連する項) + (Ξ1 , . . . , Ξr に関連する項) と分割できることを報告した. しかし, そのままでは µ̂i·Σ や Ξ̂i·Σ が求まりそうにないため, 永井 (2021)のアイデアを用いた. その結果, それぞれの項を最小にする µ̂1·Σ , . . . , µ̂r·Σ と Ξ̂1·Σ , . . . , Ξ̂r·Σを求めると, 各 ! で以下となることを報告した;



!1
!
Y
n
µ̂!·Σ = (X!! Σ−1 X! )−1 X!! Σ−1 −
Xj µ̂j·Σ  ,
n
j<!


!
Ξ̂!·Σ = (A! A)−1 A! Y −Ξ̂j·Σ Xj!  Σ−1 X! (X!! Σ−1 X! )−1 ,
j<!

ここで,

!

Xj µ̂j·Σ = 0p ,!

Ξ̂j·Σ Xj! = 0k 0!p である. これは, ! = 1 の場合の両方が陽に求まり, それを代入することで ! = 2 のときの推定量が得られることを示していることを報告し
j<1

j<1

た. これらにおいて Σ の部分を S に置き換えることで推定量が構築できる. この推定手法は
永井 (2021) と同様に, 何らかの罰則などを用いていないため最適化のための反復計算が不要である. さらに, rank(X) に関する制約なしで推定が可能となるので, 経時変動の推定に用いる関数への制約がほぼ全てなくなることを表していることを報告した. また, これらの推定量をモデル (3) へ入れることで, 経時変動が推定できる.

これらの推定量の解釈などについては, 当日の講演で報告した.

引用文献:[1] Pothoff, R. F. & Roy, S. N. (1964) A generalized multivariate analysis of variance model
useful especially for growth curve problems. Biometrika, 51, 313–326.

[2] 永井 勇 (2021) 高次元小標本における多変量線形回帰モデルでの推定法, 2021 年度統計関連学会連合大会[3] 永井 勇 (2022) 説明変数がランク落ちしている状況での多変量線形回帰における不偏推
定量, 2022 年度統計関連学会連合大会

統計科学におけるエントロピー概念の誤り訂正
得丸久文 Tokumaru Kumon

1.

意味と恣意的に結合する「記号」から，意味が一般性・普遍性をもつ「概念」への進化
パブロフの条件反射実験が示すように，視覚・聴覚・体感などの感覚から入力される記号は，その直後に経験する餌や毒物の記憶と結びつく．この記憶が言葉の意味であり，記号と意味は個体の経験にもとづいて恣意的(経験的・偶然)に結合する．
文字が生まれ法律が公布されると，法的安定性のために，言葉に一般性が求められるようになり，概念が生まれた．概念は，数学的な群の論理で言葉を例外なくすべての記憶と結合するもの
であり，例外をなくすために「法の番人」として法律家や裁判制度が生まれた．

2.

科学概念の起源論・学際性・有用性を検証する番人の必要性
科学概念の一般性は，起源を検証することで確かめられる．すべての科学概念は，ある特定の科学者が，それまで誰も気づかなかった不可視の現象を発見し，検証し，命名して誕生する．学際的に用いられる科学概念の起源は共通である．科学概念は，概念操作によって，不可視の現象が相互に生みだす複雑な現象を正しく把握することを可能にする．
法概念と同様に，科学概念の起源や有用性を検証する科学概念の番人が必要である．

3.

統計科学のエントロピー概念
熱力学・情報理論を問わず，エントロピー概念は混乱している．その原因は，概念の有用性を問わず，他の分野の概念と矛盾や齟齬があっても気にしない学会の風潮にある．個々の科学者が，概念の発生起源を調べず，間違った言語情報を鵜呑みにするため，雑音や歪が重畳し錯乱し，伝言ゲーム状態がおきている．
エントロピーを概念化するにあたってむずかしいのは，不可視の現象であることに加えて，熱力学でも定義が確立されていない．エントロピー概念とトレードオフの関係にある情報量概念が，誤ってエントロピーと呼ばれている．真逆な現象を同じ名前で呼ぶため混乱に拍車がかかる．

① エントロピーの誕生と概念化（熱力学）エントロピーは，熱力学分野でクラウジウスが発見したが，これまで起源に遡った検証が行われていなかった．八木江里氏が，クラウジウスの方程式を解き直すことで概念化し「構成要素の配置の変化量」と定義した．変化量だから⊿で示され正の値しかとらない．(八木他 2013)② エントロピー概念の学際性（情報理論）
第二次世界大戦中に生命科学や脳科学と接点をもった数学者フォン・ノイマンは，生命や知能を科学する情報理論の構築をめざした．彼は真空管式コンピュータを製造した経験から，情報理1

論は熱力学的エントロピーによって決定づけられることを直感した．ところが彼は情報理論を完成させる前に，創始者の栄誉をクロード・シャノンに譲り，1948 年時点の知見を渡した．このためシャノン理論には，ノイマンが 1950 年以降行った誤り訂正が反映されていない．

③ エントロピー概念の有用性（進化学，電波工学など）
科学概念は，それが群を構成することを確かめる．群は演算に閉じているので，概念操作の結果も群となる．不可視の現象を考える上で，おおまかな全体構造を示す参照モデルが有効である．
この状態で熱力学エントロピーを意識すると，まったく新しい仮説・知見を得ることができる．

4.

情報理論のエントロピー概念における矛盾の解明「通信の数学的理論」が掲げる一般通信モデルは，通信回線上を伝搬する信号が，回線雑音によって劣化する不可視の量子力学的現象を示す参照モデルである．この二元的関係は信号対雑音比(S/N)と呼ばれ，大気を伝搬中の信号と熱雑音という二つの不可視の変数によって信号波形の乱れ(正規値からのズレ量)が生まれる．これがエントロピーであり，それは「雑音によって生まれる信号の歪(=構成要素の配置の変化量)」である．

図 1 一般通信モデルを使った信号(S)・雑音(N)・エントロピー(⊿S)の相関図
現在通用している情報理論のエントロピーと，熱力学的なエントロピーとの矛盾は，なぜ生まれたのか．米国議会図書館に保管されているノイマンの手稿のなかに，1949 年 12 月にイリノイ大学で行った講演記録があり，エントロピーを情報量と結びつけて論じていた．おそらくノイマ
ンは出版にあたって，概念を精査して誤りを自ら訂正した．そしてそれをシャノンに伝えたのだ
が，シャノンは，ノイマンからの誤り訂正要求を無視した．だから 1982 年の IEEE オーラルヒ
ストリーの言葉があると考えられる．この矛盾解明は統計科学にとっても有益である．

図 2/3 情報理論エントロピーの誤り分析

2

Exploration of Galaxy Evolution via Manifold Learning
Tsutomu T. TAKEUCHI1,2 and Suchetha COORAY1,†
1. Division of Particle and Astrophysical Science, Nagoya University, Nagoya 464-8602, Japan
2. The Research Center for Statistical Machine Learning, the Institute of Statistical Mathematics, Tachikawa, Tokyo 190-8562, Japan† JSPS Research Fellow (DC1)

1. Galaxy Formation and Evolution
Matter in the early Universe was almost uniform, and slightly dense regions evolve by gravity, finally intoa galaxy. It was attempted to develop a theory to deal with the star formation and associated history of heavy
element synthesis, under an assumption that a galaxy has formed from a single, huge gas cloud. While theresearch in this direction was once completed in the first half of 1980s, this was not the end of the studies ofgalaxy evolution. Cosmological research that has progressed in parallel has revealed that galaxies merge andgrow. This indicates that the galaxy evolution is a very complicated process that strongly depends on thedensity of the surrounding galaxies and the gas density. In order to formulate the galaxy evolution, it isnecessary to determine such a huge system of equations. Though astrophysicists have constructed thegoverning equations from the physical laws from the first principle before, such a method is not realistic
anymore when the quantity space exceeds 10 dimensions. Galaxy surveys as of the 21st century provides
hundreds (or even thousands) of physical quantities for hundreds of millions of galaxies, typical big data inboth quality and quantity indeed. The feature space of the galaxy to be analyzed exceeds 100 dimensions.
Therefore, the characterization of the galaxy evolution is no longer possible by the traditional method relying
on phy1ical intuition.

2. Galaxy Manifold
2.1 Rise, fall, and revival of the galaxy manifold
From 1970s to the mid-1980s, classical multivariate analysis methods such as the principal component
analysis (PCA) were used to combine physical quantities of galaxies in a high-dimensional space. Various(logarithmic) linear relations, so-called galactic scaling relations, have been discovered. Research to unify thescaling relations and find the fundamental relationships has led to the concept of galaxy manifolds. However,
the galaxy manifold has once been almost forgotten because the classical PCA could treat only linear relations,and it remained a limited concept, though they are still useful for exploring (log)linear relations of galaxies.
Recently, we discovered a galaxy manifold that expresses the basics of galactic evolution by the Fisher EM
algorithm. Because of its strongly nonlinear spatial structure, it could have never been found in previous studies
based on the classical PCA. To understand the manifold, a more sophisticated method beyond a mereclassification is needed. We focused on a method known as the manifold learning, one of the latest methods ofdata science that is completely different from conventional methodologies

2.2 Galaxy manifold constructed by manifold learning
We adopt the algorithm Isomap and UMAP (Uniform Manifold Approximation and Projection). Isomap

defines the neighboring points by using input-space distance and the distant points as a sequence of “shorthops” between neighboring points. Isomap tries to find shortest paths in a graph with edges connecting
neighboring data points. By construction, Isomap preserves the “surface density” of data points in the feature
space. UMAP is based on differential geometry and algebraic topology. The algorithm is founded on threeassumptions: 1) the data are uniformly distributed on a Riemannian manifold, 2) the Riemannian metric is
locally constant (or can be approximated as such), and 3) the manifold is locally connected. From theseassumptions it is possible to model the manifold with a fuzzy topological structure. Sine it defines the manifold
so that the data points distribute as homogeneously as possible, it does not preserve the surface density of data
points. UMAP also preserves some important structural properties, and it is more robust against noise than
Isomap.

Manifold learning algorithm can “unfold” a curved and/or rolled manifold in the feature space, andprovide a local coordinate system on it. The resulting manifolds with local coordinates from Isomap and
UMAP are presented in Fig. 1. From Figure 1, we clearly see that the galaxy manifold is two-dimensional. We
also stress that two different algorithms, Isomap and UMAP yield similar two-dimensional manifolds. The
difference of the two estimated manifolds is clearly seen in Fig. 1. Since Isomap preserves the density of data
point cloud, we observe that the manifold has a density structure, i.e., dense and sparse regions on the manifold.
Figure 1: The “unfolded” galaxy manifold bya manifold learning algorithm Isomap and
UMAP. Left and right panels show the
manifolds

from

Isomapand

UMAP,

respectively. Though the global shape isslightly different from each other, they sharecommon features on the manifold.

2.3 Result
The galaxy manifold obtained with Isomap preserve thisinformation and reveal the speed of galaxy evolution at various
stages along the manifold. e.g., galaxies passes the green valley
very fast. In contrast, the galaxy manifold obtained with UMAP
is imposed uniformity on the galaxy data, leading to a morerobust and representative description of the observed galaxy
properties e.g., galaxies evolve continuously in the feature space,
without a discontinuity or “jump” on their evolutionary tracks.
Thus, the galaxy manifold provides a clue to the evolutionary
path of galaxies on the manifold. The SFR and stellar mass fields

Figure 2: The vector field of star formation
rate on the galaxy manifold.

do not show the same evolutionary path. This supports that the galaxy merger without star formation plays asignificant role in the growth of stellar mass. Next step is to fully parametrize the evolution equation of galaxies.
For example, we can construct a vector field of the star formation rate on the galaxy manifold (Fig. 2).

Fréchet ‫཭ڑ‬Λ༻͍ͨ෼ྨ໰୊͍ͭͯҏ৿ ৾ฏ 1 , ए໦ ޺จ 1

1: ޿ౡେֶ େֶӃઌਐཧ޻‫ܥ‬Պֶ‫ڀݚ‬Պ

1

͸͡Ίʹ

ຊൃදͰ͸ɼଟ‫܈‬ͷࠞ߹σʔλʹ͓͚Δਪଌ໰୊Λѻ͏ɽಛʹ֤‫܈‬ͷ฼ूஂύϥϝʔλͷਪఆʹ‫ڵ‬ຯ͕͋Δ
৔߹ɼਖ਼֬ʹ‫܈‬Λ෼ྨ͢Δ͜ͱͰਪఆ݁Ռͷ޲্͕‫ظ‬଴Ͱ͖ΔɽMercatanti et al. (2015) Ͱ͸ɼ2 ‫܈‬ͷࠞ߹ਖ਼‫ن‬෼෍Ϟσϧʹ͓͍ͯɼ‫ڵ‬ຯͷର৅Ͱ͋Δओཁม਺ʹิॿม਺ΛՃ͑ɼม਺ͷ࣍‫ݩ‬Λ֦ு͢Δ͜ͱʹΑΔɼਪଌਫ਼౓ͷ޲্ʹ͍ͭͯٞ࿦͍ͯ͠Δɽ·ͨɼ͍͔ͭ͘ͷԾఆͷԼͰɼͲͷΑ͏ͳิॿม਺ΛՃ͑Ε͹֤‫ݸ‬ମΛਖ਼͘͠෼ྨ͢Δ֬཰্͕ঢ͢Δ͔ʹ͍ͭͯɼཧ࿦తͳߟ࡯Λ༩͍͑ͯΔɽ۩ମతʹ͸ɼิॿม਺
ؒͷฏ‫ॿิͱࠩۉ‬ม਺ͱओཁม਺ͷ૬ؔͷେ͖͕͞ॏཁͰ͋Δ͜ͱΛ͍ࣔͯ͠Δɽ
ҰํͰɼ͜ͷิॿม਺ͷ༗༻ੑʹؔ͢Δ݁Ռ͸ղੳͷ໨తͱख๏ɼσʔλͳͲͷઃఆʹґଘ͓ͯ͠ΓɼҟͳΔঢ়‫͍͓ͯʹگ‬΋ɼಉ͡ಛ௃Λ࣋ͭิॿม਺͕໾ཱ͔ͭ͸Θ͔Βͳ͍ɽImori and Shimodaira (2019)Ͱ͸ɼෆ‫׬‬શσʔλղੳʹ͓͚Δ༗༻ͳิॿม਺ͷબ୒໰୊ͱͯ͠ɼ‫׬‬શσʔλʹ͓͚Δଛࣦؔ਺ΛఆΊɼͦΕʹ‫ظͮ͘ج‬଴ଛࣦؔ਺ͷਪఆྔͱͯ͠ɼAIC (Akaike, 1974) λΠϓͳͲͷม਺બ୒‫ن‬४Λಋग़͓ͯ͠
ΓɼͦΕʹΑͬͯ༗༻ͳิॿม਺ͷબ୒ΛࢼΈ͍ͯΔɽ
ຊൃදͰ͸ɼ֤‫܈‬ͷ฼ूஂ෼෍ͷ෼ྨ͓ΑͼͦͷύϥϝʔλͷਪఆΛ໨తͱ͠ɼਖ਼‫ن‬෼෍ؒͷ‫͠ͱ཭ڑ‬ͯɼFréchet ‫( ཭ڑ‬Dowson and Landau, 1982) Λ༻͍ͨ৔߹ͷ෼ྨ໰୊ʹର͠ɼ༗༻ͳิॿม਺ͷ৚݅౳ʹ͍ͭͯߟ࡯ͨ͠ɽ

2Fréchet ‫཭ڑ‬
(i)

(i)

(i)

2 ͭͷ p ࣍‫ݩ‬ଟมྔਖ਼‫ن‬෼෍ Πx : Np (µx , Σx ) (i = 1, 2) ͷ෼ྨΛߟ͑ΔɽDowson and Landau (1982)
(1)

(2)

ʹΑΓɼ2 ͭͷਖ਼‫ن‬෼෍ؒͷ Fréchet ‫ ཭ڑ‬F (Πx , Πx ) ͕࣍ͷΑ͏ʹಋग़͞Ε͍ͯΔ:
(2)
(1)
(2) 2
(1)
(2)
(1) (2) 1/2
F 2 (Π(1)).
x , Πx ) = !µx − µx ! + tr(Σx + Σx − 2(Σx Σx )
(1)

(2)

(1)

(2)

ͨͩ͠ɼ(Σx Σx )1/2 ͸ඇෛͷ‫ݻ‬༗஋Λ࣋ͭ Σx Σx

(1)

(2)

ͷฏํࠜߦྻͰ͋Γɼtr((Σx Σx )1/2 ) ͸Ұҙʹ

ఆ·Δ͜ͱʹ஫ҙ͢Δɽ
(i)͍·ɼ෼෍ Πx Λ฼ूஂ෼෍ͱ͢Δ֬཰ม਺ xi (i = 1, 2) Λ࣍ͷΑ͏ʹ෼ׂ͢Δ: 
yixi =   .
aiͨͩ͠ɼyi ∈ Rq (q < p) Λओཁม਺ɼม਺ ai ∈ Rr (r = p − q) Λิॿม਺Ͱ͋Δͱ͢Δɽ͕ͨͬͯ͠ɼyiͷपล෼෍ͷਪଌʹ‫ڵ‬ຯ͕͋Δ৔߹ʹɼai ΛՃ͑ͨ xi Λ༻͍ͯ෼ྨΛߦ͏ɽ

1ͯ͞ɼ‫ؒ܈‬ͷ‫͕޿͕཭ڑ‬Δͱɼ෼ྨਫ਼౓͕޲্͠ɼͦΕʹ൐͍ύϥϝʔλͷਪఆਫ਼౓΋޲্͢Δͱ‫ظ‬଴Ͱ͖Δɽ൓ରʹɼ͋Δม਺ΛՃ͑ͯ΋‫ؒ܈‬ͷ‫͕޿͕཭ڑ‬Βͳ͍৔߹ɼͦͷม਺͸ϊΠζʹ͔͠ͳΒͳ͍ͱߟ͑
ΒΕΔɽͦ͜ͰຊൃදͰ͸ɼFréchet ‫ʹ཭ڑ‬ର͢Δ৑௕ੑΛ΋ͱʹ༗༻Ͱͳ͍ิॿม਺ΛఆΊΔɽ͢ͳΘͪɼ
(i)

(i)

(i)

(i)

(i)

i = 1, 2 ʹର͠ɼµy = E[yi ], Σy = Cov[yi ] ͱఆΊͨͱ͖ɼq ࣍‫ݩ‬ଟมྔਖ਼‫ن‬෼෍ Πy : Nq (µy , Σy )
(i)

(i)

ͷ Fréchet ‫཭ڑ‬Λ F (Πy , Πy ) ͱ͢Δͱɼ
(2)
(1)
(2)
F (Π(1)y , Πy ) = F (Πx , Πx )Ͱ͋Δͱ͖ɼai Λ৑௕ͳʢ༗༻Ͱͳ͍ʣิॿม਺Ͱ͋ΔͱΈͳ͢ɽ࣮ࡍͷ෼ྨͰ͸ Fréchet ‫཭ڑ‬Λߏ੒͢Δύϥϝʔλ͸ະ஌Ͱ͋Δ͔Βɼ֤ύϥϝʔλͷਪఆྔΛ୅ೖ͢Δ
͜ͱͰ Fréchet ‫཭ڑ‬Λਪఆ͢Δඞཁ͕͋ΔɽຊൃදͰ͸ɼ‫؍‬ଌσʔλͱͯ͠ɼXi = (x1i , . . . , xni )! ∈ Rn×p ,
(i)

(i)

ͨͩ͠ xti ∼ Np (µx , Σx ) Λಘ͓ͯΓɼxti (i = 1, 2; t = 1, . . . , n) ͕ಠཱͰ͋Δঢ়‫گ‬Λ૝ఆ͠ɼX1 ͱ X2

ʹ‫ ͖ͮج‬Fréchet ‫཭ڑ‬Λਪఆͨ͠ࡍͷऩଋϨʔτʹ͍ͭͯߟ࡯ͨ͠ɽ۩ମతͳ݁Ռ͓Αͼ਺஋γϛϡϨʔγϣϯʹ͍ͭͯ΋౰೔ใࠂͨ͠ɽ

ࢀߟจ‫ݙ‬
1. Akaike, H. (1974). A new look at the statistical model identification. IEEE transactions on automatic control, 19, 716–723.
2. Dowson, D. C., & Landau, B. (1982). The Fréchet distance between multivariate normal distributions. Journal of multivariate analysis, 12, 450–455.
3. Imori, S., & Shimodaira, H. (2019). An information criterion for auxiliary variable selection inincomplete data analysis. Entropy, 21, 281.
4. Mercatanti, A., Li, F., & Mealli, F. (2015). Improving inference of Gaussian mixtures using auxiliary
variables. Statistical analysis and data mining, 8, 34–48. ɹ2

正方分割表における松下距離に基づく周辺同等性からの隔たり尺度
篠田
吉本
田畑

1

はじめに

2

提案した尺度覚
拓矢
耕治(横浜市立大学)(中外製薬株式会社)(東京理科大学)

行と列が同じ分類からなる正方分割表は, 医学, 教育学, 社会科学等の分野で活用されている. このような正方分割表の解析では, 多くの観測度数が主対角セルまたはその近傍に集中するため, 行と列の分類間の独立性は成り立たない. したがって, 我々は独立性ではなく対称性または周辺同等性に関心がある. たとえば, 周辺同等 (MH)モデルは行変数と列変数の各周辺分布が等しいかどうかを表す (Stuart, 1955). MH
モデルの当てはまりが悪いとき, どの程度 MH モデルからの隔たりがあるかに関心があり, これまでに複数の隔たりを測る尺度が提案された (Tomizawa et al., 2003;
Yamamoto et al., 2011; Ando et al., 2021).
近年, 質的データの視覚化に関心が高まっている. しかしながら, 量的データの可視化とは異なり, 質的データの可視化の歴史は浅く, データ解析の際に活用される機会も少ない (Friendly and Meyer, 2015).
本講演では, 複数の隔たりを測る尺度に関する先行研究を整理し, 視覚的な解釈を容易にする MH モデルからの隔たりを測る尺度を提案した.

行と列が同じ分類からなる r×r 正方分割表において, 行変数を X, 列変数を Y とし,(i, j) セル確率を pij とする (i = 1, . . . , r; j = 1, . . . , r). また, 任意の i = 1, . . . , r − 1
に対して,
i
r
i
r
!
!
!
!
G1(i) =
pst , G2(i) =pst ,s=1 t=i+1

Gc1(i) =
∆=
とする.

r−1
!
"
i=1

t=1 s=i+1

G1(i),G1(i) + G2(i)#G1(i) + G2(i) ,Gc2(i) =G∗1(i) =G2(i),G1(i) + G2(i)

G1(i)
,
∆

G∗2(i) =

G2(i)
,
∆

上記を用いて, MH モデルからの隔たりを測る尺度を次のように提案した : {G1(i) +G2(i) > 0} を仮定し,
r−1
!
" ∗
#
∗
Φ =
G1(i) + G∗2(i) Mi ,
ただし,

i=1


 12
+ ,2 )+ ,2 
√ )*

*
2+ 2
1
1
 .
Mi = 
Gc1(i) −+Gc2(i) −

2
2
2 


この提案尺度 Φ∗ は {Gc1(i) , Gc2(i) } と {1/2, 1/2} の間の松下距離 Mi の重み付き和であり, Mi は距離の定義を満たす (Matusita, 1954, 1955; Read and Cressie, 1988, p.112).
したがって, 視覚化において自然な解釈が可能となる尺度であると考えた.
さらに, 提案尺度に対して質的データの可視化方法を参考に視覚化を行なった. 提案尺度および視覚化の有用性については, 医学分野の実データ解析を用いて示した.
また, 尺度の近似信頼区間を導出し, シミュレーションによって近似信頼区間の被覆確率を評価した.

参考文献[1] Ando, S., Noguchi, T., Ishii, A. and Tomizawa, S. (2021). A two-dimensional index for
marginal homogeneity in ordinal square contingency tables. SUT Journal of Mathematics,
57, 211-224.
[2] Friendly, M. and Meyer, D. (2015). Discrete Data Analysis with R: Visualization and ModelingTechniques for Categorical and Count Data. Chapman and Hall.
[3] Matusita, K. (1954). On the estimation by the minimum distance method. Annals of theInstitute of Statistical Mathematics, 5, 59-65.
[4] Matusita, K. (1955). Decision rules based on the distance, for problems of fit, two samples,and estimation. Annals of the Institute of Statistical Mathematics, 26, 631-640.
[5] Read, T.R.C. and Cressie, N. (1988). Goodness-of-Fit Statistics for Discrete Multivariate
Data. Springer, New York.
[6] Stuart, A. (1955). A test for homogeneity of the marginal distributions in a two-way classification. Biometrika, 42, 412-416.
[7] Tomizawa, S., Miyamoto, N. and Ashihara, N. (2003). Measure of departure from marginal
homogeneity for square contingency tables having ordered categories. Behaviormetrika, 30,
173-193.
[8] Yamamoto, K., Ando, S. and Tomizawa, S. (2011). A measure of departure from averagemarginal homogeneity for square contingency tables with ordered categories. Revstat, 9, 115126.

方向データのためのカーネル密度推定量のバイアス修正
長野県立大学グローバルマネジメント学部助教

鶴田靖人

1 はじめに

方向データは，d 次元球面 Sd := {x ∈ Rd+1 : ||x|| = 1} 上に観測点が位置するデータである．方向データの例

として，次元 d = 2 のときは過去の地磁気のデータなどがあり，d ≥ 3 のときはテキストデータ（の各文書が含む単語総数を 1 に基準化したもの）がある．また，d = 1 のときは円周上のデータと呼ばれており，例えば風向などの角度データがある．カーネル密度推定法などのノンパラメトリック密度推定法の利点は柔軟な推定が可能である．方向データのためのカーネル密度推定量の初期の研究は，Beran (1979), Hall et al. (1987) や Bai et al.

(1989) が挙げられる．方向データ X1 , . . . , Xn は独立同一分布 f (x) (x ∈ Sd ) に従うとする．f についてのカー
ネル密度推定量は，n1!
fˆκ (x) :=Kκ (xT Xi ),n i=1

ただし，Kκ (xT Xi ) : Sd → R は球面上のカーネル関数，κ := κ(n) は limn→∞ κ = ∞ かつ limn→∞ n−1 κd/2 = 0
を満たす平滑化パラメータとする．ここで，η ∈ Sd とおく． カーネル関数を Kκ (xT η) := Cκ (L)−1 L(κ(1−xT η))と定義する．ただし，L(r) : [0, ∞) → R は微分可能とし，Cκ (L) :=L(κ(1 − xT η))wd (x) は基準化定数を表"す．よく用いられる誤差基準は，平均積分二乗誤差 (MISE) MISE[fˆκ (x)] := Sd E[(fˆκ (x) − f (x))2 ]wd (dx) であ"Sd

る（ただし，wd は d 次元球面上のルベーグ測度を表す）．f が十分に滑らかであるという仮定の下で，MISE の収束レートは O(n−4/(4+d) ) であることが知られている (Hall et al. 1987)．つまり，次元 d が増加することで MISEの収束レートは悪化する（これを次元の呪いという）．

本研究の目的は，次元の呪いの影響を減らすためにカーネル密度推定量 fˆ のバイアスを修正し，MISE の収束レートを改良することである．そのために，実数空間上のバイアス修正法である Jones and Foster 型修正法

(Jones and Foster 1993) と Terrell and Scott 型修正法 (Terrell and Scott 1980) を d ≥ 2 次元の球面上のカーネ
ル密度推定量に適用し，バイアスを修正可能なことを示す．本研究の提案手法は，Tsuruta and Sagae (2017) が円周上 (d = 1) のカーネル密度推定量のために提案したバイアス修正法を d ≥ 2 次元に拡張したものである．

2 提案手法

期待値 E[[fˆκ (x)] は，f をテイラー展開すると κ のべき乗に L のモーメント µl (L) :="∞0

L(r)r(l+d−2)/2 dr をかけた項の線形和になる．したがって，次式のように低次のモーメントが 0 になるようなカーネル関数はバイアスを修正する．

µ0 (L) &= 0,µl (L) = 0,l = 2, 4, 6, . . . , p − 2,

µp (L) &= 0,(1)

ただし，p ≥ 2 は偶数である．(1) を満たすカーネル関数を p 次オーダーカーネルと呼ぶ (Tsuruta and Sagae

2017)．p 次オーダーカーネルは以下のような性質を持つ．

定理 1. 緩やかな仮定の下で，p 次オーダーカーネルを採用したカーネル密度推定量のバイアスと分散はそれぞれ，

Bias[fˆκ (x)] = O(κ−p/2 ) と Var[fˆκ (x)] = O(n−1 κd/2 ) となる．また，平滑化パラメータとして κ∗ = O(n2/(2p+d) )を選択すると MISE[fˆκ (x)] の収束レートは O(n−2p/(2p+d) ) となる．

p 次オーダーカーネルの構成法として，以下のような Jones and Foster 型修正法を提案する．

p 次オーダーカーネルの関数 L を L[p] と表す．このとき，p + 2 次オーダーカーネルの関数 L を次式のように定義する．

L[p+2] (r) :=p+d2
L[p] (r) + rL![p] (r).
ppJones and Foster 型修正法は低次のモーメントが 0 となるように２つの関数 L を足し合わせたものである．低次のモーメントが 0 となるように２つの異なるカーネル密度推定量をかけ合わせることでもバイアスを修正することができる．そのようなバイアス修正法として，Terell and Scot 型のバイアス修正法を提案する．ここで，fˆκ,[p] (x) は p 次オーダーカーネルを採用したカーネル密度推定量を表す．
Terell and Scot 型のカーネル密度推定量を次式のように定義する．

ただし， a ∈ (0, 1) とする.

!
"1/(1−a) !
"−a/(1−a)f˜κ (x) = fˆκ,[2] (x)fˆaκ,[2] (x),Terell and Scot 型のカーネル密度推定量は以下の性質を持つ．

定理 2. 緩やかな仮定の下で，Terell and Scot 型カーネル密度推定量のバイアスと分散はそれぞれ，Bias[f˜κ (x)] =O(κ−2 ) と Var[f˜κ (x)] = O(n−1 κd/2 ) となる．また，平滑化パラメータとして κ∗ = O(n2/(8+d) ) を選択するとMISE[f˜κ (x)] の収束レートは O(n−8/(8+d) ) となる．
したがって，Terell and Scot 型のカーネル密度推定量は MISE の収束レートが 4 次オーダーカーネルに対応するカーネル密度推定量であると言える．Terell and Scot 型のカーネル密度推定量はパラメータ a の最適値を求める必要がある．そこで a に依存しないバイアス修正法として，指数型のバイアス修正法を提案する．
指数型のカーネル密度推定量を次式で定義する．

fˇκ (x) = fˆκ,[2] (x) exp#$fˆκ,[4] (x)−1fˆκ,[2] (x)なお，指数型のカーネル密度推定量は，緩やかな仮定の下でバイアス，分散と最適な MISE のオーダーはそれぞれ定理 2 と同じオーダーとなる．

参考文献
ZD Bai, C Radhakrishna Rao, and LC Zhao. Kernel estimators of density function of directional data.
Multivariate Statistics and Probability, 27:24–39, 1989.
Rudolf Beran. Exponential models for directional data. The Annals of Statistics, 7(6):1162–1178, 1979.
Peter Hall, GS Watson, and Javier Cabrera. Kernel density estimation with spherical data. Biometrika, 74(4):751–762, 1987.
MC Jones and PJ Foster. Generalized jackknifing and higher order kernels. Journal of Nonparametric
Statistics, 3(1):81–94, 1993.
George R Terrell and David W Scott. On improving convergence rates for nonnegative kernel density estimators. The Annals of Statistics, 8(5):1160–1163, 1980.
Yasuhito Tsuruta and Masahiko Sagae. Higher order kernel density estimation on the circle. Statistics &Probability Letters, 131:46–50, 2017.

Outcome-adaptive lasso と adaptive lasso を用いた
AIPW 推定量の性能評価
本江渡 1,2 ，安藤宗司 3 ，土田潤 4 ，寒水孝司 51

東京理科大学大学院工学研究科情報工学専攻，2 ノバルティスファーマ株式会社，3

東京理科大学理工学部情報科学科，4 同志社大学文化情報学部文化情報学科，
5

東京理科大学工学部情報工学科

近年，医療情報データベースを用いて因果効果を推定する研究が実施できるようになった．データには多くの変数が含まれているため，交絡変数の候補となる変数の数が多くなる場合がある．交絡による因果効果の推定量のバイアスを避けるためには，興味のあるアウトカム・暴露・患者背景等の変数間の因果関係に基づいて交絡の調整に用いる変数を研究計画時に特定することが望ましいが，変数の数が多いとすべての交絡変数を研究計画時に特定することは難しい．このような場合，データ駆動型の変数選択法が一つの解決策となりうる [1, 2]．
交絡変数の調整に加えてアウトカムのみに関連する変数を調整すると因果効果の推定量の分散を小さくすることができる．アウトカムに関連する変数（交絡変数とアウトカムのみに関連する変数）を選択する傾向があるデー
タ駆動型の変数選択法として，outcome-adaptive lasso がある [1]．Outcome-adaptive lasso は傾向スコアを推定するためのロジスティック回帰モデルのパラメータ推定のために提案されている．Outcome-adaptive lasso により算出される傾向スコアに基づく inverse probability weighted（以下，IPW）推定量は，アウトカムに関連する
変数を選択する傾向がない既存のデータ駆動型の変数選択法により算出される傾向スコアに基づく IPW 推定量
よりも，分散が小さいことが数値実験によって示されている．しかしながら，IPW 推定量は，傾向スコアモデルを誤特定した場合にバイアスが生じ，傾向スコアの推定値が極端な値になると推定が不安定になる．このような問題に対処するために，IPW 法にアウトカム回帰を組み合わせた二重頑健推定量を用いることができる [3].
二重頑健推定量として，傾向スコアモデルとアウトカム回帰モデルを用いる augmented inverse propensity

weighted（以下，AIPW）推定量 [5, 6] がよく使用される．通常，二重頑健推定量に用いる傾向スコアモデルとアウトカム回帰モデルには同じ変数が使用され [3]，異なる変数を使用するとバイアスが生じる可能性がある [4]．

outcome-adaptive lasso と adaptive lasso は，アウトカムに関連する変数に対するパラメータの推定量がオラク
ル性を有する．そのため，それぞれ傾向スコアモデルとアウトカム回帰モデルの構築のために用いることにより，両モデルには同じ変数が選択される傾向があると考えられる．
本研究の目的は，傾向スコアモデルとアウトカム回帰モデルのパラメータをそれぞれ outcome-adaptive lassoと adaptive lasso で推定する AIPW 推定量が傾向スコアモデルのパラメータを outcome-adaptive lasso で推定す
る IPW 推定量よりも性能が向上するかを評価することである．AIPW 推定量を構築する際に outcome-adaptive

lasso と adaptive lasso を組み合わせることによる良さを評価するため，比較対象として通常の lasso を用いた単純な AIPW 推定量と Farrell の推定量 [7] を含めた．Farrell の推定量は，傾向スコアモデルとアウトカム回帰モ
デルの変数選択のためだけに lasso を適用し，lasso により選択された変数のパラメータを最尤法により推定した
傾向スコアモデルとアウトカム回帰モデルに基づく AIPW 推定量である．比較対象には，adaptive lasso を用いたアウトカム回帰のみに基づく推定量も含めた．実際の臨床研究データの解析と数値実験を通して，これらの推定量の性能を比較した．
実際の臨床研究データの解析として，Vanderbilt 大学が公開しているデータ（https://hbiostat.org/data/）[8]から，右心カテーテル法の集中治療室滞在日数に対する平均因果効果を推定した．これより，outcome-adaptive

lasso を用いた IPW 推定値は他の推定値より高かった．
数値実験では，処置とアウトカムの生成に用いる変数をそのまま傾向スコアモデルとアウトカム回帰モデルの

説明変数に含めるシナリオ（傾向スコアモデルとアウトカム回帰モデルの両方を正しく特定可能なシナリオ）と処置とアウトカムの生成に用いる変数を誤変換して傾向スコアモデルまたはアウトカム回帰モデルの説明変数に含めるシナリオ（傾向スコアモデルまたはアウトカム回帰モデル，あるいはその両方を常に誤特定するシナリオ）を設定した．outcome-adaptive lasso を用いた IPW 推定量は，傾向スコアモデルを正しく特定可能な場合でも，共変量間の相関係数が大きくなるにつれて因果効果の推定量のバイアスと分散が大きくなった．この原因として，
傾向スコアの逆数による重みが極端に大きくなることが挙げられる．傾向スコアまたはアウトカム回帰を単独で用いる平均因果効果の推定量は，そのモデルを誤特定する場合，偏りと分散が大きくなった．単純に lasso を用い
た AIPW 推定量は全てのシナリオにおいてバイアスが大きくなったにも関わらず，Farrell の推定量 [7] は傾向ス
コアモデルまたはアウトカム回帰モデルを正しく特定可能な場合，単純に lasso を用いた AIPW 推定量よりもバ
イアスが小さくなる傾向があった．このことから，単純に lasso を用いると縮小推定のために十分に交絡が調整できないためバイアスが生じたと考えた．outcome-adaptive lasso と adaptive lasso を用いた AIPW 推定量は，傾向スコアモデルまたはアウトカム回帰モデルを正しく特定可能な場合，他の推定量と比べてバイアスと分散が小さくなった．これは，outcome-adaptive lasso と adaptive lasso がオラクル性を有し，傾向スコアモデルとアウト
カム回帰モデルの両方でアウトカムに関連する変数を選択でき，通常の lasso よりも縮小推定が起きにくいこと
で十分に交絡が調整されたためであると考えた．傾向スコアモデルとアウトカム回帰モデルの両方を誤特定するシナリオでは，すべての推定量のバイアスと分散が同程度に大きくなった．

outcome-adaptive lasso を用いた IPW 推定量は，実データ解析では他の推定量による推定値よりも高い値をとり，数値実験ではバイアスと分散が大きくなる傾向があった．このことから，outcome-adaptive lasso を用いた

IPW 推定量は他の推定量よりもバイアスと分散が大きくなりやすいことが示唆された．この理由として傾向スコ
アの逆数による重みが極端に大きくなることが考えられるため，極端な重みが生じにくい stabilized weight を用いることで今回使用した IPW 推定量が改善する可能性がある．傾向スコアモデルまたはアウトカム回帰モデル
を正しく特定可能な場合は，outcome-adaptive lasso と adaptive lasso を用いた AIPW 推定量のバイアスと分散
は，outcome-adaptive lasso を用いた IPW 推定量よりも小さくなった．一方で，両モデルを誤特定する場合は今回検討したすべての AIPW 推定量のバイアスや分散が IPW 推定量やアウトカム回帰による推定量と同程度に大きくなった．

主要参考文献[1] Shortreed S and Ertefaie A. Outcome-Adaptive Lasso: Variable selection for causal inference. Biometrics
2017; 73: 1111–1122.
[2] Antonelli J, Cefalu M, Palmer N. Doubly robust matching estimators for high dimensional confounding
adjustment. Biometrics 2018; 72: 1–9.
[3] Lunceford K and Davidian M. Stratification and weighting via the propensity score in estimation of causal
treatment effects: a comparative study. Statistics in Medicine 2004; 23: 2937–2960.
[4] Shinozaki T and Nojima M. Misuse of regression adjustment for additional confounders following insufficient propensity score balancing. Epidemiology 2019; 30: 541–548.
[5] Robins JM, Rotnitzky A, and Zhao LP. Estimation of regression coefficients when some regressors are notalways observed. Journal of the American Statistical Association 1994; 89: 846–866.
[6] Glynn AN and Quinn KM. An introduction to the augmented inverse propensity weighted estimator.
Political Analysis 2010; 18: 36–56.
[7] Farrell M. Robust inference on average treatment effects with possibly more covariates than observations.
Journal of Econometrics 2015; 189: 1–23.
[8] Connors AF Jr, Speroff T, Dawson NV, Thomas C, Harrell FE Jr, Wagner D, Desbiens N, Goldman L, Wu
AW, Califf RM, Fulkerson WJ Jr, Vidaillet H, Broste S, Bellamy P, Lynn J, Knaus WA. The effectivenessof right heart catheterization in the initial care of critically ill patients. SUPPORT Investigators. Journalof American Medical Association 1996; 276: 889–97.

特異ウィシャート行列の固有値の正確分布論
清水 康希 1,2 橋口 博樹 1

1 東京理科大学

2 日本学術振興会

特別研究員 (PD)

本報告では，変量の次元がサンプルサイズよりも大き場合の標本分散共分散行列，その定数倍である特異ウィシャート行列の正確な分布論について議論した．特異ウィシャート行列の密度関数は Uhlig [7] によって与えられた．Srivastava [4] は特異ウィシャート行列の固有値の同時分布を
Stiefel 多様体上の積分で表現した．Shimizu and Hashiguchi [1] は，その同時分布を関数形で表示するために非斉次超幾何関数を導入し，Sugiyama [5] の積分公式を利用して最大固有値の正確な分布を求めた．最大固有値の分布関数は，特異と非特異を区別することなく統一的に表現できる．
また，Shimizu and Hashiguchi [2] は，特異ウィシャート行列とそれとは独立の非特異ウィシャー
ト行列の比の密度関数，固有値の同時密度関数，最大固有値分布などを求めた．この最大固有値
分布も特異と非特異の区別なく統一的に表現できる．さらに，Shimizu and Hashiguchi [3] では，
球形検定に必要な最大固有値と最小固有値の比の正確な分布を導出した．また，実際に分布の数値計算をするためにゾーナル多項式の積をゾーナル多項式の線形結合で表現するアルゴリズムを提案した．なお，Shimizu and Hashiguchi [1, 2, 3] は，実数を β = 1 として含む β-特異ウィシャー
ト行列に関するものであり，β = 2 の複素数，β = 4 の四元数の場合も含む統一的な議論である．
本報告では，Shimizu and Hashiguchi [1, 3] の成果を実数の場合（β = 1）で紹介した．
多変量正規分布 Nm (0, Σ) に独立に従う m 次元ベクトルを x1 , x2 , . . . , xn とし X = (x1 , x2 , . . . , xn )とする．ただし Σ は m 次の正定値行列とする．このとき，ウィシャート行列を W = XX ! と定
義し，ウィシャート分布を Wm (n, Σ) と表す. n ≥ m である場合の W を非特異ウィシャート行
列，m > n では特異ウィシャート行列と呼ぶことにする．特異の場合の W は H1! H1 = In とな
る m × n 行列 H1 を用いることで，W = H1 L1 H1! とスペクトル分解することができる．ただし，
L1 = diag("1 , . . . , "n ) である．Uhlig [7] は，特異ウィシャート行列 W の密度関数を2

f (W ) =π (−mn+n )/2exp(−trΣ−1 W/2)(detL1 )(n−m−1)/2
2mn/2 Γn (n/2)(detΣ)n/2と与えた．ただし，Γn (·) は多変量ガンマ関数である．m 次対称行列 A に対して，行列変数の超幾
何関数を
p Fq (α; β; X) =∞ !
!
(α1 )κ · · · (αp )κ Cκ (X)k
k=0 κ∈Pm(β1 )κ · · · (βq )κ

k!

(1)

"
と定義する．ただし，α = (α1 , . . . , αp )，β = (β1 , . . . , βq )，(a)κ = m
i=1 (a − (i − 1)/2)κi は一般" κi化ポッホハマー記号であり，(a)κi = j=1 (a + j − 1) である．
Shimizu and Hashiguchi [1, 2] は，特異ウィシャート行列に関する固有値分布の導出で重要となる非斉次超幾何関数を次のように定義した．
定義. A を m 次の対称行列，B を n 次の対称行列とし，m ≥ n とする．このとき，非斉次超幾何関数を以下で定義する．
(m,n)(α; β; A, B) =
p Fq

∞ !
!
(α1 )κ · · · (αp )κ Cκ (A)Cκ (B)
k=0 κ∈Pnk(β1 )κ · · · (βq )κ

k!Cκ (Im )(2)

定理 1. W ∼ Wm (n, Σ)，m > n とする．このとき，W の最大固有値 !1 の分布関数は以下で表される．
!
"
Γn ((n + 1)/2)( x2 )nm/2
1 −1(m,n) m n + m + 1
Pr(!1 < x) =;; − xΣ , In
(3)
1 F1
2
2
2
Γn ((n + m + 1)/2)|Σ|n/2
特異と非特異の場合を区別することなく，ウィシャート行列の最大固有値の正確分布は以下の定理 2 のように統一的に表現できる．

定理 2. 特異・非特異ウィシャート行列 W において，その最大固有値 !1 の分布関数は以下で表される．
!
"
x nm/2Γt ( t+1
n n + m + 1 1 −1
2 )( 2 )Pr(!1 < x) =;; − xΣ
1 F1
2
2
2
Γt ( m+n+1)|Σ|n/2
2!
"!
"
t+1 x nm/2
Γt ( 2 )( 2 )x
m + 1 n + m + 1 1 −1
=
exp
−
trΣ
F
;
;
xΣ
1 1
2
2
2
2
Γt ( m+n+1)|Σ|n/22
ただし，t = min(n, m) である．

Sugiyama [6] は非特異の場合に，ウィシャート行列の最大最小固有値の比 1 − !m /!1 の正確な分布を求めた．m > n の場合は，Shimizu and Hashiguchi [3] で次のように与えられた．
定理 3 (Shimizu and Hashiguchi [3]). W ∼ Wm (n, σ 2 Im )，m > n とする．このとき，x = 1−!n /!1の密度関数は以下のように表される．

f (x) = C∞∞## Γ(mn/2 + k) #(n − 1)(n + 2)/2 + k + tx(n−1)(n+2)/2+k+t−1
kt!
(n k!)kk=0 κ∈Pn−1

×

#

t=0

# {(n − m + 1)/2}τ (n/2 + 1)δ
δgκ,τ
Cδ (In−1 )
(n
+
1)
δ
k+t

t
τ ∈Pn−1δ∈Pn−1

δ はーナル多項式の積に現れる係数である．
ただし，gκ,τ

参考文献[1] Shimizu, K. and Hashiguchi, H. (2021). Heterogeneous hypergeometric functions with two
matrix arguments and the exact distribution of the largest eigenvalue of a singular betaWishart matrix. Journal of Multivariate Analysis 183, 104714.
[2] Shimizu, K. and Hashiguchi, H. (2022). Expressing the largest eigenvalue of a singularbeta F -matrix with heterogeneous hypergeometric functions. Random Matrices: Theory
and Applications.
[3] Shimizu, K. and Hashiguchi, H. (2022). Algorithm for the product of Jack polynomials andits application to the sphericity test, Statistics & Probability Letters. 187, 109505.
[4] Srivastava, M. S. (2003). Singular Wishart and multivariate beta distributions. Annals of
Statistics. 31, 1537–1560.
[5] Sugiyama, T. (1967). On the distribution of the largest latent root of the covariance matrix.
The Annals of Mathematical Statistics 38, 1148–1151.
[6] Sugiyama, T. (1970). Joint distribution of the extreme roots of a covariance matrix. Annalsof Mathematical Statistics. 41, 655–657.
[7] Uhlig, H. (1994). On singular Wishart and singular multivariate beta distributions. Annalsof Statistics 22, 395–405.

ロジスティック回帰モデルにおける安定な非凸スパース正則化法
東京大，国立精神・神経医療研究センター
東京大

小川光紀

東京大

大庭幸治

政策研究大学院大，東京女子医大
東京大

塘由惟

片井みゆき

松山裕

はじめに

ロジスティック回帰モデルでは分離とよばれる状況で最尤推定量が存在しないことが知られて

おり，医学研究などの実証研究においてしばしば問題となる [1, 2]．分離の状況においても推定値が得られるような推定法として Jeffreys 事前分布による罰則項を対数尤度に課す推定法が提案されており，本稿ではこの推定法を Firth の罰則付き最尤法とよぶ [3, 4, 5]．

分離の状況における推定法は検討されている一方で，モデル選択の方法はほとんど提案されていない．本研究では Lasso 法に代表されるスパース正則化法に着目する [6]．Fan と Li は，スパー

ス推定における推定量の望ましい性質としてオラクル性（oracle property）を提案した [7]．推定量
がオラクル性をもつようなスパース正則化法として，2 段階推定によるアプローチである Adaptive

Lasso 法や，非凸の正則化項を用いる SCAD 法，MCP 法などが提案されている [8, 7, 9]．分離の状況においてこれらの手法の利用を検討することもできるが，調整パラメータの選択に伴って推定が不安定となる問題が生じる [10]．また，SCAD 法や MCP 法では絶対値の大きなパラメータ

に対して罰則を課さないことから計算が収束しないことがある．分離の状況においても安定的に推定値を計算することができ，かつ推定量がオラクル性をもつようなスパース推定法は考案されていない．

そこで本研究では，Firth の罰則付き最尤法と非凸スパース正則化法の性質に着目し，Jeffreys 事

前分布に基づく罰則項と非凸スパース正則化項を組み合わせた正則化項を課す推定法を提案する．

提案法

p 次元の説明変数ベクトルを xi = (xi,0 , . . . , xi,p−1 )" ∈ Rp と表す (i = 1, . . . , n)．ただし，切片に対応する成分は xi,0 = 1 である．n × p のデザイン行列を X = (x1 , . . . , xn )" と表す．X は列フ
ルランクであると仮定する．p 次元の回帰係数ベクトルを β = (β0 , . . . , βp−1 )" ∈ Rp と表す．ただし，β0 は切片を表す．標本 (x1 , y1 ), . . . , (xn , yn ) が得られた場合のロジスティック回帰モデルの対#$%!n " "
数尤度を l(β) := i=1
yi xi β − log 1 + exp(x"で表す．また，対角行列 W(β) を W(β) =
i β)
"
%
"
diag {w1 (β), . . . , wn (β)} , wi (β) = πi (β) {1 − πi (β)} , πi (β) = exp(x"i β)/ 1 + exp(xi β) で定1

義する．提案法では罰則付き対数尤度l† (β) := l(β) − n

p−1
!
r=1

pλn (|βr |) +

"
"
1
"
"
log "X" W(β)X"
2の最大化により推定量を定める．pλn (θ) は λn > 0 を用いて θ ∈ [0, ∞) で定義される罰則関数であり，非凸スパース正則化項を構成する．

講演当日は，提案法に関する詳細な条件，提案推定量の存在性と漸近性質に関する議論を紹介したのち，数値実験を通した有限標本下での提案法の性能と実際の医学研究データへの適用例について報告した．

References[1] Mervyn J. Silvapulle. On the existence of maximum likelihood estimators for the binomial
response models. Journal of the Royal Statistical Society. Series B (Methodological), Vol. 43,
No. 3, pp. 310–313, 1981.
[2] Adelin Albert. and J. A. Anderson. On the existence of maximum likelihood estimates inlogistic regression models. Biometrika, Vol. 71, No. 1, pp. 1–10, 1984.
[3] David Firth. Bias reduction of maximum likelihood estimates. Biometrika, Vol. 80, No. 1,
pp. 27–38, 1993.
[4] Georg Heinze and Michael Schemper. A solution to the problem of separation in logistic
regression. Statistics in Medicine, Vol. 21, No. 16, pp. 2409–2419, 2002.
[5] Ioannis Kosmidis and David Firth.

Jeffreys-prior penalty, finiteness and shrinkage inbinomial-response generalized linear models. Biometrika, Vol. 108, No. 1, pp. 71–82, 2020.
[6] Robert Tibshirani. Regression shrinkage and selection via the lasso. Journal of the Royal
Statistical Society. Series B (Methodological), Vol. 58, No. 1, pp. 267–288, 1996.
[7] Jianqing Fan and Runze Li. Variable selection via nonconcave penalized likelihood and its
oracle properties. Journal of the American Statistical Association, Vol. 96, pp. 1348–1360,
2001.
[8] Hui Zou. The adaptive lasso and its oracle properties. Journal of the American Statistical
Association, Vol. 101, No. 476, pp. 1418–1429, 2006.
[9] Cun Hui Zhang. Nearly unbiased variable selection under minimax concave penalty. The
Annals of Statistics, Vol. 38, No. 2, pp. 894–942, 2010.
[10] Mohammad Ali Mansournia, Angelika Geroldinger, Sander Greenland, and Georg Heinze.
Separation in logistic regression - causes, consequences, and control. American journal of
epidemiology, Vol. 187, , 08 2017.

2

ブリッジ推定量を用いた BIC の妥当性とその周辺
高崎経済大学・経済学部

宮田 庸一

科研費シンポジウム（東京理科大学野田キャンパス）10 月 27 日 ∼10 月 29 日
" 記号は, 行列の転置を表すものとする. 観測ベクトル (Yi , x!i )! は, 以下のスパースな線形モデルに従うとする.

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

(j = 1, ..., kn ), および βj∗ = 0 (j = kn + 1, .., pn ) であるとする. ui は期待値 E(ui ) = 0, 分散

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

Yi = β0 + βj1 xj1 i + · · · + βjk xjk i + # i ,(i = 1, ...n),

ただし j1 , ..., jk ∈ {1, ..., pn }, k = 1, ..., pn とする. このとき，どの説明変数をモデルに組み込

むかについては, pn の次元がそれほど大きくない場合には, AIC, BIC などの情報量規準を用いて総当たり法を行うのが標準的なアプローチであろう. しかしその一方で, pn の次元が高くなるにと, 評価すべきモデルの数が指数的に増加するため, 総当たり法を行うのが困難になる. このような問題に対する一つのアプローチとしては, Tibshirani (1996) による LASSO(Least absoluteshrinkage and selection operator) に代表される罰則付き最小二乗推定量がある. 具体的には,!
"
y = (Y1 , ..., Yn )! , X1 = (xji ):n × pn 行列, 1n = (1, ..., 1)! , X = 1n
L0,n (β) = &y − Xβ&2 + λn &β&1を最小にする β の値 β̂L を求める. ただし &β&1 =X1 とおき, 損失関数

(1)j=1 |βj | とする. この推定量は LASSO 推定#pn量と呼ばれ, λn → ∞ (n → ∞) とその他のいくつかの条件のもとで漸近正規性を持つことが知られている. 一方で LASSO 推定量をモデル選択手法の観点から見たときには, いくつかの問題が生じる. それは収束レートが√n である漸近正規性が成り立つためのチューニングパラメーター λnのオーダーの仮定の下では, 真のモデルを選ぶ確率は漸近的に 1 にならないことが知られている.
この問題を改善するために, 式 (1) における罰則項を, 非凸の形のものに置き換えた以下の損失関数を考え, これを最小にする推定量 β̂ を求める:

2

Ln (β) = "y − Xβ" + λn

pn
!
j=1

|βj |γ ,(2)

ただし 0 < γ < 1 とする. この推定量はブリッジ推定量と呼ばれている. 当然のことながら, ブ
リッジ推定においても適切な λn の選び方が重要になる. Huang et al. (2008) では, オラクル性と呼ばれる, 推定量 β̂ が漸近正規性と, ある種の漸近有効性を持ち, なおかつ説明変数の選択に関しても一致性を持つような条件を与えている. 例えば γ = 1/2 として, kn = k(定数), pn = log n とおいた場合, λn のオーダーに関する条件は, λn /(n1/4 (log n)3/4 ) → ∞ (n → ∞), λn = o(n1/2 ) となるが, そのような λn の選択肢は, λn = 3n3/8 でも良いし, λn = 10n3/8 でもよいことになる. すなわち, このような漸近理論に関する結果は λn の選択に関しては, 何も述べていないことになる.

Huang et al. (2009), Wang et al. (2010) では, λn を選ぶために, ある BIC 型の規準を提案しているが, この導出に関する妥当性に関しては, 現時点では明らかになっていない. 本報告では, ある疑似周辺尤度が Yamanishi (1998) により提案された Extended stochastic complexity とみなせることを紹介し, また上述した BIC が，この疑似周辺尤度に対するラプラス近似に対数を取ったものの主要項として正当化できることを報告した.

参考文献
Huang, J., J. L. Horowitz, and S. Ma (2008). Asymptotic properties of bridge estimators insparse high-dimensional regression models. Ann. Statist. 36 (2), 587–613.
Huang, J., S. Ma, H. Xie, and C.-H. Zhang (2009). A group bridge approach for variableselection. Biometrika 96 (2), 339–355.
Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. J. Roy. Statist. Soc.
Ser. B 58 (1), 267–288.
Wang, M., L. Song, and X. Wang (2010). Bridge estimation for generalized linear models witha diverging number of parameters. Statistics & probability letters 80 (21-22), 1584–1596.
Yamanishi, K. (1998). A decision-theoretic extension of stochastic complexity and its applications to learning. IEEE Trans. Inform. Theory 44 (4), 1424–1439.

ロジスティック分布における母数推定について
作村建紀 ∗

柳本武美 †の事後平均は，

1 はじめにµτ"

ロジスティック分布の密度は次で与えられる．

p(x | µ, τ ) =τ exp(τ (x − µ))2.

{1 + exp(τ (x − µ))}τ̂= E[τ µ; πR (µ, τ |x)](2)

= E[τ ; πR (µ, τ |x)]で得られる．このとき，µ の推定量は，(1)

ここで，−∞ < x < ∞, −∞ < µ < ∞, τ > 0 であ

µ̂ =

τ"
µ
,
τ&

(3)

る．この分布の詳しい解説は，[Balakrishnan, 1992]

として求める．

その分布のシンプルな構造にも関わらず，推定量をシド (&µmode , τ&mode ) を考えることもできる．n = 2 のにある．この分布のパラメータ (µ, τ ) の推定では，

また，参照事前分布 πR (µ, τ ) のもとでの事後モー

ンプルな形式で表せないことで知られている．標本とき，µ̂mode = x̄ が成り立つ．また，πR (µ, τ |x) =x = (x1 , . . . , xn ) をこのロジスティック分布からの!n
独立同一な無作為標本とする．また，x̄ = n1i=1 xiτ n−1exp(τ (xi −µ))であるから，尤度を改善し{1+exp(τ (xi −µ))}2た推定量という見方もできる．

とする．E[p; q] は分布 q のもとでの p についての期

2.1

待値を表すとする．

n = 2 における陽な表現

n = 2 のときには，提案推定量は陽な表現が可能で本研究の目的は，パラメータ (µ, τ ) のベイズ推定ある．すなわち，x = (x1 , x2 ) として，式 (2) より，量を導入することと，その性能を評価することである．この目的のために，無情報事前分布として参照τ"µ=事前分布を仮定する [Bernardo, 1979]．

2x̄
2, τ̂ =|x1 − x2 ||x1 − x2 |(4)が成り立つ．さらに，式 (3), (4) より，µ̂ = x̄ が2 提案推定量得られる．これは µ に関して不偏である．また，

1/τ̂ = |x1 − x2 |/2 であるから，E[1/τ̂ ] = 1/τ が本研究でのベイズ的アプローチでは，(τ µ, τ ) の事成り立つ．よって，n = 2 においては，(µ̂, 1/τ̂ ) は

後平均 ("
τ µ, τ̂ ) を扱う．正規分布の場合，この推定量(µ, 1/τ ) の不偏推定量になる．

をプラグインした予測分布はカルバックライブラー

損失のもとで最適な予測子と一致する．このとき，密

2.2

線形性

a, b を任意の定数，1 = (1, . . . , 1)" として，y =度 p(x|µ, τ ) = (2π)−1/2 τ 1/2 exp(−τ (x − µ)2 /2) として，("
τ µ, τ̂ ) = (τ̂ x̄, (n − 1)/'ax + b1 (a $= 0) とすると，

2i=1 (xi − x̄) ) とな

!nり，また，µ̂ = x̄ を得る．この結果から，(µ̂, 1/τ̂ ) はµ̂(y) = aµ̂(x) + b, τ̂ (y) =(µ, 1/τ ) の不偏推定量になる．ロジスティック分布1
τ̂ (x)
a(5)

においてパラメータ (τ µ, τ ) を推定対象とすることが成り立つ．これは積分での変数変換によって証用できると思われる．事前分布は，πR (µ, τ ) ∝ 1/τ

ラメータ µ について左右対称であることと，等式明できる．また，ロジスティック分布の密度はパについての理論的な裏付けはないが，形式的には応(5) から，E[µ̂(x); p(x|0, τ )] = 0 を得る．すると，で与えられる [Ghosh, 2011, Example 1]．

E[µ̂(x); p(x|µ, τ )] = µ が成り立つ．つまり，µ̂ はす

事 前 分 布 πR (µ, τ ) の 下 で の 事 後 分 布 をべての n について µ の不偏推定量である．

πR (µ, τ |x) ∝ p(x|µ, τ )πR (µ, τ ) とすれば，(τ µ, τ )3 リスク比較∗ 法政大学理工学部：〒184-8584 東京都小金井市梶野町

本節では，提案推定量の性能を調査する．比較す3-7-2.

る既存推定量として，最尤推定量 (µ̂ml , τ̂ml ) を考え† 統計数理研究所：〒190-8862 東京都立川市緑町 10-3.

1

結果となったが，しかしこれらのリスク値の差は小表 1: 1/τ̌ のバイアスn1/τ̂

1/τ̂mode

1/τ̂ml

2
3
4
5
10
20

(0)
0.00340.00300.00180.00490.0017−0.03020.01800.0241
0.02140.01620.0076−0.3442−0.2109−0.1529−0.1211−0.0549−0.0277

さかった．

3.4

パ ラ メ ー タ (µ, τ ) の リ ス ク 比 較 を カ ル バ ッ ク

ラ イ ブ ラ ー 損 失 の も と で 行 っ た ．p = p(x|µ, τ ),p̌ = p(x|µ̌, τ̌ ) として，カルバックライブラー損失
E [log (p̌/p) ; p̌] においては，n = 2 で提案推定量が良い結果であり，それ以外は提案推定量と事後モードがほぼ同等の性能であった．また，E[log (p/p̌) ; p]

では n = 2, 3 で提案推定量が優勢な結果を示していた [Johnson et al., 1995]．性能はいくつかの損失関1数のもとで，その平均 m(µ, τ ) についての比較た．一方，(µ̂ml , τ̂ml ) は値が大きい結果となった．

i=1 L(θ̌i , θ) をリスクとし!mて評価した．ここで，m はシミュレーション回数，θ̌i

4 考察

定値，θ は真値，L(θ̌i , θ) は損失関数である．すべて

本研究で提案したベイズ推定量の性能を見る上で，回行った．(µ, τ ) をパラメータの真値，(µ̌, τ̌ ) をそ考える上でその根拠の一つである正規分布においては i 番目のサンプルから得られるパラメータ θ の推

ロジスティック分布のパラメータ推定において，のケースにおいて，シミュレーションは m = 10, 000

τ よりもまず 1/τ に注目した．これは，推定対象をの推定値とするが，パラメータ (µ, τ ) の各コンポーは不偏であることから，ロジスティック分布におい

ネント µ と τ についても個別に比較した．式 (5) かてもその不偏性を期待したことが理由である．実際，ら，µ = 0, τ = 1 のケースのみを扱った．

3.1

シミュレーションの結果から，1/τ についてはほぼ

1/τ についての比較

不偏であることが観測された．ベイズ推定量では，

まず，1/τ̌ に注目し，そのバイアス，二乗損失，変一般に不偏性を考慮していない．本研究で提案した

動係数についてリスク比較した．n = 2 のとき 1/τ̂

推定量も不偏性から誘導されたものではなく，指数は式 (4) で計算される．表の丸括弧の表記は理論値型分布族における最適性を有する予測子に着想を得の結果であることを表している．提案推定量の 1/τ̂て考えたものである．それにも関わらず，近似的不はほぼ不偏であることが観察された．1/τ̂mode もほ

偏性が観測されることは興味深い．この非指数型分

ぼ不偏の傾向があるが，1/τ̂ のほうがより不偏の傾布族においての推定対象の設定と，その不偏性につ向が強い．ここで，τ̂ は式 (2) の複雑な形で計算さいての理論的根拠については，今後明らかにされるれることに注意されたい．1/τ̂ は近似的不偏性を有ことを期待したい．

しており，また変動係数の値も小さいことが観察さ

参考文献

れた．

3.2

τ についての比較[Balakrishnan, 1992] Balakrishnan, N. (1992). Handbook of the Logistic Distribution. CRC Press.
[Bernardo, 1979] Bernardo, J. M. (1979).
Reference posterior distributions for bayesian inference.
Journal of the Royal Statistical Society. Series B(Methodological), 41(2):113–147.
[Ghosh, 2011] Ghosh, M. (2011). Objective Priors:
An Introduction for Frequentists. Statistical Science, 26(2):187 – 202.
[Johnson et al., 1995] Johnson, N. L., Kotz, S., andBalakrishnan, N. (1995). Continuous Univariate
Distributions, Volume 2, volume 289. John wiley& sons.

次に，τ について調べた．τ についての提案推定量の性能は，バイアス，二乗損失，変動係数に加えて，
ガンマ分布から得られる近似カルバックライブラー
損失によって評価した．二乗損失と双対な近似 KL
損失では n = 2, 3 で良い性能を示しており，また近似 KL 損失では選択したすべての n において良好な性能を示していた．

3.3µ についての比較

次に，µ の推定量についてのリスク比較を行った．

まず，すべての n のおいて，提案推定量 µ̂ は不偏である．二乗損失で比較すると，n = 3, 4, 5 で µ̂ml のリスクが小さく，n = 10, 20 で µ̂ のリスクが小さい2対数オッズ比の推定量と検定を改善するためのベイズ法
三重大学
統計数理研究所

小椋 透
柳本 武美

1. はじめに独立した 2 群間における 2 値データは表 1 に示す 2 × 2 分割表に要約される. 臨床試験においては, 2 群間における予後の違いを評価する指標の一つとして対数オッズ比がしばしば用いられる. 対数オッズ比は最尤法 (MLE) で推定されることが多いが, 一つ以上のセルに 0 がある場合
は −∞, ∞ 又は推定不能となる. これは明らかに過小推定又は過大推定である. 条件付き MLE(CMLE) は, さまざまなモデルで MLE よりも優れていることが報告 (Yanagimoto and Anraku,
1989; Yanagimoto, 1991) されているが, 対数オッズ比における CMLE は MLE と同様に一つ以上のセルに 0 がある場合は −∞, ∞ 又は推定不能となる. 臨床試験ではサンプルサイズが小さい場合や稀な発症確率の場合に一つ以上のセルに 0 があることは珍しくないが, その場合においてもリーズナブルな対数オッズ比の推定が必要とされる. 研究者が対数オッズ比を使用する場合, 群ごとの割合や群ごとのオッズに対して関心が低いため, 群ごとの割合や群ごとのオッズを改善する必要がない場合がある. 本研究はベイズ法を用いて対数オッズ比を直接推定する方法を提案した.
観測値
処理群
対照群

Positivexy

表 1. 2 × 2 分割表
Negative Total
確率
n−x
n
処理群
m−y
m
対照群

Positivepq

Negative1−p1−q2. 提案推定量
二項分布 Bi(n, p) 及び Bi(m, q) からの観測値をそれぞれ x 及び y とした. ここで, n と m はサン
プルサイズとし, p と q は確率とした. 研究者が Fisher’s exact test を用いる場合は p と q の推定にも関心があると考えられる. 一方, 研究者が対数オッズ比を用いる場合は p と q の推定に関心が低い
p(1−q)qと考えられる. そこで, 本研究は二つのパラメータ p と q の代わりに α = log 1−qと β = log q(1−p)
を用いてベイズ推定を行った. このとき, β は対数オッズ比になっており, β の推定は対数オッズ比の直接推定であった. ベイズ法により, Jeffreys 事前分布を用いた事後密度はπM (α, β|x, y) ∝ !min(n,t)

exp(xβ)"n#" m #z=max(0,t−m) z×t−zexp((t + 1)α + β/2)
(1
+
exp(α
+ β))n+1 (1 + exp(α))m+1
exp(zβ)(1 + exp (α(t, β) + β))n (1 + exp (α(t, β)))m
exp(tα(t, β))として表された. ただし, t = x + y であり, α(t, β) は t = n 1+exp(α(t,β)+β) + m 1+exp(α(t,β)) をexp(α(t,β)+β)exp(α(t,β))解くことによって求められた. よって, 二つのパラメータ α と β は (α̂, β̂) = −∞ −∞ (α, β) ×πM (α, β|x, y)dαdβ として推定された.
提案推定量に対応してオッズ比の信用区間に基づいた検定は次のように構成された. 帰無仮説
H0 : β = 0, 対立仮説 H1 : β > 0 の片側検定とするとき, H1 の credibility (Ghosh et al., 2006) は
$∞ $∞

$∞$∞

πM (α, β|x, y)dαdβ−∞ −∞ πM (α, β|x, y)dαdβ0$−∞C(x, y) = Pr{β > 0; πM (α, β|x, y)} = $ ∞∞として表されて, 検定統計量を T (x, y) = C(x, y) とした. 次に, 領域 S(x, y) = {(w, z)|T (w, z) ≥T (x, y)} を設定すると, P 値は次の式で算出された.

Pr{S(x, y)|H0 } = max{S(x, y); p(x, y|α, 0)} = max
α

α

1

%

(w,z)∈S(x,y)& '& 'nw

m exp((w + z)α)z (1 + exp(α))N3. リスク比較
提案推定量と既存の推定量の比較に次のロスが用いられた. 一つ目のロスは kullback-Leibler
divergence (KLD) に基づいた e-divergence 損失でLe ((p̌, q̌), (p, q)) = np̌ logp̌
1 − p̌q̌
1 − q̌
+ n(1 − p̌) log+ mq̌ log + m(1 − q̌) log
p
1−p
q
1−q

として表された. ここで, p̌ と q̌ は推定量を表した. 二つ目のロスは Le の双対である m-divergence
損失 (dual KLD; dKLD) で Lm ((p, q), (p̌, q̌)) = Le ((p̌, q̌), (p, q)) として表された. 三つ目のロスは平均二乗誤差 (mean squared error; MSE) で LMSE (β̌, β) = (β̌ − β)2 として表された. ここで, β̌ は推定量を表した. 様々な設定でリスク比較を行い, 提案推定量は既存の推定量 (CMLE,
MLE, Jeffreys 事前分布を仮定した自然パラメータの事後平均 (posterior mean of the canonicalparameter; PMCP), 修正 MLE (MMLE), 修正中央値の不偏推定量 (MMUE)) と比べて, 広い範
囲でリスクが低いことが確認された.

4. 実例
進行性大腸癌患者における二つの新しい化学療法 (ホモハリングトニンとカラセミド) のランダ
ム化された第 II 相臨床試験のデータ (Parzen, 2002) を使用した. この研究の結果は 2 × 2 分割表
として表 2 のように要約された. 提案推定量と既存の推定量で対数オッズ比の推定を行った. 検定の比較には, exact 法が用いられた. 結果は表 3 にまとめられた.
表 2. 第 II 相臨床試験のデータにおける 2 × 2 分割表
主要評価
副次評価
Positive Negative Total
Positive Negative

ホモハリングトニン群
カラセミド群

主要評価
副次評価

0
0

14
11

14
11

0.49110.5108−0.2318
0.3365

2
1

12
10

Total
14
11

表 3. 第 II 相臨床試験のデータにおける対数オッズ比の推定および検定
推定
検定 P 値提案法
PMCP CMLE
MLE
MMLE MMUE
提案法
Exact

−0.2458
0.4780−0.24100.4845−0.24090.4271

0.999
0.272

> 0.999
0.593

5. まとめ臨床試験においてサンプルサイズが小さい場合や稀な発症確率の場合には, 実例のような偏りが生じることは珍しくなかった. その場合に, よく知られている MLE と CMLE による推定は −∞,∞ 又は算出不能となり, これは過小推定又は過大推定であった. 対数オッズ比を用いる場合, 群ごとの割合の改善を必要としないことが多いことから, 本研究は対数オッズ比を直接推定する方法を提案した. リスク比較において, 提案推定量のリスクは既存の推定量より安定して小さいことが確認された. 実例においても, 提案推定量の有効性が検証された.
参考文献
Gart, J. J. and Zweifel, J. R. (1967). Biometrika, 54(1/2), 181–187.
Ghosh, J. K., Delampady, M. and Samanta, T. (2006). An introduction to Bayesian analysis:
theory and methods. New York: Springer.
Haldane, J. B. (1956). Ann Hum Genet, 20(4), 309–311.
Parzen, M., Lipsitz, S., Ibrahim, J. and Klar, N. (2002). J Comput Graph Stat, 11(2), 420–436.
Yanagimoto, T. (1987). Ann I Stat Math, 39(2), 247–261.
Yanagimoto, T. and Anraku, K. (1989). Ann I Stat Math, 41(2), 269–278.
Yanagimoto, T. (1991). Ann I Stat Math, 43(4), 735–746.

2‫ڞ‬໾ղੳͷ࠶ߏ੒ͱ֦ுͷࢼΈ༄ ຊ ෢ ඒ: ౷‫਺ܭ‬ཧ‫ॴڀݚ‬1. ং‫ڞ‬໾ղੳ͸ඪຊ෼෍͕ࢦ਺෼෍ʹै͏৔߹ʹɺཧղ͕༰қͳࣄલ෼෍ΛԾఆͯ͠ɺ؆୯ͳઢ‫ܗ‬ਪఆ஋͕ಘΒΕΔਪଌͷ࿮૊Έͱͯ͠ొ৔ͨ͠ɻBernardo and Smith (2000) Ͱ΋‫ج‬ຊత࿮૊Έͱͯ͠঺հ͞Ε͍ͯΔɻ͔͠͠ɺࠓ೔ͷࢹ఺Ͱ‫ݟ‬ΔͱԾఆ͕ແࡶ࡞Ͱ͋ΔҰํͰɺ‫ࣜܗ‬తͰൃలੑʹ๡͍͠࿮૊ΈͰ͋Δɻͦ͜Ͱɺ֦ுΛҙࣝͨ͠‫ڞ‬໾ղੳͷ࠶ߏ੒ʹऔΓ૊ΜͰ͍Δɻຊ‫ڀݚ‬ͷ΋͏Ұͭͷૂ͍͸ɺ‫ڞ‬໾෼෍ͷ‫ͯ͠ͱݶۃ‬ͷແ৘ใࣄલ෼෍Ͱ͋Δɻແ৘ใͱ͸ར༻Ͱ
͖Δ৘ใ͕‫ݶ‬Γͳ͘খ͍͞ɺͱཧղͰ͖Δɻͦ͏͢Δͱɺແ৘ใࣄલ෼෍ΛҰൠԽͯ͠޿ཱ͍৔͔
Β‫͢ڀݚ‬Δ͜ͱʹͳΔɻ‫ڵ‬ຯਂ͍݁ՌΛಘͭͭ͋ΔͷͰใࠂͨ͠ɻ

2. ‫ط‬ଘͷ‫ڞ‬໾ղੳ!
"ࢦ਺෼෍଒ʹଐ͢Δେ͖͞ n ͷඪຊαΠζͷඪຊͷີ౓Λ p(x|θ) = exp{n x̄θ − M (θ) } a(x), θ ∈Θ. ͱ͢Δɻ‫ڞ‬໾ղੳͰ͸ɺඪຊີ౓ʹରͯ͠૒ରͳߏ଄Λ࣋ͭ‫ڞ‬໾ࣄલີ౓Λ!
"
πT (θ; m, δ) = exp{δ mθ − M (θ) }h(m, δ)(1)

ͱఆٛ͢ΔɻओཁͳҼࢠ͸ɺx̄, n ͕ m, δ ͱஔ͖‫͑׵‬ΒΕ͍ͯΔɻ‫ڞ‬໾ղੳͰ͸ m ͱ δ > 0 ͷ஋͸ओ‫؍‬తʹ༩͑Δ͜ͱ͕‫ج‬ຊͰ͋Δɻࣄ‫౓ີޙ‬͸ɺ µ∗ = (nx̄ + δm)/(n + δ) ͱ͓͍ͯ!
"πT (θ| x) = exp{(n + δ) µ∗ θ − M (θ) }h(µ∗ , n + δ)

(2)ͱද͞ΕΔɻ͜ͷ෼෍͸ࣄલີ౓ͱಉ͡‫ܗ‬Λͯ͠ɺµ, δ Λ µ∗ , n + δ ͱೖΕସ͑ͯಘΒΕΔɻ͜ͷੑ࣭Λ closed under sampling ͱ‫ݺ‬ΜͰ͍Δɻ฼਺ µ = M " (θ) ͸ x̄ ͷฏ‫͋Ͱۉ‬Δɻ͜ͷ฼਺͸ Argmaxµ π(θ| x) = µ∗ ɺ͋Δ͍͸

Epost [µ] = µ∗(3)Ͱ༩͑ΒΕΔɻ‫ऀޙ‬ͷද‫͕ݱ‬੒ཱ͢ΔͨΊͷे෼৚݅͸ɺ฼਺ۭؒ Θ ͕։ू߹Ͱ͋Δ͜ͱͰ͋Δ(Diaconis and Ylvisker, 1979)ɻ͜ͷਪఆ஋͸؆୯ͳ‫ܗ‬Λ͍ͯ͠Δ্ʹɺ n ͱ δ ͰॏΈ෇͚ͨͳ͍จలͰ͋Δɻղऍ͕༰қͰ͋Δɻ

3. ‫ڞ‬໾ղੳ΁ͷջٙ‫ط‬ଘͷ‫ڞ‬໾ղੳΛ࢔ࡉʹ఺‫͢ݕ‬Δͱೋͭͷॏཁͳܽ఺ʹ‫͘෇ؾ‬ɻ1) ࣄલີ౓ͷબ୒͕҆қͰ͋Δɻ
ࣗવ฼਺্ͷҰ༷෼෍͕ӄʹԾఆ͞Ε͍ͯΔ2) ଟ࣍‫ݩ‬฼਺Ͱࣄલີ౓ΛԾఆ͢ΔͨΊʹ͸େ෯ͳमਖ਼͕ඞཁʹͳΔɻਖ਼‫ن‬෼෍ N(µ, σ 2 ) ͷ৔߹΋ѻ͍ʹ͍͘ߋʹ͸ɺ࣍ͷೋ఺Λ΋ࢦఠ͢Δɻ3) ࣄલີ౓ͷબ୒Ͱ͸ࣗવ฼਺͕ɺਪఆͰ͸‫ظ‬଴஋฼਺͕༻͍ΒΕ͍ͯΔɻԿΒ͔ͷઆ໌͕‫ٻ‬ΊΒΕΔ4) ࣄલີ౓ͷ଒ͷ௒฼਺ (µ, δ) ͸௚ަ͢ΔΑ͏ʹબͼ͍ͨɻ֓೦্શ͘ҟͳΔ໾ׂΛՌͨ͢௒฼਺Ͱ͋Δ‫ط‬ଘͷ‫ڞ‬໾ղੳ͸ɺࣄ‫ͱ౓ີޙ‬ਪఆ஋ͷ؆ૉͳ‫͕ܗ‬ॏࢹ͞Εͯɺࣄલີ౓ͷબ୒͸͓͟ͳΓͰ͋Δɻ͜ͷ‫ݱ‬ঢ়͸ϕΠζਪ࿦ͷ‫ج‬ຊతͳ࿮૊Έͱဃ཭͍ͯ͠ΔɻϕΠζ๏͸ɺࣄલ෼෍Λ஫ҙਂ͘બ΂͹ɺ‫ޙ‬ͷਪ࿦͸‫ۃ‬Ίͯ؆ศͰ͋Δɻ

4. ‫ڞ‬໾ղੳͷ࠶ߏ੒‫ط‬ଘͷ‫ڞ‬໾ղੳͰ͸ɺࣄલີ౓ πT (θ; m, δ) ʹඪຊີ౓ͷ a(x) ʹରԠ͢Δ୆ଌ౓ b(θ) ͕ͳ͍ɻ‫ݴ‬͍‫͑׵‬Δͱɺb(θ) ∝ 1 ͕ӄʹԾఆ͞Ε͍ͯΔɻࣄલີ౓ͷબ୒Ͱ͸‫ࣜܗ‬తʹ δ = 0 ͱ͓͘ͱɺີ౓
Ͱ͸ͳ͘ͳΔ͜ͱ͕ଟ͍ɻ͔͠͠ɺࣄ‫Ͱ౓ີޙ‬͸‫ࣜܗ‬తʹ δ = 0 ͱ͓͍ͯ΋௨ৗ͸ີ౓ʹͳΔɻ͜ͷࣄ࣮͸٬‫؍‬తϕΠζ๏ͷ‫ج‬ຊతͳ‫ڌ‬ΓॴͰ͋Δɻࣄલີ౓ͷԾఆͰ҉໧ͷԾఆΛઃ͚Δͷ͸ॏେͳܽ఺Ͱ͋Δɻ͔͠΋ɺ۩ମతͳྫΛௐ΂Δͱɺ͜ͷԾఆ͸௨ৗશ͘ັྗతͰͳ͍ɻ
ߋʹɺඪຊ෼෍͕ࢦ਺෼෍଒ʹଐ͢Δͱ͖ͷࣄલີ౓ͷද‫Ͱݱ‬͸ M (θ) ʹରͯ͠‫ڞ‬໾ͳತؔ਺
N (m) Λಋೖ͢Δɻ݁‫ࣄ͍͠৽ہ‬લີ౓͸!
"π(θ; m, δ) = exp{δ mθ − M (θ) − N (m) }b(θ) exp(−k(m, δ)

(4)

!
"
ͱද͞ΕΔɻલऀͷ exponent ͸ Kullback-Leibler divergence Λ༻͍ͯ D p(x|N ! (m)), p(x|θ) ͱ΋ද͞ΕΔɻଈͪ θ = N ! (m) ͷ࣌ʹ 0 ʹͳΔɻ͜ͷࣄ࣮͸஋ m ͷผͷઆ໌Λ༩͑Δɻ୆ଌ౓ b(θ) ͸໌ࣔతʹԾఆ͢Δɻ͜ͷଌ౓ͷબ୒ͱ٬‫؍‬తϕΠζ๏ͷ‫͢ڀݚ͕ऀڀݚ‬Δແ৘ใࣄલ෼෍ͱ͸ຆͲಉ͡Ͱ͋ΔɻDefault ͱͯ͠ͷબ୒͸ Jeffreys prior Ͱ͋ΔͱࢥΘΕΔɻෳ਺࣍‫ݩ‬ͷ฼਺৔߹ʹ͸ɺBernardo (1979) ͕ఏএͨ͠ reference prior ͷఆٛ͸આಘྗ͕͋ΓୈҰબ୒Ͱ͋Δ
ͱߟ͑ΒΕΔɻࣄ‫౓ີޙ‬͸ (2) ͱಉ༷ʹಘΒΕΔɻ͜ͷఆٛͰ΋ੑ࣭ closed under sampling ͸੒Γཱͭɻ‫ط‬ଘͷࣄ‫ޙ‬Ϟʔυɺ‫ظ‬଴஋฼਺ͷࣄ‫ޙ‬ฏ‫ؔ͢ʹۉ‬Δੑ࣭͸੒Γཱͨͳ͍ɻͦΕͰ΋ɺਪఆ஋͸ µ∗ Λ௨ͯ͠ͷΈ x̄, m ɹʹґଘ͢Δɻ
ਪఆ஋͸ࣗવ฼਺ͷࣄ‫ޙ‬ฏ‫͕ۉ‬ਪ঑͞ΕΔɻ࠷େͷཧ༝͸ඪຊ෼෍͕ࢦ਺෼෍଒ʹଐ͢Δ৔߹ʹ͸plug-in ༧ଌࢠ p(y|θ̂) ͕޿͍ൣғͰ࠷దੑΛຬͨ͢͜ͱͰ͋Δ (Yanagimoto and Ohnishi, 2009)ɻ
ฏ‫ۉ‬஋฼਺ͷࣄ‫ޙ‬ฏ‫Ͱۉ‬͸ͦͷΑ͏ͳੑ࣭͸ຬͨ͞ͳ͍ɻ
·ͨࣗવ฼਺ͷࣄ‫ޙ‬ฏ‫ۉ‬͸
Epost [θ] = N ! (µ∗ ) + ∂k(µ∗ , n + δ)/∂mͱද‫͖Ͱݱ‬Δɻ͜ͷੑ࣭͸‫ۃ‬Ίͯऑ͍৚݅ͷԼʹ੒ཱ͢Δɻ͜ͷࣄ࣮͸ද‫( ݱ‬3) ͕‫ٯ‬Ψ΢ε෼෍ɺ‫ٯ‬ೋ߲෼෍ɺ von-Mises ෼෍Ͱ੒Γཱͨͳ͍͜ͱͱରൺ͞ΕΔɻ۩ମతʹ͸ɺb(θ) ͷબ୒Ͱ͸ɺ‫ن‬֨Խఆ਺͕ δ ͷΈͷؔ਺ͱͳΔ̺ͱ͕Ұͭͷ‫ن‬४ʹͳΔɻ͜ͷ‫ن‬४͕ຬͨ͞ΕΔͱɺm ͱ δ ͕௚
ަ͢Δͱ‫ʹڞ‬ɺθ ͷࣄ‫ޙ‬ฏ‫͔ۉ‬Β༠ಋ͞ΕΔ µ ͷਪఆ஋͕ µ∗ ʹͳΔɻจ‫ ݙ‬1) Bernardo, J.M. (1979). J. Roy. Statist. Soc.. Ser. B, 41, 113-147. 2) Bernardo, J.M.,
Smith, A.F.M., (2000). Bayesian Theory. Wiley: Chichester. 3) Diaconis, P and Ylvisker,
D. (1979). Ann. Statist., 7, 269-281. 4) Yanagimoto, T. and Ohnishi, T., (2009). J. Statist.
Plann.Inf., 139, 3064-3075.

大規模時空間データに対するベイズモデル
東京大学経済学研究科 若山智哉

1

はじめに

携帯端末の普及に伴い膨大な人口データが取得され、その活用が求められている。少なくとも東

京都では、過去のある時間、ある場所にどれだけの人数が滞在しているかが分かるようになっている。このようなデータは、実に様々な場面で役に立つ。例えば、交通計画による混雑や交通渋滞の緩和、タクシーや宅配便の効率化、区画ごとの消費促進、災害時の避難誘導あるいは犠牲者数推定などに活用することができる。本報告ではこうしたデータから知見を抽出し、将来の意思決定につながるモデルの開発を紹介した。

2

提案モデル

都市の人口データの分析では、都市の構造を考えてモデルを構築する必要があった。例えば住宅

街では時間ごとの滞在人口の変動は小さいが、オフィス街では日中は大きく増加するが、夜中や朝は住宅街の半分以下の人数になることもある。これらの傾向は近隣の地区にも波及し、結果として全体の人口分布図に表れる。つまり、空間依存的なデータである。

2.1

観測誤差モデルyts (τ ) を t ∈ {1, . . . , T } 日に地点 s ∈ {1, . . . , N } で観測されたデータとする。τ ∈ {τ1 , . . . , τK }は観測点であり、毎日全ての地点で K 回観測されるとする。

このデータに対して、次のような観測誤差モデルを考えた。

yts (τ ) = zts (τ ) + εts ,εts ∼ N (0, e2s ),t = 1, . . . , T,s = 1 . . . , N,

ここで、εts は t や s に対して独立な誤差項であり、e2s は地点ごとの未知分散である。また、zts は次のガウス過程に従うと仮定した。

zts ∼ GP(ft , ηs2 R(φs )),!s を (i, j) 成分が
ここで、ft は平均関数、R(φ) = ρφs (d) は共分散関数である。つまり、Rηs2 ρφs (|τi − τj |) な K × K 行列とすると、のように表せる。

"#!s ,(zts (τ1 ), . . . , zts (τK )) ∼ N (fts (τ1 ), . . . , fts (τK )), R

K 次元のベクトルを確率過程のパスの有限部分集合とすることには実用上の利点がある。もし

通常の多変量解析のように K × K 行列をそのまま推定する場合は、K × (K − 1)/2 個のパラメー
タが各地点で必要になってしまう。しかし、確率過程のパスの有限部分集合と見ることで、共分散カーネルのパラメータ（今回のケースでは ηs と φs ）のみを推定すればよく、計算上有益である。

1

2.2

因子モデルの導入

さらに、次の因子モデルを導入した。

!s ))νt ∼ N (0, (Rzt = (B ⊗ IK )xt + νt ,xt = Gxt−1 + ωt ,

ωt ∼ N (0, Λ),

ここで、zt := (zt1 , zt2 , . . . , ztN ) := (zt1 (τ1 ), . . . , zt1 (τK ), zt2 (τ1 ), . . . , ztN (τK )) は N K 次元ベク
トル、νt はその誤差項、xt := (xt1 , xt2 , . . . , xtM ) := (xt1 (τ1 ), . . . , xt1 (τK ), xt2 (τ1 ), . . . , xtM (τK ))は M K 次元ベクトル、B は次の N × M 行列である。(M < N )



1
b21
b31
..
.







B=
 bM 1 bM +1,1


..

.
bN,1

0
1
b32
..
.

0
0
1

bM 2
bM +1,2..
.

bM 3
bM +1,3

bN,2

bN,3

...
...
...
..
.

0
0
0
..
.









.

...
1

. . . bM +1,M 


..
..

.
.
...
bN,M

これは、N 地点のトレンドを、より少数の M 地点のトレンドで説明する因子モデルになっている。
また、因子負荷行列 B について、対角成分に 1 が並ぶようにデザインし識別可能性を保証した。

ここで重要なのは、因子をどのように決定するかという問題であった。一変量問題でも因子自体や因子の数を決めるのは難しく、領域知識に基づいて重要そうなものを揃えるという方法が一般的である (Prado et al., 2021)。この選択の問題点は、必要そうな因子を主観的に選択する必要があり、肝心な因子が選択されているかどうかが不明な点である。そこで我々は縮小事前分布を用いた
因子選択の方法も提案した。また、因子負荷行列に空間依存構造を導入した。

3

数値実験とデータ解析

数値実験を通して、時空間データに対する提案手法の有用性を確かめるために。また、実際に人流データに対して提案手法に休日効果を加えらものを実装し、何が読み取れるかについて議論
した。

参考文献
Prado, R., M. A. Ferreira, and M. West (2021). Time series: Modeling, computation, andinference.

2‫ڀݚ࡯؍‬ͷޮՌਪఆ஋ΛඪతूஂʹҰൠԽʗҠૹ͢Δํ๏ͷ‫ݕ‬౼౦‫ژ‬ཧՊେֶ େֶӃ޻ֶ‫ڀݚ‬Պ ৘ใ޻ֶઐ߈ɹງߐ ༔ੜ౦‫ژ‬ཧՊେֶ ޻ֶ෦ ৘ใ޻ֶՊɹࣰ࡚ ஐେ͋ΔΞ΢τΧϜʹର͢Δ࣏ྍޮՌͷ‫ʹূݕ‬͸ɼ಺తଥ౰ੑ͕‫ظ‬଴͞ΕΔϥϯμϜԽൺֱࢼ‫ٻ͕ݧ‬ΊΒΕΔ͜ͱ͕ଟ͍ɽ͔͠͠ɼϥϯμϜԽΛ൐͏Α͏ͳհೖࢼ‫ݧ‬ͷࢀՃऀ͸‫͍͠ݫ‬ద֨‫ج‬४ʹΑͬͯબ୒͞Εɼ͞Βʹ༷ʑͳཧ༝Ͱ‫ࢀڀݚ‬ՃʹࢸΔͨΊɼ‫ࢪ࣮ྍ࣏ʹ࣮ݱ‬Λ૝ఆ͢Δूஂʢඪతूஂʣͱ༷ʑͳ
఺ͰҟͳΔ͜ͱ͕ଟ͍ɽಛʹɼΞ΢τΧϜͷϦεΫҼࢠʹؔͯ͠‫ࢀڀݚ‬Ճूஂͱඪతूஂ͕ҟͳͬͯ
͍Δͱɼ͍ͣΕ͔ͷޮՌࢦඪͰϦεΫҼࢠʹΑΔ࣏ྍޮՌͷम০͕ੜ͡Δɽͦͷ݁Ռɼࢼ‫ஂूݧ‬શମʹର͢Δ࣏ྍޮՌͷਪఆ஋͕ɼඪతूஂͰͷ࣏ྍޮՌͱେ͖͘ҟͳΔ͜ͱ͕͋Δɽ͜ͷΑ͏ͳঢ়‫گ‬͸ɼ֎తଥ౰ੑ͕ܽ೗͍ͯ͠ΔҰྫͰ͋Δɽۙ೥ɼޮՌͷҰൠԽՄೳੑʢgeneralizabilityʣ͋Δ͍͸ҠૹՄೳੑʢtransportabilityʣͷୡ੒Λ໨ࢦ͢ख๏͕જࡏΞ΢τΧϜͷ࿮૊ΈͰٞ࿦͞Εɼ֎తଥ౰ੑͷͨΊͷҰൠԽʗҠૹΛҙਤ͢Δࡍʹ΋ɼ಺తଥ౰ੑΛߴΊΔ͜ͱΛ໨తͱ͢Δަབྷௐ੔Ͱ༻͍ΒΕΔηϛύϥϝτϦοΫਪఆྔΛద༻Ͱ͖Δ͜ͱͱɼͦͷͨΊͷࣝผԾఆ͕੔ཧ͞Εͭͭ͋ΔʢDahabreh

et al., 2019a, 2019b, 2020; Lesko et al., 2017ʣɽɹ্‫ه‬ͷऔΓ૊Έ͸ɼϥϯμϜԽൺֱࢼ‫ݧ‬ͷ݁ՌΛɼҟͳΔඪతूஂʹҰൠԽʗҠૹ͢Δ࿮૊ΈͰߟ͑ΒΕ͖ͯͨɽ͔͠͠ɼ‫ʹ࣮ݱ‬͸‫ަͰڀݚ࡯؍‬བྷௐ੔Λ௨ͯ͠ಘ࣏ͨྍޮՌͷਪఆ஋ΛɼҟͳΔඪతूஂʹ൓ө͍ͨ͠ঢ়‫گ‬΋ଟ͘ଘࡏ͢ΔɽJASTIS ‫ڀݚ‬ʢࣾձͱ৽‫ܕ‬λόίʹؔ͢ΔΠϯλʔωοτௐࠪ‫ڀݚ‬ʣͰ͸Πϯλʔωοτௐࠪूஂͷ݁ՌΛ೔ຊਓશମूஂͷ݁ՌʹҰൠԽ͢Δ͜ͱΛࢼΈ͍ͯΔ
ʢTabuchi et al., 2021ʣɽ·ͨɼWeb Ξϯέʔτ͔Βੜ‫׆‬श‫ؔ͢ʹ׳‬Δ༗ප཰Λਪఆ͠ɺαϯϓϦϯά
֬཰ͷ‫཰֬ٯ‬Λ༻͍ͯ Web Ξϯέʔτͷର৅ूஂ͔ΒΑΓ޿ൣͳඪతूஂʹ݁ՌͷҰൠԽΛࢼΈͨ‫ݚ‬‫ڀ‬΋͋ΔʢFerri-Garcı́a et al., 2021ʣɽ͜ͷͱ͖ɼϥϯμϜԽൺֱࢼ‫Ͱݧ‬ͷҰൠԽʗҠૹͷํ๏Λɼ‫؍‬࡯‫Ͱڀݚ‬ͷҰൠԽʗҠૹͷํ๏ͱͯͦ͠ͷ··ద༻͢Δ͜ͱ͕Ͱ͖Ε͹؆୯Ͱ͋Δ͕ɼ‫Ͱڀݚ࡯؍‬ͷҰൠԽʗҠૹΛ೉͘͢͠Δٕज़తͳ໰୊͕গͳ͘ͱ΋ 3 ͭଘࡏ͢Δɽ1 ͭ໨͸‫Ͱ಺ڀݚ࡯؍‬ͷަབྷௐ੔ͷඞཁੑɼ2 ͭ໨͸ަབྷௐ੔ʹඞཁͳม਺ͱҰൠԽʗҠૹʹඞཁͳม਺͕ҟͳΓಘΔՄೳੑɼ3 ͭ໨͸‫ڀݚ‬։࢝‫ׂʹޙ‬Γ෇͚ΒΕΔϥϯμϜԽൺֱࢼ‫ͱݧ‬ҟͳΓɼ࣏ྍ͕‫ࢀڀݚ‬ՃʹӨ‫͢ڹ‬ΔՄೳੑͰ͋
Δɽ͜ͷΑ͏ͳٕज़తͳোน͸͜Ε·Ͱͷ‫Ͱڀݚ‬͸શ͘ࢦఠ͞Ε͓ͯΒͣɼํ๏࿦Λ౷Ұతʹ੔ཧͨ͠‫ڀݚ‬΋ͳ͞Ε͍ͯͳ͍ɽ
ɹຊ‫Ͱڀݚ‬͸ɼ͜Ε·ͰϥϯμϜԽൺֱࢼ‫ݧ‬ͷҰൠԽʗҠૹͰ‫ݕ‬౼͞ΕͨਪఆྔͷதͰ΋‫཰֬ٯ‬ॏΈ෇͚ਪఆྔʢinverse probability weightingʀҎԼɼIPW ਪఆྔʣʹண໨͠ɼ্‫ ه‬3 ͭͷٕज़తͳ໰୊
Λߟྀͨ͠‫ڀݚ࡯؍‬ͷ࣏ྍޮՌͷਪఆ஋ΛඪతूஂʹҰൠԽʗҠૹ͢Δ IPW ਪఆྔΛఏҊ͢Δɽ൓ࣄ࣮‫ه‬๏Λ༻͍ͨ୅਺ʹΑ্ͬͯ‫ه‬໰୊఺ͷঢ়‫گ‬Λ੔ཧ͠ɼҰൠԽʗҠૹ͞ΕͨޮՌͷࣝผՄೳ৚݅Λಘͯ IPW ਪఆྔΛಋग़ͨ͠ɽຊൃදͰ͸ɺҠૹʹ͍ͭͯఏҊ͢Δ IPW ਪఆྔͷৄࡉ΍਺஋࣮‫ʹݧ‬ΑΔੑೳධՁͷ݁ՌΛใࠂͨ͠ɽ

ࢀߟจ‫ݙ‬[1] Dahabreh IJ, Robertson SE, Tchetgen EJ, Stuart EA, Hernán MA. Generalizing causal inferences fromindividuals in randomized trials to all trial-eligible individuals. Biometrics. 2019a; 75: 685–694.

1[2] Dahabreh IJ, Robertson SE, Hernán MA. On the relation between g-formula and Inverse probability
weighting estimators for generalizing trial results. Epidemiology. 2019b; 30: 807-812.

[3] Dahabreh IJ, Robertson SE, Steingrimsson JA, Stuart EA, Hernán MA. Extending inferences from arandomized trial to a new target population. Statistics in Medicine. 2020; 39: 1999–2014.

[4] Lesko CR, Buchanan AL, Westreich D, Edwards JK, Hudgens MG, Cole SR. Generalizing study results:
a potential outcomes perspective. Epidemiology. 2017; 28: 553–561.

[5] Tabuchi T, Kiyohara K, Hoshino T, Bekki K, Inaba Y, Kunugita N. Awareness and use of electronic
cigarettes and heat-not-burn tobacco products in Japan. Addiction. 2016; 111: 706–713.

[6] Ferri-Garcı́a R, Rueda M, Cabrera-León A. Self-perceived health, life satisfaction and related factors
among healthcare professionals and the general population: analysis of an online survey, with propensity
score adjustment. Mathematics. 2021; 9: 791; Available at: https://doi.org/10.3390/math9070791.

2

トーリックモデルからの直接抽出の代数的アルゴリ
ズム

Z

間野

修平

(統計数理研)∗1高山信毅(神戸大理)∗2

状態空間を [m] := {1, 2, . . . , m} とし，行空間に (1, . . . , 1) を含む行列 B = (aij ) ∈d×mと h ∈ Rm
>0 が定める確率分布の集合

MB,h := {p ∈ int(∆m−1 ) : log p ∈ log h + rowspan(B)}をトーリックモデル，もしくは離散指数型分布族，対数アフィンモデル，という．
!

母数 ψ ∈ Rd>0 と φ ∈ Rd>0 があり，φ を局外母数とする．

B=

!

A
Ã

"

!

!

∈ Z(d+d )×m , A = (aij ) ∈ Zd×m , Ã = (ãij ) ∈ Zd ×m ,(1, . . . , 1) ∈ rowspan(A) に対し，MB,h を
pj =# a # ãhj
φi ij
ψk kj ,Z(φ, ψ)!
i∈[d]k∈[d ]j ∈ [m]とパラメトライズする．多項抽出，もしくはポアソン抽出された標本における状態 j をとる観察の数を uj j ∈ [m] で表すとき，A, b が定める集合

Fb (A) := {u : Au = b, u ∈ Nm },を φ に対する最小十分統計量 b ∈ NA :=$N := {0, 1, 2, . . .}j∈[m] Naj に付随する b ファイバーという．

ここで aj は A の j 列ベクトルである．b による条件付き確率分布はP(U = u|AU = b) =である．母数を yj := hj

1yu,
ZA (b; y) u!
%u ∈ Fb (A),y u :=#uyj j , u! :=

j∈[m]

#

uj !

(1)

j∈[m]

ãkj
∈ R>0 とした．正規化定数k∈[d! ] ψk

ZA (b; y) :=

& yu
u!

u∈Fb (A)は GKZ 超幾何多項式，もしくは A 超幾何多項式と呼ばれる．整数格子 NA の凸包は多面体錘を成し，b ∈/ NA のとき ZA (b; y) = 0 と規約する．(1) の確率函数が表す分布を A
超幾何分布と呼ぶ．
∗1〒190-8562 東京都立川市緑町 10-3∗2e-mail: smano@ism.ac.jp〒657-8501 兵庫県神戸市灘区六甲台町 1-1
e-mail: takayama@math.kobe-u.ac.jp

統計数理研究所 数理・推論研究系
神戸大学 大学院理学研究科

NA に埋め込まれた有界な整数格子で半順序v ∈ NA and v − aj ∈ NA⇒v − aj ≺ v

を備え，最大元 b と最小元 0 を持つものをマルコフ束 LA (b) と呼ぶ．以下のアルゴリズ
ムは状態空間を LA (b) とするマルコフ連鎖である．状態は元に 1 対 1 に対応し，推移は近傍，つまり u ≺ v を満たす元の組において v から u に起こる．

アルゴリズム 1 (Gauss–Manin ベクトル q(b; y), 行列 T (b; y) は講演で与えた).
入力: 行列 A ∈ Zd×m , 十分統計量 b ∈ N0 A, パラメタ y ∈ Rd>0 .
出力: A, b が定める A 超幾何分布に従うランダムベクトル u.

Step 1: Gauss–Manin ベクトル q(b; y) を初期化．
Step 2: t = 1, n = deg(b) とする．
Step 3: 規格化しない推移確率 ẽ(b, b − aj ; y) = {T (b; y)q(b; y)}j , ∀j ∈ [m] を計算．
Step 4: [0, 1] を比 ẽ(b, b − a1 ; y) : ẽ(b, b − a2 ; y) : · · · : ẽ(b, b − am ; y) に分割．
b − aj ∈/ LA (b) のとき j 番目は 0.

Step 5: [0, 1] の一様乱数をとり，区間 ẽ(b, b − aj ; y) に入れば jt = j とする．

Step 6: deg(b − ajt ) = 0 であれば Step 7, そうでなければ qk (b − ajt ; y) =yj−1
{P̃jt (b; y)q(b; y)}k , k ∈ 0 ∪ [r − 1] を計算，t ← t + 1, b ← b − ajt とし，Step 3.
t
Step 7: (u1 , . . . , um ), uj := |{t : jt = j, t ∈ [n]}| を出力．

行列 P̃j (b; y) はいかなる行列 A に対してもある微分作用素環のイデアルの Gröbner
基底に対する標準形を計算することで得られる．{0, a1 , . . . , am } の凸包を単位単体{0, e1 , . . . , ed } で測った体積を r = vol(A) で表した．

定理 1. アルゴリズム 1 による A 超幾何分布に従うランダムベクトルの直接抽出の時間
計算複雑性は O(max{m, r}rn) である．ここで n = deg(ZA (b; y)), r = vol(A). ただし
有理数演算のコストを O(1) とした．
本講演ではアルゴリズム 1 とその導出について説明し，二元分割表の独立性のないモデルについて計算機実験の結果を報告した．結果をまとめる．メトロポリス連鎖の有効サンプルサイズと比較して，分割表 1 個の生成にかかる時間は 800 倍ほどであった．計算に要した時間は，現れる推移確率の計算を尽くした後は，分割表の個数にほぼ比例した．モジュラ算法に基づく並列化の有効性を検証した．推移確率の計算はコア数を増やすだけでは短縮できないが，それを終了した後の計算時間はコア数にほぼ反比例した．

参考文献[1] Mano, S., Takayama, N., Algebraic algorithm for direct sampling from toric models,arXiv: 2110.14992