ISSN

0385-2148

研究所報
No.54

統計的モデリング

２０２１年 １１月

法政大学

日 本 統 計 研 究 所

ISSN

0385-2148

研究所報
No.54

統計的モデリング

２０２１年 １１月

法政大学

日 本 統 計 研 究 所

001

002

はじめに
第 1 回目は統計的モデリングについての研究集会を行った。統計的モデリング、と言っても関連する分野は幅広くあるので、本研究集会では特に角度統計やコピュラを用いたモデリングに関連する研究者らに最近の研究内容の紹介をしていただくことにした。
最初に、小方浩明先生（東京都立大学）に直線上の多次元高次マルコフデータに対するヴァイ
ンコピュラモデリングを紹介していただき、それの円周上のデータに対する方法についての提案やその課題について講演していただいた。次に、シリンダーモデリングの分野で世界的に活躍している気鋭の研究者である井本智明先生（静岡県立大学）にシリンダーモデルやトーラス上のモデリン
グについて最新の統計モデルの提案、サッカーなどのスポーツのデータに対するモデルの適用例を講演していただいた。塩濱敬之先生（南山大学）には高次の角度マルコフ過程の紹介や mixture
transition distribution modeling について講演していただいた。宮田庸一先生（高崎経済大学）には非凸な罰則項によるスパース推定、モデルの拡張型確率的コンプレキシティに対するラプラス近似に関する研究について講演していただいた。最後に江村剛志先生（久留米大学）にはコピュラとフレイルティ混合効果から導かれる多変量故障時間分布について講演していただいた。
本研究集会では講演中にも活発にディスカッションがされ、また、発表後も研究情報の意見交換ができ、大変有意義な研究集会となった。

2021 年 1１月 日本統計研究所

i
i

目次
講演

円周上の高次マルコフ過程におけるペアサーキュラモデリング
1
小方浩明
講演

シリンダー分布構成法の提案について
6
井本智明
講演

Mixture transition distribution modeling for higher order circular Markov processes
8
塩濱敬之
講演

非凸な罰則項によるスパース推定における拡張型確率的コンプレキシティについて
10
宮田庸一

講演

コピュラと混合効果から生成される多変量故障時間分布
19
江村剛志

ii
ii

ԁप্ͷߴ࣍Ϛϧίϑաఔʹ͓͚ΔϖΞαʔΩϡϥϞσϦϯά*౦‫ཱ౎ژ‬େֶɹখํߒ໌

͸͡Ίʹ

1

ํ޲౷‫ֶܭ‬͸ԁप্ͷ఺ʹΑͬͯද͞ΕΔप‫ظ‬తͳσʔλΛѻ͏౷‫ֶܭ‬ͷҰ෼໺Ͱ͋Δɽ෩޲ͳͲͷ֯౓Λද͢σʔλ΍ɼපӃʹӡͼࠐ·ΕΔ‫ऀױ‬ͷ౸ண࣌ࠁͳͲͷɼ24 ࣌ؒද‫Ͱه‬ද͞ΕΔ࣌ࠁͷσʔλͳͲ͕୅දతͳྫͰ͋Δɽํ޲౷‫ؔ͢ʹֶܭ‬ΔςΩετͱͯ͠͸ɼྫ͑͹ Jammalamadaka
and SenGupta (2001) ΍ Mardia and Jupp (1999) ౳͕‫͛ڍ‬ΒΕΔɽઌʹ‫ͨ͛ڍ‬෩޲σʔλͳͲ͸ɼ͋ΔҰఆ‫ूʹؒظ‬ΊΒΕͨ࣌‫ྻܥ‬σʔλͰ͋Δ͜ͱ͕ଟ͘ɼΏ͑ʹ‫ྻܥ‬૬ؔΛ࣋ͭͱߟ͑Δͷ͕ࣗવͰ͋ΔɽຊใࠂͰ͸ɼԁप্ͷଟ࣍‫ྻܥ࣌ݩ‬σʔλͷϞσϦϯάΛߟ࡯͢Δɽ·ͣ͸ಉ࣌෼෍Λɼ֤पล෼෍ͱɼίϐϡϥͷԁप্൛ʹ૬౰͢Δଟ࣍‫ݩ‬αʔΩϡϥεͷੵʹ෼ղͯ͠ද‫͢ݱ‬Δɽଟ࣍‫ݩ‬ͷαʔΩϡϥεΛߋʹϖΞαʔΩϡϥεͷੵʹ෼ղ͢Δ͜ͱʹΑΓɼ֤֬཰ม਺ͷैଐੑ͕ղऍ͠΍͍֨͢޷ͱͳΔɽ͜Ε͸௚ઢ্ͷग़࣌‫ྻܥ‬σʔλʹର͢ΔɼϰΝΠϯίϐϡϥʹΑΔϞσϦϯάʹྨࣅ͍ͯ͠Δɽ࣮ࡍʹ͸ɼ࣌‫ྻܥ‬σʔλʹ‫ڧ‬ఆৗੑͱߴ࣍ͷϚϧίϑੑΛ՝͢͜ͱʹΑΓɼϞσϧͷେ͖͞Λ੍‫͠ݶ‬ɼ‫࣮ݱ‬తʹਪఆՄೳͳϨϕϧ·Ͱམͱ͠ࠐΉɽຊใࠂͷߏ੒͸ҎԼͰ͋Δɽ2 ষͰɼSmith (2015) ʹΑΔ௚ઢ্ͷଟ࣍‫࣍ߴݩ‬Ϛϧίϑσʔλʹର͢ΔϰΝΠϯίϐϡϥϞσϦϯάΛ঺հ͠ɼ3 ষͰԁप্ͷσʔλʹର͢Δํ๏ΛఏҊ͢Δɽ
4 ষͰ؆୯ͳγϛϡϨʔγϣϯ݁ՌΛఏࣔ͠ɼ5 ষʹ·ͱΊͱ՝୊Λड़΂ΔɽϰΝΠϯίϐϡϥʹΑΔଟ࣍‫ݩ‬ϚϧνΦʔμʔϚϧίϑաఔ2ຊষͰ͸ɼSmith (2015) ʹΑΔɼଟ࣍‫ྻܥ࣌ݩ‬σʔλʹ͓͚ΔैଐੑΛϰΝΠϯίϐϡϥʹΑͬͯද‫͢ݱ‬ΔϞσϧΛ঺հ͢ΔɽҎԼͷΑ͏ͳ R ্ͷ m ࣍‫ྻܥ࣌ݩ‬σʔλ Yt = (Y1,t , . . . , Ym,t )′ ,Y = (Y1′ , . . . , YT′ )′ Λߟ͑ΔɽY ͸ N = T m ࣍‫ݩ‬ͷ֬཰ϕΫτϧͰ͋Γ, ఴ͑ࣈͷ؆ུԽͷͨΊҎ߱Ͱ͸ Y = (Y1 , . . . , YN )′ ͱ͍͏Α͏ͳද‫ݱ‬΋༻͍Δɽͭ·Γɼఴ͑ࣈ͕ΧϯϚͰ෼͚ΒΕ͍ͯΕ͹ཁૉͱ࣌ࠁͷΠϯσοΫεΛ໌ࣔͨ͠‫͋Ͱܗ‬ΓɼΧϯϚ͕ͳ͚Ε͹શͯͷཁૉΛҰྻʹฒ΂ͯ·ͱΊͯ൪߸Λৼ͍ͬͯΔͱཧղ͞Ε͍ͨɽ͜ͷͱ͖ɼY ͷಉ࣌ີ౓ؔ਺͸ҎԼͷΑ͏ʹ༩͑ΒΕΔɽ

f (y) = c(u)

T ∏
m
∏

f (yj,t ).

t=1 j=1͜͜ʹɼu′ = (u1 , . . . , uN ) = (F (x1 ), . . . , F (xN )), f (yj,t ) ͸ yj,t ͷपลີ౓ؔ਺ɼF (yi ) = F (yj,t )͸ରԠ͢Δ෼෍ؔ਺ɼͦͯ͠ c(u) ͸ U = (U1 , . . . , UN )′ = (F (Y1 ), . . . , F (YN ))′ ͷಉ࣌ີ౓ؔ਺
Ͱ͋ΓɼҰൠʹίϐϡϥີ౓ؔ਺ͱ‫ݺ‬͹ΕΔ΋ͷͰ͋Δɽίϐϡϥີ౓ؔ਺͸ [0, 1]N ্ͷؔ਺Ͱ
͋Δ͕ɼҎԼͷΑ͏ʹ N (N − 1)/2 ‫ݸ‬ͷϖΞίϐϡϥʢD-ϰΝΠϯίϐϡϥͱ‫ݺ‬͹ΕΔʣʹ෼ղ͢Δ͜ͱΛߟ͑Δɽ

c(u) =

N ∏i−1∏ci,j (ui|j+1 , uj|i−1 ).

i=2 j=1*ຊ‫ڀݚ‬͸ JSPS Պ‫ݚ‬අ 18K11193 ͷॿ੒Λड͚ͨ΋ͷͰ͢ɻ

1

1

(1)

͜͜ʹɼui|i ≡ ui Ͱ͋Γ, i > j ʹରͯ͠͸ ui|j+1 ≡ F (ui |ui−1 , . . . , uj+1 ), uj|i−1 ≡ F (uj |ui−1 , . . . , uj+1 )Ͱ͋Δɽc(u) ʹൺ΂ ci,j (ui|j+1 , uj|i−1 ) ͸ϖΞ͝ͱͰ͋ΔͷͰɼ૬ؔߏ଄͕ΑΓղऍ͠΍͘͢ͳΔɽ·ͨɼϖΞίϐϡϥΛϒϩοΫ͝ͱʹάϧʔϐϯά͢Δ͜ͱʹΑͬͯɼཁૉؒͰͷ૬ؔͱ‫ྻܥ‬૬ؔΛ෼͚ͯද‫͢ݱ‬Δ͜ͱ͕Ͱ͖Δɽut = (u1,t , . . . , um,t )′ ͱ͠ɼҟͳΔ࣌ࠁؒ t1 > t2 ʹରͯ͠b(t1 )

Kt1 ,t2 (ut2 , . . . , ut1 ) =∏b(t2 )∏

ci,j (ui|j+1 , uj|i−1 )i=a(t1 ) j=a(t2 )ͱఆٛ͠ɼಉ࣌ࠁ t ʹରͯ͠͸

Kt,t (ut ) =

b(t)
∏

i−1
∏

ci,j (ui|j+1 , uj|i−1 )i=a(t)+1 j=a(t)ͱఆٛ͢Δɽ͜͜ʹ a(t) = (t − 1)m + 1, b(t) = tm Ͱ͋Δɽ͢Δͱ
{T} { T t−1}∏∏∏c(u) = c(u1 , . . . , uT ) =Kt,t (ut )Kt,t−i (ut−i , . . . , ut )
t=1

(2)

t=2 i=1ͱද‫͖Ͱݱ‬ΔɽҰ൪໨ͷத‫͕಺ހׅ‬ಉ࣌ࠁʹ͓͚ΔཁૉؒͷैଐੑΛද‫͢ݱ‬Δ෦෼Ͱ͋Γɼೋ൪໨ͷத‫͕಺ހׅ‬ҟͳΔ࣌ࠁؒͰͷैଐੑɼ͢ͳΘͪ‫ྻܥ‬૬ؔΛද‫͢ݱ‬Δ෦෼Ͱ͋Δɽ
͜͜Ͱ࣌‫ڧʹྻܥ‬ఆৗաఔͷߏ଄ΛೖΕΔͱɼ(2) ࣜதͷؔ਺ Kt,t ฒͼʹ Kt,t−i ͸ t ʹґଘͤͣɼՃ͑ͯ p ࣍Ϛϧίϑੑʢա‫ڈ‬શͯͷσʔλΛॴ༩ͱͨ͠ͱ͖ͷ‫Ͱ఺࣌ݱ‬ͷ৚݅෇෼෍͕ɼա‫ ڈ‬p ࣌఺·ͰΛॴ༩ͱͨ͠ͱ͖ͷ‫Ͱ఺࣌ݱ‬ͷ৚݅෇෼෍ʹ౳͍͠ʣͷߏ଄ΛೖΕΔͱɼi > p ʹ͓͍ͯ Kk,k−i = 1 ͱͳΔͨΊɼ(2) ࣜ͸{T}  T min(t−1,p)∏ ∏∏c(u) = c(u1 , . . . , uT ) =K0 (ut )Ki (ut−i , . . . , ut )


t=1

t=2

i=1

ͱॻ͖‫͑׵‬ΒΕɼϞσϦϯά͢Δࡍʹ࢖༻͢ΔϖΞίϐϡϥؔ਺ͷ਺͕཈͑ΒΕΔɽ

3ԁप্ͷϚϧνΦʔμʔϚϧίϑաఔϰΝΠϯίϐϡϥͷ֓೦Λ༻͍ͯɼԁप্ͷϚϧνΦʔμʔϚϧίϑաఔΛ‫ن‬ఆ͢Δ͜ͱΛߟ͑Δɽԁप্ͷ෼෍ʹ͓͚Δίϐϡϥີ౓ؔ਺ʹ૬౰͢Δ΋ͷ͍ͭͯ͸ɼWehrly and Johnson (1980)
ͳͲ͞·͟·ͳ࿦จͰఏࣔ͞Ε͓ͯΓɼJones et al. (2015) ͕αʔΩϡϥε (circulas) ͱ͍͏‫Ͱܗ‬౷Ұతʹ·ͱΊ͍ͯΔɽԁप্ (Π Ͱද͢͜ͱͱ͢Δ) ͷ m ࣍‫ྻܥ࣌ݩ‬σʔλΛߟ͑, 2 ষʹ฿͍ Θt = (Θ1,t , . . . , Θm,t )′ ,
Θ = (Θ′1 , . . . , Θ′T )′ = (Θ1 , . . . , ΘN )′ ͷΑ͏ʹද‫͢ه‬Δɽͭ·Γ Θ ͸ ΠN ্Ͱͷ֬཰ϕΫτϧͱ
ͳΔɽΠ ্ͷ֤पลີ౓ؔ਺ɼपล෼෍ؔ਺ΛͦΕͧΕ fi (θi ), Fi (θi ) (i = 1, . . . , N ) Ͱද͢ɽ͜ͷͱ͖ɼΘ ͷಉ࣌ີ౓ؔ਺͸ҎԼͷΑ͏ʹॻ͚Δɽ

fΘ (θ) = (2π)N c(2πF1 (θ1 ), . . . , 2πFN (θN ))N∏fi (θi ).

(3)

i=1

͜͜ʹ c(θ1 , . . . , θN ) ͸ɼαʔΩϡϥεͱΑ͹ΕΔɼΠN ্Ͱఆٛ͞Εͨಉ࣌ີ౓ؔ਺Ͱ͋Γɼ࣍ͷೋͭͷੑ࣭Λ࣋ͭɽ22

(i) ֤पล෼෍͕ Π ্ͰͷҰ༷෼෍ɼ͢ͳΘ֤ͪपลີ౓ؔ਺͕ ci (θi ) = (2π)−1ͱͳΔɽ(i = 1, . . . , N )(ii) ֤ม਺ʹ͓͍ͯप‫ظ‬ੑc(θ1 + 2k1 π, . . . , θN + 2kN π) = c(θ1 , . . . , θN ), (θ1 , . . . , θN ) ∈ ΠN , k1 , . . . , kN = 0, ±1, ±2, . . .
Λຬͨ͢ɽ(ii) ͷप‫ظ‬ੑ͕ཁ੥͞Ε͍ͯΔ͜ͱʹΑΓɼαʔΩϡϥε͸ɼ௨ৗͷίϐϡϥΛ୯ʹϦεέʔϧͨ͠΋ͷͱ͸ҟͳΔɽ(1) ʹ฿͍ɼN ࣍‫ݩ‬αʔΩϡϥεΛϖΞαʔΩϡϥεʹ෼ղ͢Δͱc(θ1 , . . . , θN ) = (2π)−N

N ∏
i−1
∏

(2π)2 ci,j (2πθi|j+1 , 2πθj|i−1 )(4)

i=2 j=1

ͱͳΔɽ͜͜ʹ θi|i ≡ (2π)−1 θi Ͱ͋Γ, i > j ʹରͯ͠͸ θi|j+1 ≡ F (θi |θi−1 , . . . , θj+1 ), θj|i−1 ≡F (θj |θi−1 , . . . , uj+1 ) Ͱ͋Δɽ·ͨ ci,j (θi , θj ) ͸ Π2 ্ͷαʔΩϡϥεͰ͋Γɼ۩ମతʹci,j (θi , θj ) = (2π)−1 gi,j (θj − qi,j θi )(5)

ͱ͢Δ͜ͱʹΑΓαʔΩϡϥεͷ৚݅Λຬ଍͢Δɽ͜͜ʹ gi,j ͸ Π ্ͷ೚ҙͷີ౓ؔ਺Ͱ͋Γɼqi,j ∈ {−1, 1} ͸ɼθi ͷූ߸ΛܾΊΔඇ֬཰తͳཁૉͰ͋Δɽશͯͷ i, j ʹ͓͍ͯ gi,j (θ) = (2π)−1 (Π্ͷҰ༷෼෍) ͱઃఆ͢Δͱɼ(3)–(5) ΑΓfΘ (θ) =N∏fi (θi )i=1ͱͳΓɼಠཱϞσϧͱͳΔɽ·ͨ Jones et al. (2015) ͷ 2.3 ষʹ͓͍ͯɼgi,j ͷूத౓͕ߴ͚Ε͹Θi ͱ Θj ͷ૬͕ؔߴ͘ͳΔ͜ͱ͕‫͞ٴݴ‬Ε͓ͯΓɼgi,j ͷूத౓ύϥϝʔλΛͦͷ·· Θi ͱ Θj ͷ૬ؔ౓ͱͱΒ͑Δ͜ͱ͕Ͱ͖Δɽ·ͨɼϒϩοΩϯάʹ͍ͭͯ΋ 2 ষͱಉ༷ʹɼҟͳΔ࣌ࠁؒ t1 > t2 ʹରͯ͠b(t1 )

Kt1 ,t2 (θt2 , . . . , θt1 ) =∏b(t2 )∏(2π)2 ci,j (2πθi|j+1 , 2πθj|i−1 ),i=a(t1 ) j=a(t2 )·ͨಉ࣌ࠁ t ʹରͯ͠͸

Kt,t (θt ) =

b(t)
∏

i−1
∏

(2π)2 ci,j (2πθi|j+1 , 2πθj|i−1 )i=a(t)+1 j=a(t)ͱఆٛ͠ɼc(θ) = c(θ1 , . . . , θT ) = (2π)

−N

{T
∏

Kt,t (θt )t=1} { T t−1∏∏

Kt,t−i (θt−i , . . . , θt )t=2 i=1}ͱද‫͖Ͱݱ‬Δɽ·ͨɼ‫ڧ‬ఆৗੑͱ p ࣍ϚϧίϑੑΛԾఆͨ͠ͱ͖ʹɼ2 ষͱಉ༷ʹ}  T min(t−1,p)
{T
∏ ∏

∏
c(θ) = c(θ1 , . . . , θT ) = (2π)−NK0 (θt )Ki (θt−i , . . . , θt )


t=1

t=2

ͱද‫͖Ͱݱ‬Δɽ

3

3

i=1

4

γϛϡϨʔγϣϯ

γϛϡϨʔγϣϯʹΑͬͯɼm = 2 ͷ֯౓࣌‫ྻܥ‬σʔλΛ T = 50 ‫(( ݸ‬θ1,1 , θ2,1 )′ , . . . , (θ1,50 , θ2,50 )′ )ൃੜͤ͞Δɽཁૉؒʹ͓͍ͯ΋࣌఺ؒʹ͓͍ͯ΋ಠཱ͔ͭಉҰͷ෼෍ʹै͓ͬͯΓɼ֤पล෼෍
͸Ґஔύϥϝʔλ µ = 0, ूத౓ύϥϝʔλ κ = 3 ͷϑΥϯɾϛʔθε෼෍ʢvM(µ = 0, κ = 3) ͱද‫͢ه‬Δʣͱ͍͏ઃఆͰ͋ΔɽϑΥϯɾϛʔθε෼෍ͷ֬཰ີ౓ؔ਺͸ҎԼͰ༩͑ΒΕΔɽ

fvM (θ; µ, κ) =1
exp{κ cos(θ − µ)}.
2πI0 (κ)͜͜ʹɼIr (κ) ͸ୈ 1 छͷ r ࣍मਖ਼ϕοηϧؔ਺Ͱ͋Γɼ1
Ir (κ) =2π∫ 2πcos(rθ) exp(κ cos θ) dθ,0

r = 0, ±1, . . .

Ͱఆٛ͞ΕΔɽूத౓ύϥϝʔλ͕ κ = 0 ͷͱ͖ɼԁप্ͷҰ༷෼෍ʹҰக͢Δɽ্‫ه‬ͷΑ͏ʹͯ͠ൃੜͤ͞ΒΕͨγϛϡϨʔγϣϯσʔλʹରͯ͠ɼ3 ষͰߟ࡯ͨ͠ϞσϧΛ͋ͯ͸ΊΔɽ۩ମతʹ͸ɼ·ͣ p = 1 ࣍ͷϚϧίϑ‫ڧ‬ఆৗաఔΛ͋ͯ͸Ίɼ֤पล෼෍ʹ fi ∼ vM(µ = 0, κi )(i = 1, 2) Λ͋ͯ͸ΊΔɽ·ͨɼϖΞαʔΩϡϥεΛ (5) ࣜͰߏங͠ɼಉ࣌ࠁʹ͓͚Δ෦෼ʹ͓͍
(1)
(0)
(0)
(1)
ͯ͸ g2,1 ∼ vM(µ = 0, κ2,1 )ɼϥά 1 ͷ࣌ࠁؒʹ͓͚Δ෦෼ʹ͓͍ͯ͸ gi,j ∼ vM(µ = 0, κi,j )(i, j = 1, 2) Λ͋ͯ͸ΊΔɽ(5) ࣜ಺ͷූ߸ΛܾΊΔඇ֬཰తͳ෦෼ʹ͍ͭͯ͸ɼ͢΂ͯ q = 1 ͱ
(0)
(1)
(1)
(1)
(1)
ઃఆ͢Δɽσʔλੜ੒ͷ࢓ํʹΑΓɼਅ஋͸ κ1 = κ2 = 3, κ2,1 = κ1,1 = κ1,2 = κ2,1 = κ2,2 = 0 ͱͳΔɽ࣮ࡍͷύϥϝʔλਪఆ͸ Rstan ʹΑΔϚϧίϑ࿈࠯ϞϯςΧϧϩ๏ (MCMC) Ͱߦ͏ɽMCMC
ͷઃఆ͸ chain = 2, iter = 1500, warmup = 150, iter = 1 Ͱ͋Γɼશͯͷύϥϝʔλʹ͓͍ͯແ৘ใࣄલ෼෍ͱ͢Δɽ֤ύϥϝʔλʹ͓͚Δࣄ‫ޙ‬෼෍ͷཁ໿Λද 1 ʹࣔ͢ɽ͓͓ΉͶྑ޷ͳਪఆ
݁Ռ͕ಘΒΕͨɽද 1: MCMC ʹΑΔ֤ύϥϝʔλʔͷࣄ‫ޙ‬෼෍ཁ໿ฏ‫ ۉ‬ඪ४ภࠩ 2.5% ఺ 50% ఺ 97.5% ఺

κ1
κ2
(0)
κ21
(1)
κ11
(1)
κ12
(1)
κ21
(1)
κ22

5

3.28
3.25
0.28
0.16
0.16
0.22
0.15

0.58
0.58
0.18
0.12
0.12
0.15
0.12

2.26
2.24
0.02
0.01
0.01
0.01
0.01

3.26
3.21
0.25
0.14
0.13
0.19
0.13

4.48
4.46
0.70
0.44
0.46
0.57
0.45

·ͱΊͱ՝୊

ຊ‫Ͱڀݚ‬͸ɼԁप্ͷϚϧνΦʔμʔϚϧίϑաఔʹରͯ͠αʔΩϡϥεΛ༻͍ͨϞσϦϯά
ΛఏҊͨ͠ɽαʔΩϡϥε͸ԁप্൛ͷίϐϡϥʹରԠ͢Δ֓೦Ͱ͋Δɽଟ࣍‫ݩ‬ͷαʔΩϡϥε͸ϖΞαʔΩϡϥεʹ෼ղ͢Δ͜ͱʹΑΓɼ֬཰ม਺ͷैଐੑ͕ΑΓղऍ͠΍͍͢ϞσϦϯάͱͳΔɽϖΞαʔΩϡϥε͸೚ҙͷԁप্֬཰ີ౓ؔ਺ʹΑͬͯߏஙͰ͖ɼͦͷ֬཰ີ౓ͷूத౓͕ͦͷ··ैଐੑΛද͢ࢦඪͱͳΔɽ44·ͨɼγϛϡϨʔγϣϯʹΑͬͯൃੜͤ͞ΒΕͨಠཱಉҰ෼෍ʹै͏ 2 ࣍‫ݩ‬ͷ֯౓࣌‫ྻܥ‬σʔλʹରͯ͠ɼԁप্ͷ 1 ࣍Ϛϧίϑ‫ڧ‬ఆৗաఔΛ͋ͯ͸ΊɼMCMC ʹΑͬͯύϥϝʔλʔͷࣄ‫ޙ‬෼෍Λ‫ͨݟ‬ɽ͓͓ΉͶྑ޷ͳਪఆ݁Ռ͕ಘΒΕͨɽ‫Ͱ఺࣌ݱ‬͸͘͝؆୯ͳγϛϡϨʔγϣϯ͔࣮͠ߦͰ͖͓ͯΒͣɼࠓ‫ޙ‬͸࣮σʔλͳͲʹରͯ͠ΑΓߴ࣍ͷԁप্ϚϧίϑաఔΛ͋ͯ͸ΊɼύϥϝʔλʔΛਪఆ͢Δ༧ఆͰ͋ΔɽͦͷͨΊʹ͸‫ࢉܭ‬ෛՙͷܰ‫͕ݮ‬ඞਢͰ͋Δɽ۩ମతʹ͸ɼ(4) ࣜͷҾ਺Ͱ͋Δ৚݅෇෼෍ؔ਺ͷ஋Λ‫ٻ‬ΊΔ෦෼ʹ‫ࢉܭ‬ෛՙ͕͔͔ΔͱࢥΘΕɼͦͷ෦෼ͷ‫ࢉܭ‬ͷ޻෉͕ඞཁͱͳΔɽࢀߟจ‫ݙ‬
Jammalamadaka, S. A., SenGupta, A. S. (2001). Topics in circular statistics. World Scientific
Publishing Co., New York.
Jones, M.C., Pewsey, A., Kato, S. (2015). On a class of circulas: copulas for circular distributions.
Annals of the Institute of Statistical Mathematics, 67, 843–862.
Mardia, K. V., Jupp, P. E. (1999). Directional statistics. Wiley, Chichester.
Smith, M. S. (2015). Copula modelling of dependence in multivariate time series. International
Journal of Forecasting, 31, 815–833.
Wehrly, T. E., Johnson, R. A. (1980). Bivariate models for dependence of angular observations
and a related Markov process. Biometrika, 67, 255–256.

55

シリンダー分布構成法の提案について
井本智明
静岡県立大学 経営情報学部ある地点での風向や樹木の生長方向のような角度に関するデータは様々な社会科学分野で現れる。こうしたデータを適切に分析するためには，正規分布やポアソン分布のように直線上で定義される確率分布ではなく，角度の持つ周期性を考慮に入れた円周上の確率分布を用いたモデリン
グが必要となる。また，ある地点の温度と風向の組合せからなるデータや樹木の配置地点と生長方向の組合せからなるデータのように，直線上変量と円周上変量の組合せからなるデータを分析する際にはシリンダー上分布が必要となる。
本報告では，平均 0，分散 1 の実数値全体を台とする直線上対称確率密度関数 gL (·) と平均方向
0 の (−π, π] を台とする円周上確率密度関数 gC (·) から構成される確率密度関数
)}
)
(
]
[
{
(
x − µL1x − µLsin(θ − µC ) gL
gC (θ − µC ),f (x, θ) = 1 + λ sin 2Arctan
τ
τ
τ
[
]
(
)
)
(
L
2λ x−µsin(θ − µC ) 1
x − µLτ= 1+
gL
gC (θ − µC ),
(1)
(
)2
τ
τ
1 + x−µLτ− ∞ < x < ∞, −π < θ ≤ πを持つシリンダー上確率変数 (X, Θ) について考える。この関数に含まれる −1 ≤ λ ≤ 1 は直線上
変量 X と円周上変量 Θ の相関を操るパラメータ，−∞ < µL < ∞ と τ > 0 はそれぞれ X の平均と標準偏差を表すパラメータ，−π < µC ≤ π は Θ の平均方向を表すパラメータである。GL (·) とGC (·) をそれぞれ gL (·) と gC (·) の分布関数とすると，分布 (1) は接合関数をc(x, θ) = 1 +

−1
2λG−1
L (x) sin GC (θ)
2
1 + {G−1
L (x)}としたときのコピュラ分布としてみることもできる。
この分布の X と Θ の周辺確率密度関数 fX (x) と fΘ (θ) は)(x − µL1,fΘ (θ) = gC (θ − µC ),fX (x) = gLττであることがわかる。つまり，直線上変量が実数値全体で定義されており，その分布が対称であることを仮定するならば，提案法 (1) を用いることで周辺分布を任意に設定してシリンダー上分布の構成が行える。しかも，この方法によって構成された分布には複雑な関数が追加的に含まれることがないので統計推測のための計算コストを抑えることができる。また，Θ を所与としたと
きの X の条件付き確率密度関数 fX|Θ (x | θ) や X を所与としたときの Θ の条件付き確率密度関数
fΘ|X (θ | x) も容易に
]
[
)
(
)
(
L
sin(θ − µC ) 1
2λ x−µx − µLτg,
fX|Θ (x | θ) = 1 +
L
)2
(
L
τ
τ
1 + x−µτ][()
L
2λ x−µsin(θ − µC )
τfΘ|X (θ | x) = 1 +gC (θ − µC )()2L
1 + x−µτ66となることがわかる。特に後者は，[1] による sine-skewed 分布族に属する分布となっていることから回帰分析を行う際にも役立つ性質を有している。
√
確率変数 X と Θ についてのモーメントを，i = −1 を虚数単位として，][
X − µL= µp ,
E
E[eip(Θ−µC ) ] = αp + iβpτとおくと，分布 (1) の同時モーメントは] {[X − µL iq(Θ−µC )µp (αq + iβq ),p が偶数のとき,
e=EτλAp {(βq+1 − βq−1 ) − i(αq+1 − iαq−1 )}, p が奇数のとき,∫ ∞xp
f (x)dx である。これより，X と (cos Θ, sin Θ)
2 X−∞ 1 + xの正準相関係数として定義される X と Θ に関するシリンダー上相関係数 [2, 3] は，と表されることがわかる。ここで，Ap =2rXΘ
= 2λ2 A21 (1 − α2 )と表されることがわかる。このようなシンプルさから，モーメントによって定義される特徴量に対する各パラメータの役割が確認しやすくなる利点もある。
こうした分布 (1) のシンプルな特徴から，スポーツにおけるのパス地点とパス方向からなるデー
タやある地域での樹木配置と生長方向からなるデータのようなシリンダー上データの分析を容易に実行できる可能性を秘めている。今後は，提案手法の多変量拡張や変量間の相関を強める拡張を考えることで，シリンダー上データの統計分析に役立たせていく予定である。

参考文献[1] T. Abe and A. Pewsey. Sine-skewed circular distributions. Statistical Papers, 52:683–707,
2011.
[2] K. V. Mardia. Linear–circular correlation coeﬃcients and rhythmometry.
63:403–405, 1976.

Biometrika,[3] T. E. Wehrly and R. A. Johnson. Bivariate models for dependence of angular observations
and a related Markov process. Biometrika, 66:255–256, 1979.

7

7

Mixture Transition Distribution Modeling forHigher-Order Circular Markov Processes
Takayuki Shiohama
Department of Data Science, Nanzan University
Hiroaki Ogata
Faculty of Economics and Business Administration
Tokyo Metropolitan University

Abstract
The stationary higher-order Markov process for circular data is considered. We employ
the mixture transition distribution (MTD) model to incorporate higher-order dependency
in the circular time series. The underlying circular transition distribution is based on the
Wehrly and Johnson’s bivariate circular models. The structure of the circular autocorrelation
function is found to be similar to the autocorrelation function of the autoregressive process
on the line. The validity of the model is assessed by applying it to a series of real directionaldata.
keywords : circular statistics, Markov process, mixture transition distribution.

1

Introduction

Statistical analysis related to the data taking values on unit vectors is called directional statistics, where the direction is more important than the magnitude. Directional data occurs inmany areas, namely geostatistics, natural sciences, environmental sciences, biology, information
science, genetics, among others. The basic manifold of the high dimensional sphere is the unit
circle which is embedded in the Euclidean space R2 , and the statistical methods for the unit
circle are known as circular statistics.
While many data taking values on the unit circle or sphere are characterized with time series
structure, the statistical analysis for circular time series are not fully developed so far. We canobserve many time series data for the circle, such as the wind or wave directions, time records
of a certain event, and animal movement trajectories, which show typical periodic patterns. As
for the correlation coeﬃcient of two consecutive events in circular data, Fisher and Lee (1983)
proposed circular correlation coeﬃcients, and considered basic circular time series models as
linked autoregressive and moving average models, and wrapping process for real valued timeseries (Fisher and Lee (1994)). The comprehensive analysis of circular time series modeling canbe found in Breckling (1989).
One of the major approaches for modeling circular time series is to apply the Markov models
for the circular transition densities. Wehrly and Johnson (1980) proposed joint distribution forthe circular random variables. Holzmann et al. (2006) applied it for time series modeling using
Hidden Markov models. Abe et al. (2017) showed the circular autocorrelation structure of thecircular Markov process proposed by Wehrly and Johnson (1980). Some other Markov based
circular time series are considered by Kato (2010) where he considered the Möbius transformationof the circular random variables to produce the transition densities. His models are extended

88

by Jones et al. (2015), where they considered the circular joint distributions called ‘circulas‘ foranalyzing correlated structures for circular random variables.
Whilst all these approaches are based on the ﬁrst order Markov modeling for circular timeseries, investigating higher-order Markov models on the circle is primary importance. It is
worthwhile to mention that, as far as the authors knowledge, there exist no approaches for higherorder circular Markov models in literature. Modeling discrete valued Higher-order Markov chainsare very useful for the analysis of complex temporal relationships. Raftery (1985) considered the
mixture transition distribution (MTD) model. Since then, the MTD model has been developedand improved in various ways, see for details, Raftery and Tavaré (1994), Le et al. (1996), andBerchtold et al. (2002).
In this study, we develop the higher-order Markov process on the circle using the transition
density of Wehrly and Johnson (1980)’s model together with the MTD approaches. The autocorrelation structure of the proposed models are investigated. Since the probabilistic structureof the proposed models are quite complex and is not included in the class of regular parametricmodels, special attentions must be paid for parameter estimation. For this, we proposed themethods of maximum likelihood estimation, and provide some information criteria for model
selections. Since the MTD model is analogous to the standard autoregressive model in thatthe autocorrelations satisfy a system of linear equations similar to the Yule-Walker equations,the proposed higher-order circular Markov models have similar autocorrelation structures. The
asymptotic properties for the proposed models are considered.

References
Abe, T., H. Ogata, T. Shiohama, and H. Taniai (2017). Circular autocorrelation of stationary circularmarkov processes. Statistical Inference for Stochastic Processes 20 (3), 275–290.
Berchtold, A., A. Raftery, et al. (2002). The mixture transition distribution model for high-order markov
chains and non-gaussian time series. Statistical Science 17 (3), 328–356.
Breckling, J. (1989). The analysis of directional time series: applications to wind speed and direction,
Volume 61 of Lecture Notes in Statistics. Springer Science & Business Media.
Fisher, N. and A. Lee (1994). Time series analysis of circular data. Journal of the Royal StatisticalSociety: Series B (Methodological) 56 (2), 327–339.
Fisher, N. I. and A. Lee (1983). A correlation coeﬃcient for circular data. Biometrika 70 (2), 327–332.
Holzmann, H., A. Munk, M. Suster, and W. Zucchini (2006). Hidden markov models for circular andlinear-circular time series. Environmental and Ecological Statistics 13 (3), 325–347.
Jones, M., A. Pewsey, and S. Kato (2015). On a class of circulas: copulas for circular distributions.
Annals of the Institute of Statistical Mathematics 67 (5), 843–862.
Kato, S. (2010). A markov process for circular data. Journal of the Royal Statistical Society: Series B(Statistical Methodology) 72 (5), 655–672.
Le, N. D., R. D. Martin, and A. E. Raftery (1996). Modeling ﬂat stretches, bursts outliers in time series
using mixture transition distribution models. Journal of the American Statistical Association 91 (436),
1504–1515.
Raftery, A. and S. Tavaré (1994). Estimation and modelling repeated patterns in high order markov
chains with the mixture transition distribution model. Journal of the Royal Statistical Society: Series
C (Applied Statistics) 43 (1), 179–199.
Raftery, A. E. (1985). A model for high-order markov chains. Journal of the Royal Statistical Society:
Series B (Methodological) 47 (3), 528–539.
Wehrly, T. E. and R. A. Johnson (1980). Bivariate models for dependence of angular observations and a
related markov process. Biometrika 67 (1), 255–256.

99

非凸な罰則項によるスパース推定における拡張型確率的コンプレキシティについて
宮田 庸一

第一回 日本統計研究所

研究集会

統計的モデリングと統計推測理論

2021 年 8 月 12 日
概要

非凸な罰則関数を用いたスパース推定の一つであるブリッジ推定におけるチューニングパラ

メーターの選択に関して考える. Wang et al. (2009) によりベイズ型情報量規準によるチュー

ニングパラメーターの選択方法が提案されているが, モデルの拡張型確率的コンプレキシティ
に対するラプラス近似の主要項として正当化できることについて報告を行った.

1 はじめに記号 ′ は, 行列の転置を表すものとする. 観測ベクトル (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 + β j 1 x j 1 i + · · · + β j k x j k i + ϵ i ,高崎経済大学・経済学部, 〒370-0801 高崎市上並榎町 1300

1
10

10

(i = 1, ...n),ただし, j1 , ..., jk ∈ {1, ..., pn }, k = 1, ..., pn とし, ϵi (i = 1, ..., n) は独立に正規分布 N (0, σ 2 ) に従うとする. このとき, どのようにして適切な説明変数を選ぶのかが問題となるが, pn の次元がそれほど大きくない場合には, AIC(Akaike, 1973), BIC(Schwarz, 1978) などの情報量規準を用いて総当たり法を行うのが標準的なアプローチでとなる. しかしその一方で, pn の次元が高くなるにつれ,評価すべきモデルの数が指数的に増加するため, 総当たり法を行うのが困難になる. このような問題に対する一つの標準的なアプローチは, Tibshirani (1996) の LASSO(Least absolute shrinkage. 具体的には,and selection operator) に代表される罰則付き最小二乗推定量を用いることである

y = (Y1 , ..., Yn )′ , X1 = (xji ):n × pn 行列, 1n = (1, ..., 1)′ , X = 1n

X1 とおき, 損失関数

L0,n (β) = ∥y − Xβ∥2 + λn ∥β∥1

(1)を最小にする β = (β1 , ..., βpn )′ の値 β̂L を求めればよい. ただし, ∥β∥1 =j=1 |βj | とする.

 pn

こ

の推定量は LASSO 推定量と呼ばれ, λn → ∞ (n → ∞) とその他の適切な条件のもとで漸近正規

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

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

pn

j=1

|βj |γ ,(2)

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

なわち Huang et al. (2008) の結果は λn の選択に関しては, 一意に定まらないことになる. 一方で, Huang et al. (2008) で与えられた条件の下では, ブリッジ推定量の漸近的な良さが保証されていることもわかる. このため, 本報告では, Huang et al. (2008) の条件の下で λn を選ぶための規

準を, 拡張型確率的コンプレキシティの近似として与えられ, それが Wang et al. (2009) の規準と同様の形になることを示す. なお, 今回の話は Miyata (2021) の紹介となる.

2 ベイジアンアプローチ

ここでは, 前章の (2) 式で与えた, ブリッジ推定量をベイズ統計学の観点からの解釈するとともに，チューニングパラメーター λn のベイズ型の選択方法を提案する.

2
11

11

まずパラメーター β1 , ..., βpn は, 以下の指数ベキ (exponential-power) 分布 EP (0, λ0,n , σ 2 ) の確率密度関数に独立に従うとする (Box and Tiao, 1973, p.157).



1/γ
γλ0,n
λ0,n
γ
π(βj |λ0,n ) =exp − 2 |βj | .
2σ2Γ(1/γ)(2σ 2 )1/γ

λ0,n は定数でもよいが, サンプルサイズ n に依存してもよい形にする. このため, β の事前密度関数はπ(β|λ0,n ) = π(β0 )γ(λ0,n )1/γ
λ0,n
γ
exp
−
|β
|
.
j
2σ 2
2Γ(1/γ)(2σ 2 )1/γ
j=1
pn


で与えられる. π(β0 ) は β0 に対する事前分布で, 十分に滑らかなものとする. また記号を簡略化するため, 切片に対する事前分布を除いたものを π1 (β|λ0,n ) =および β を与えたときの Y

j=1 π(βj |λ0,n ) とする. そして, Xの条件付き確率密度関数は多変量正規分布 Nn (Xβ, σ 2 In ),p(y|X, β) =

p n



1
1
2
||Y
−
Xβ||
exp
−
.
2σ 2
(2πσ 2 )n/2に従うものとする. 通常, β の事後密度関数は p(y|X, β)π(β|λ0,n ) に比例する形になるが, ここではこの事後密度関数をさらに一般化して, β の事後密度関数は p(y|X, β)ϵn π(β|λ0,n ) に比例するものを考える. この疑似事後密度関数を pQ (β|X, y) とする. ただし ϵn は正の定数, もしく

は n → ∞ としたときに減少する列とする. なお ϵn が減少するということは, β に関する尤度

p(y|X, β) の情報を弱めることを意味している.

ここで, β の損失関数を Hn (β) = − log {p(y|X, β)ϵn π1 (β|λ0,n )} とおくと, この疑似事後密度関数は pQ (β|X, y) ∝ p(y|X, β)ϵn π(β|λ0,n ) = exp {−Hn (β)} π(β0 ) となる. これより, 切片の事

前分布 π(β0 ) を無視して exp {−Hn (β)} を最大化する疑似的な MAP(maximum a posteriori) 推

定量 β̂λn を考えると,

β̂λn = argmax {−Hn (β)}β= argmax {ϵn log p(y|X, β) + log π1 (β|λ0,n )}
β

pn



= argmin ||y − Xβ||2 + λn|βj |γ .
β

j=1となる. ただし λn = λ0,n /ϵn と置いた. つまり β̂λn は, 非凸な罰則項 λnリッジ推定量は同等なことがわかる.

pn

j=1 |βj |γを持つブ

次に, ブリッジ推定量により選ばれたモデルの良さを評価する. Ŝ := {j ∈ {1, ..., pn }|β̂λn ,j ̸= 0}は，β̂λn = (β̂λn ,0 , β̂λn ,1 , · · · , β̂λn ,pn )′ において非ゼロな推定量の添え字の集合とする. |Ŝ| は集合Ŝ における要素の個数を表す記号とする.

3
12

12

X(Ŝ) は, 集合 Ŝ の添え字に対応する説明変数の行列 X の部分行列とする. 例えば, pn = 3,β̂λ = (β̂λ,0 , β̂λ,1 , 0, β̂λ,3 )′ とすると,


x11
..
.

x21
..
.

1

x1n

x2n

1
 ..
X = .

x31..  ,
. 



1
 ..
X(Ŝ) =  .

x3n′β̂λ (Ŝ) = (β̂λ,0 , β̂λ,1 , β̂λ,3 ) ,

1

x11
..
.
x1n

β(Ŝ) = (β0 , β1 , β3 )

′


x31
.. 
. 

x3n



x11
 ..
X1 (Ŝ) =  .

x1n


x31
.. 
. 

x3n

となる.

3 拡張型確率的コンプレキシティ

ここでは, ブリッジ推定により選ばれたモデルの良さを評価するために, L(Yi , xi (Ŝ)′ β(Ŝ)) :=− log p(Yi |xi (Ŝ), β(Ŝ)) = (1/2) log(2πσ 2 )+{Yi −xi (Ŝ)′ β(Ŝ)}2 /(2σ 2 ) とおき, Yamanishi (1998)

により提案された拡張型確率的コンプレキシティ (Extended Stochastic Complexity, ESC) を考える:

1
ESC(Ŝ, λ0,n , ϵn ) := −
log
ϵn





exp −ϵnn
i=1

L(Yi , x′i (Ŝ)β(Ŝ))

π(β(Ŝ)|λ0,n )dβ(Ŝ).

(3)

ただし π(β(Ŝ)|λ0,n ) は, 事前密度関数 π(β|λ0,n ) における β(Ŝ) の周辺確率密度関数とする.

ESC を用いるときには, 関数 L(Yi , xi (Ŝ)′ β(Ŝ)) が, 観測される確率変数 Yi と, その予測値

Ŷi (Ŝ) := xi (Ŝ)′ β(Ŝ) の 損 失 関 数 の 形 に な っ て い る 必 要 が あ る が, 我 々 は 正 規 線 形 モ デ ル

を考えているためこの条件は見たされる. ここで (3) 式に対する解釈を述べる. 積分内のnexp −ϵn i=1 L(Yi , Ŷi (Ŝ)) であるが, これは大きな値であるほど, モデル Ŝ により与えられた予測量 Ŷi (Ŝ) で, Yi をうまく予測できることを表している. すなわち, (3) 式の積分は, 予測の良さ

を表す指標 exp −ϵn


L(Y
,
Ŷ
(
Ŝ))
を, 事前分布の下で平均を取ったものと解釈することが
i
i
i=1

n

できる. ここでは (3) 式に対する直接的な解釈を与えたが, 逐次型確率予測問題からの解釈が山西健司 (2010) の補題 16(p.114) で与えらえている. こちらの方が, ESC が事後平均の関数で与えられるという意味で, ベイズ的な解釈ができる.
なお (3) 式は

−

1
log
ϵn



R|Ŝ|+1

p(y|X(Ŝ), β(Ŝ))ϵn π(β(Ŝ)|λ0,n )dβ(Ŝ)(4)と表すことができるため, 対数周辺尤度を拡張した形と関連があることもわかる. もし ϵn = 1 のときには,

ESC(Ŝ, λ0,n , 1) := − log

R|Ŝ|+1

p(y|X(Ŝ), β(Ŝ))π(β(Ŝ)|λ0,n )dβ(Ŝ),

4
13

13

(5)

となるため, Ŝ に対応するモデルの周辺尤度に対数を取り, (−1) を乗じたものになることがわかる* 1 .

4 拡張型確率的コンプレキシティに対するラプラス近似

記号を簡略化するために, Σn = n−1 X1′ X1 , Σ1n = n−1 X1 (S0 )′ X1 (S0 ) とする. また ρ1n =λmin (Σn ), ρ2n = λmax (Σn ) とし, τ1n = λmin (Σ1n ), τ2n = λmax (Σ1n ) とおく. ただし,λmin (A)は行列 A の最小固有値を表す記号とし, λmax (A) は A の最大固有値を表す記号とする. このとき,以下の条件を仮定する:

1 nx2 = 1 (j = 1, ..., pn )n i=1 ji
(A2) ϵi (i = 1, 2, ...) は独立で同一分布に従い, E(ϵi ) = 0, V ar(ϵi ) = σ 2 , (0 < σ 2 < ∞) とする.
pn + λn kn→ 0 (n → ∞).
(A3) (a) ρ1n > 0 for all n, (b)nρ1n1/2λn k nλn ρ2−γ1n(A4) (a)→ 0 (n → ∞) (b)→ ∞ (n → ∞)1/2(2−γ)/2
nnγ/2 pn(A5) ある定数 0 < b0 < b1 < ∞ が存在して,
(A1)

n

i=1 xji = 0,b0 ≤ min{|β0,j | |1 ≤ j ≤ kn } ≤ max{|β0,j | |1 ≤ j ≤ kn } ≤ b1 .

(A6) X1 (S0 ) = (w1 , ..., wn ) とおく. このとき, ある正の定数 0 < τ1 < τ2 < ∞ が存在して,τ1 ≤ τ1n ≤ τ2n ≤ τ2 for ∀ n.

(b)

n−1/2 max wi′ wi → 0,1≤i≤nが成り立つ.

(n → ∞)(B1) 以下のことが成り立つ.
kn log n
→ 0,
ϵn n

and

ϵn n1/2 → ∞(n → ∞)(6)

条件 (A1)–(A6) は Huang et al. (2008) とまったく同じものである. 条件 (A2) は，線形回帰モデ

ルでは標準的な仮定である. 条件 (A3) (a) は ρ1n → 0 (n → ∞) であることは許していることから, フルモデルにおいては, Σn が漸近的に退化する状況, 例えば漸近的に多重共線性が起こるよ

うな状況は許していることになる. (A3) (b) は, ブリッジ推定量の一致性を示すときに必要になる. 条件 (A4) は, ブリッジ推定量がオラクル性，および漸近正規性を示すために必要になる. 条件(A6)(a) は, Σ1n が，n に関して一様に正値定符号であることを仮定している. また (A5)(b) は非ゼロの回帰係数に対するブリッジ推定量の漸近正規性を示す際に必要になる. 条件 (B1) は, 一般化

周辺尤度において漸近展開するために必要になるが, 真のモデルにおける説明変数の数 kn の増加のスピードはかなり遅いものと考えるのが一般的であるため，これはかなり弱い条件となる.
*1 小さな値を取るほど良いモデル.

5
14

14

ここで, 記号を簡略化するため, 標準化された疑似 log posterior を

h(β(Ŝ)) := −1
log p(y|X(Ŝ), β(Ŝ))ϵn π(β(Ŝ)|λ0,n )nとする. またその β(Ŝ) に関する 1 次の導関数を β(Ŝ) = β̂λn (Ŝ) で評価したものをDh : =
h(β(Ŝ))
∂β(Ŝ)∂(7)β(Ŝ)=β̂λn (Ŝ)=−ϵnλ0,n γ +ψ ,X(Ŝ)′ (y − X(Ŝ)β(Ŝ)) +σ2
2σ 2 n 1n−1で定義する. ただし, Σ̂+
X(Ŝ)′ X(Ŝ), k̂n := |Ŝ|,
1n = n+ψ1n:=0, |β̂λn ,1 |γ−1 sgn(β̂λn ,1 ), ..., |β̂λn ,k̂n |γ−1 sgn(β̂λn ,k̂n )′と お い た.

本 来, Ŝ の 集 合 は

ĵ1 , ĵ2 , ..., ĵk̂n と書くべきであるが, 記号を簡略化するために, Ŝ = {1, 2, ..., k̂n } とした. 実際, このような簡略化は一般性を失わずにできる. また



0





+
Ψ2n := γ(γ − 1) |β̂λn ,1 |γ−2

..



0
.

0

..

.
|β̂λn ,k̂n |γ−2とし, その β(Ŝ) に関する 2 次の導関数を β(Ŝ) = β̂λn (Ŝ) で評価したものを

D2 h :=

∂2










(8)



h(β(Ŝ))

∂β(Ŝ)∂β(Ŝ)′λ0,n +ϵnΨ+= 2 Σ̂+σ 1n 2σ 2 n 2nβ(Ŝ)=β̂λn (Ŝ)

とする. このとき, 以下が成り立つ.

定理 1 条件 (A1)–(A6), (B1) の下で, Ŝ, λ0,n , ϵn ) {1 + op (1)} ,exp −ϵn ESC(Ŝ, λ0,n , ϵn ) = exp −ϵn ESC(+が成り立つ. ただし |M | は行列 M の行列式を表すものとし, k̂n= k̂n + 1, Ŝ, λ0,n , ϵn ) = (2π)k̂n+ /2 nD2 h−1/2 p(y|X(Ŝ), β̂λ (Ŝ))ϵn π(β̂λ (Ŝ)|λ0,n )exp −ϵn ESC(
n
n

とする.

6
15

15

(9)

5 ベイズ型情報量規準

これまでは, σ 2 を既知として扱ってきたが, 実際にデータ解析を行うためには何らかの方法で推定する必要がある. ここでは, Wang et al. (2009) と同様にして

σ̂λ2 n =1||y − X β̂λ0,n /ϵn ||2 .
nを用いて推定を行う. ここで定理 1 で得られた近似式 (9) に対して, 対数をとり, (−2)/ϵn をかけるとk̂n + 1λ0,n 
log(nϵn ) + Op
2ESC(Ŝ, λ0,n , ϵn ) = n log σ̂λ2 n +|β̂λn ,q |γ +ϵn σ̂λ2 n
ϵn
q∈Ŝ



k̂n
ϵn



+ op (1).

となる. λ0,n と ϵn に対しては様々な取り方が考えられるが, ここでは λ0,n = λ0 (定数), ϵn ↓ 0(n → ∞), |Ŝ| = op (log n) とおき, Opk̂n log nϵn

ベイズ型情報量規準を得ることができる.



より低いオーダーの項を取り除くと, 以下の形のBtIC(Ŝ) = n log(σ̂λ2 n ) + k̂n+ log(nϵn )/ϵn= n log(σ̂λ2 n ) + Cn∗ k̂n+ log(n),(10)

1 + (log ϵn )/ log n
ϵn

(11)

ただし

Cn∗ =

となる. Cn∗ は, 条件 (A1)–(A6) の下では多項式オーダーで増加する項, 即ち Cn∗ = O(nαn ) (0 <αn < 1/2) となる. これに関しては, 次の章で述べる. 通常の BIC であれば, n log(σ̂λ2 n )+ k̂n+ log(n)となるので, 式 (10) は異なる形になっているが, これは Wang et al. (2009) で与えらえた情報量規準と同じ形になっている. Wang et al. (2009) の結果は, n log(σ̂λ2 n ) + Cn∗ |Ŝ| log(n) の形の情報量規準を形式的に作っておいて, それがモデル選択に関して一致性を持つような Cn∗ に対する条件を求めたものであるが, 今回, Huang et al. (2008) で与えられた条件の下で, 拡張型確率的コンプレキシティの近似として, Wang et al. (2009) の情報量規準を導出できたことにより, Wang et al.

(2009) の結果に対するある種の妥当性を与えたといえる.

6 Wang et al. (2009) の情報量規準との相違点

拡張型確率的コンプレキシティの近似から情報量規準 (10) が導出されたわけであるが, 厳密に考えた場合, Wang et al. (2009) の情報量規準との相違点がある. Wang et al. (2009) では, Cn∗ は多項式オーダーよりも低い任意のオーダーで増加することも認めている. またその論文のシミュレー
ションにおいては, Cn∗ = log(log pn ) が採用されている. 一方で, 我々の導出した結果においては,Cn∗ がどのようなオーダーで増加するかを考えてみる.
7
16

16

最初に ϵn = ϵ0 n−αn とおき, αn の取りうる範囲を求めよう. ただし, 話を簡単にするため,

kn = O(1), ρ1n = O(1) とする. このとき, 条件 (A4)(a), (b) は
λn
→ 0,
n1/2

λn
(2−γ)/2
nγ/2 pnとなる. これより,となる. もし γ = 1/2 とおくとき,

2−γ
1
γ
+
logn pn < αn <
2
2
2

1
1 3
+ logn pn < αn <
4 4
2

(12)

となる. この結果と (11) 式の形から, 我々の仮定の下では, Cn∗ は多項式オーダーで増加することがわかる. これは Wang et al. (2009) で提案されたものより速いオーダーで無限大にいくため, 実

際には定数 ϵ0 > 0 を調整する必要がある.

また Wang et al. (2009) においては, 線形モデルの説明変数は切片を持たない確率変数を考えているのに対して, Miyata (2021) においては, 説明変数は切片を持ち, 固定されたもの (すなわち計画行列) を考えている. この点においても相違点がある.

参考文献
Akaike, H. (1973). Information theory and an extension of the maximum likelihood principle.
In Second International Symposium on Information Theory (Tsahkadsor, 1971), pp. 267–
281.
Box, G. E. P. and G. C. Tiao (1973). Bayesian inference in statistical analysis. AddisonWesley Series in Behavioral Science: Quantitative Methods. Addison-Wesley PublishingCo., Reading, Mass.-London-Don Mills, Ont.
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.
Knight, K. and W. Fu (2000). Asymptotics for lasso-type estimators. Ann. Statist. 28 (5),1356–1378.
Miyata, Y. (2021). A Bayesian approach for selecting a tuning parameter in sparse linear
models (for submission).
Schwarz, G. (1978). Estimating the dimension of a model. Ann. Statist. 6 (2), 461–464.
Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. J. Roy. Statist. Soc.
Ser. B 58 (1), 267–288.
Wang, H., B. Li, and C. Leng (2009). Shrinkage tuning parameter selection with a diverging
number of parameters. J. R. Stat. Soc. Ser. B Stat. Methodol. 71 (3), 671–683.
8
17

17

Yamanishi, K. (1998). A decision-theoretic extension of stochastic complexity and its applications to learning. IEEE Trans. Inform. Theory 44 (4), 1424–1439.
Zou, H. (2006). The adaptive lasso and its oracle properties. J. Amer. Statist. Assoc. 101 (476),
1418–1429.

山西健司 (2010). 情報論的学習理論. 共立出版.

9
18

18

コピュラとフレイルティ混合効果から導かれる多変量故障時間分布
王尹辰（Yin-Chen Wang）
国立中央大学・統計研究所; peacebubu@gmail.com
江村剛志（Takeshi Emura）
久留米大学・バイオ統計センター; takeshiemura@gmail.com
要旨: コピュラモデルとフレイルティ（混合効果）モデルは, 多変量故障時間分布における相関構
造をモデリングするための一般的なモデル族である. 本稿は, コピュラモデルとフレイルティモデ
ルを含むさらに一般的な多変量故障時間モデル族を調査することを目的とする. この目的のために, 先行研究である Marshall and Olkin(1988)の多変量故障時間モデルの族を再訪する. 彼らの研究は, 特殊なコピュラのみ, また二変量モデルに限定されていた. これに対して, 我々の研究は, 応用上よく用いられるコピュラに焦点をあて, また３次元以上の多変量モデルも考察する. 部
分族として, 共有フレイルティモデルにより生成されるモデル族「フレイルティ-コピュラ」を定義・提案する. また, この族に属する個々の分布の名称を系統的に与える方法を提案する. また, 先行研究では考慮されてこなかった複雑なフレイルティ分布(対数正規, 切断正規, および折り畳み正規)の有用性も示す. 本稿は同著者による Wang and Emura(2021)を邦訳し, やや簡潔にまとめたものである.
キーワード: 二変量分布, コピュラ, FGM コピュラ, フレイルティ, 信頼性, 生存時間解析

1. はじめに寿命試験や臨床試験などでは, 観測対象から複数の「故障」や「死亡」に関する値が得られる. 観
測値間の相関構造をモデリングには, 多変量分布が必要となる(Duchateau and Janssen 2008;
Balakrishnan and Lai 2009; Crowder 2012; Emura et al. 2019; 江村‧道前 2020; 塚原 2021). 「コピュラ」と「フレイルティ」は, 多変量故障時間モデルの構築のための良く知られたツールである.
フレイルティモデルでは, 複数の故障時間に影響を与える未観測要因を考慮することで相関構造をモデリングする(Vaupel et al. 1979; Aalen 1994; Hougaard 1995). フレイルティモデルは, 様々なタイプの多変量故障時間の相関構造をモデリングするための一般的なツールとして重要な役割を担ってきた(Hougaard 1995; Duchateau et al. 2002; Duchateau and Janssen 2008; Rondeau etal. 2015; Ha et al. 2017; Ha and Lee 2021). フレイルティモデルでは, 未観測因子の分散は故障リスクの「不均一性」の程度を表すとともに, 故障時間の間での相関構造の強さも表す.
コピュラは, 不均一性の概念なしに純粋な相関構造を特徴づけるモデルを与える. コピュラ関
19

19

数(Nelsen 2006)は多数の周辺分布を結合し, 同時分布を生成する関数のことを指す. コピュラ関数の選択は, フレイルティモデルとは異なり, 周辺分布には影響しない. この性質により, 多変量
またはクラスター化された故障時間をモデリングするための明解な技法となっている (Shih 2014;
Prenen et al. 2018; Emura et al. 2019; Campos et al. 2021; Sofeu et al. 2021; kwon et al. 2021;
塚原 2021). スクラーの定理によって (Sklar 1959; Nelsen 2006), 任意の多変量故障時間分布はコピュラとして表現されるが, 一部の分布はフレイルティモデルとして表現されない場合がある.
コ ピ ュ ラ モ デ ル と フ レ イ ル テ ィ モ デ ル は 仕 組 み が 概 念 的 に 異 な る [see Section 3.3.4 ofDuchateau and Janssen (2008)]. 例えば, フレイルティモデルで定義された区分的定数ハザードモ
デル(e.g. Lo et al. 2017; Schneider al. 2020)は, コピュラモデルで定義された区分的定数ハザード
モデルと異なる(e.g. Emura and Michimae 2017; Lipowski et al. 2021)ため, 2 つのモデル異なる解
釈を持つ.
上述の議論は, 相関構造と不均一性という 2 つの異なる側面を統一する多変量故障時間モデルを構築する動機付けとなる. このようなモデルの研究は, コピュラとフレイルティが別々に議論されてきたため, 極めて限られている. Marshall and Olkin(1988)がコピュラとフレイルティの性質を同時に持つモデルを提案したが, 彼らの研究は特殊な二変量コピュラ(独立コピュラ, フレシェ-ヘフ
ディング限界コピュラ, FGM コピュラ)に焦点を当てた. 一方本稿では, クレイトンとグンベルコピュラなどの応用上重要なものや, 多変量 FGM コピュラに焦点を当てる. 部分族として, 共有フレイル
ティモデルにより生成されるモデル族「フレイルティ-コピュラ」を定義・提案する. また, この族に属する個々の分布の名称を系統的に与える方法を提案する. また, 先行研究では考慮されてこなかった複雑なフレイルティ分布(対数正規, 切断正規, および折り畳み正規)の有用性も示す.
本稿は以下のように構成されている. ２節で我々が提案するモデルを与える. ３節～５節では,クレイトンコピュラ, グンベルコピュラ, および多変量 FGM コピュラによってそれぞれ生成されたモデルを紹介する. ６節で結論と今後の課題を述べる.

2. 提案するモデル
本節では, 不均一性と相関構造という２つの異なる側面を持つ多変量故障時間モデルの基本的な枠組みを提案する.

2.1 節では, フレイルティでモデリングされる不均一性の概念を紹介する.

2.2 節では, フレイルティ項を所与とした条件付きコピュラモデルを紹介し, 「フレイルティ-コピュラモデル」という我々の提案する多変量故障時間分布族を定義する.
2.1 フレイルティ分布
故障のリスクに影響を与える未観測な因子を記述する「フレイルティ」の概念(Vaupel et al. 1979)を復習しよう. 部品の母集団を考えたとき, ある部品は壊れやすい一方で, 他の部品は頑健であろう. この不均一性をモデリングするために, リスクに影響を与える正の確率変数𝑍𝑍を考え, これをフ
20

20

レイルティ項と呼ぶ(Hougaard 1995; Duchateau and Janssen 2008; Ha et al. 2017). フレイルティ項
𝑍𝑍はハザード関数のスケーリング係数として扱い, 大きな𝑍𝑍の値は高いハザード(脆弱な部品の場合)に関連付けられる. 一方, 小さな𝑍𝑍の値は低いハザード (頑健な部品の場合)に関連付けられる.

いま𝑓𝑓𝑍𝑍 (𝑧𝑧)を𝑍𝑍の確率密度関数とし, パラメタベクトル𝜼𝜼によって特徴づけられるとする. 一般的によく用いられる𝑍𝑍のモデルは, ガンマ分布, 対数正規分布, 安定分布, 逆ガウス分布である(Whitmore
and Lee 1991; Joe 1993; Hougaard 1995; Duchateau and Janssen 2008; Ha et al. 2017; Piancastelli
et al. 2020). 一方, 数学的な利便性から我々が提案するのは, 次の密度関数である.

1.

２パラメタガンマ分布； 𝑍𝑍~𝐺𝐺𝐺𝐺𝑚𝑚𝑚𝑚𝑎𝑎(𝛼𝛼, 𝛽𝛽)𝑓𝑓𝑍𝑍 (𝑧𝑧) =2.

対数正規分布； 𝑍𝑍~𝐿𝐿𝐿𝐿(𝜇𝜇, 𝜎𝜎 2 )𝑓𝑓𝑍𝑍 (𝑧𝑧) =

3.

4.

1
𝑧𝑧
𝑧𝑧 𝛼𝛼−1 exp (− ) ,𝛤𝛤(𝛼𝛼)𝛽𝛽 𝛼𝛼
𝛽𝛽1𝑧𝑧𝑧𝑧√2𝜋𝜋

exp {−

切断正規分布； 𝑍𝑍~𝑇𝑇𝑇𝑇(𝜇𝜇, 𝜎𝜎 2 )𝑓𝑓𝑍𝑍 (𝑧𝑧) =

√2

𝜎𝜎√𝜋𝜋

(log⁡(𝑧𝑧) − 𝜇𝜇)2
},
2𝜎𝜎 2

exp {−(𝑧𝑧 − 𝜇𝜇)2
},
2𝜎𝜎 2

𝑧𝑧 > 0,𝛼𝛼 > 0,𝑧𝑧 > 0,− ∞ < 𝜇𝜇 < ∞,𝑧𝑧 > 𝜇𝜇,

折り畳み正規分布(Leone et al. 1961)； 𝑍𝑍~𝐹𝐹𝐹𝐹(𝜇𝜇, 𝜎𝜎 2 )𝑓𝑓𝑍𝑍 (𝑧𝑧) =1𝜎𝜎√2𝜋𝜋[ exp {−

𝛽𝛽 > 0− ∞ < 𝜇𝜇 < ∞,(𝑧𝑧 + 𝜇𝜇)2
(𝑧𝑧 − 𝜇𝜇)2
} + exp {−
} ],
2
2𝜎𝜎
2𝜎𝜎 2

𝜎𝜎 > 0𝜎𝜎 > 0𝑧𝑧 > 0, 𝜇𝜇 ≥ 0, 𝜎𝜎 > 0

上記４つのモデルは次の分布特性を持つ. 期待値は𝐸𝐸 𝐺𝐺𝐺𝐺𝐺𝐺𝐺𝐺𝐺𝐺 (𝑍𝑍) = 𝛼𝛼𝛼𝛼, 𝐸𝐸 𝐿𝐿𝐿𝐿 (𝑍𝑍) = exp(𝜇𝜇 +𝜎𝜎 2 /2), 𝐸𝐸 𝑇𝑇𝑇𝑇 (𝑍𝑍) = 𝜇𝜇 + (2/𝜋𝜋)1/2 𝜎𝜎, 𝐸𝐸 𝐹𝐹𝐹𝐹 (𝑍𝑍) = (2/𝜋𝜋)1/2 𝜎𝜎 × exp{−𝜇𝜇 2 /(2𝜎𝜎 2 )} + 𝜇𝜇{1 − 2Φ(−𝜇𝜇/𝜎𝜎)}となる. 切断正規分布と折り畳み正規分布は, 𝜇𝜇 = 0場合にのみ同等である. 従って, この２つのモデルが一致するのは, ゼロ切断モデルとゼロ折り畳みモデルの場合のみである. 折り畳み正規分布で𝜇𝜇 > 3𝜎𝜎の場合, 正規分布とほぼ同等であるが, 切断正規分布と大きく異なる. しかし, 正
規分布自体は負の値をとるため, フレイルティ分布とは考えられない. 切断正規分布と折り畳み正規分布は多くの文献で考慮されている(Cohen 1951; 1961; 1991; Leone et al. 1961; Sundberg 1974;
Lodhi, et al. 2021; Zeng and Gui 2021)にもかかわらず, フレイルティモデルの文献においてほとん
ど考慮されていない. また４つのモデルは Marshall and Olkin(1988)では考慮されていない. 彼ら
は 1 パラメタガンマ分布, 安定分布, 負の二項分布, および切断二項分布を考慮した.

21

21

2.2 多変量フレイルティ-コピュラモデル
本節では多変量故障時間(𝑇𝑇1 , . . . , 𝑇𝑇𝑘𝑘 )に対して「フレイルティ-コピュラモデル」を定義・提案する.

Oakes (1989)に従い, すべての故障時間はフレイルティ値𝑍𝑍 = 𝑧𝑧の影響を受けるとする. すな

わち, 𝑇𝑇𝑗𝑗 の𝑍𝑍 = 𝑧𝑧を所与としたときの条件生存関数を次のように定義する
𝑆𝑆𝑇𝑇𝑗𝑗|𝑍𝑍 (𝑡𝑡𝑗𝑗 |𝑧𝑧) ≡ Pr(𝑇𝑇𝑗𝑗 > 𝑡𝑡𝑗𝑗 |𝑧𝑧) = 𝑆𝑆0𝑗𝑗 (𝑡𝑡𝑗𝑗 ) 𝑧𝑧 ,ここで𝑆𝑆0𝑗𝑗 (⋅)は任意の連続型生存関数である. フレイルティ項が正の確率変数であるから, 任意の𝑧𝑧 > 0に対して𝑆𝑆𝑇𝑇𝑗𝑗|𝑍𝑍 (𝑡𝑡𝑗𝑗 |𝑧𝑧)は生存関数となる.

Marshall and Olkin(1988)に従い, 条件付き生存関数をコピュラ𝐶𝐶𝜃𝜃 : [0, 1]𝑘𝑘 → [0, 1]を用いて𝑆𝑆𝑇𝑇1,...,𝑇𝑇𝑘𝑘|𝑍𝑍 (𝑡𝑡1 , . . . , 𝑡𝑡𝑘𝑘 |𝑧𝑧) ≡ Pr(𝑇𝑇1 > 𝑡𝑡1 , . . . , 𝑇𝑇𝑘𝑘 > 𝑡𝑡𝑘𝑘 |𝑧𝑧) = 𝐶𝐶𝜃𝜃 {𝑆𝑆𝑇𝑇1|𝑍𝑍 (𝑡𝑡1 |𝑧𝑧), . . . , 𝑆𝑆𝑇𝑇𝑘𝑘|𝑍𝑍 (𝑡𝑡𝑘𝑘 |𝑧𝑧)},(1)と定義する. ここで𝜃𝜃はコピュラのパラメタである. 式 (1) は, 任意のコピュラに対して, 多変量生存関数となることは明らかである. Marshall and Olkin(1988)は, 次のより一般的なモデル
𝐶𝐶𝜃𝜃 {𝑆𝑆𝑇𝑇1|𝑍𝑍 (𝑡𝑡1 |𝑧𝑧1 ), . . . , 𝑆𝑆𝑇𝑇𝑘𝑘|𝑍𝑍 (𝑡𝑡𝑘𝑘 |𝑧𝑧𝑘𝑘 )},を考えた. ここで(𝑧𝑧1 , . . . , 𝑧𝑧𝑘𝑘 )多変量フレイルティ項である. 我々はこの一般的なモデルではなく,
式 (1) の 「 共 有 」 フ レ イ ル テ ィ 𝑧𝑧1 = ⋯ = ⁡ 𝑧𝑧𝑘𝑘 = 𝑧𝑧 の 場 合 に 焦 点 を 当 て る . 独 立 コ ピ ュ ラ𝐶𝐶𝜃𝜃 (𝑢𝑢1 , . . . , 𝑢𝑢𝑘𝑘 ) = ∏𝑘𝑘𝑗𝑗=1 𝑢𝑢𝑗𝑗 の場合, モデル(1)は通常のフレイルティモデルに帰着する (Oakes
1989).

Emura et al.(2017, 2019)は𝑍𝑍 = 𝑧𝑧をクラスター内の全ての部品（対象者）で共有される値とし, ２次元コピュラ𝐶𝐶𝜃𝜃 に対してモデル(1)を考え, これを結合フレイルティ-コピュラ（joint frailty-copula）モデルと呼んでいる. 彼らは−d{log𝑆𝑆0𝑗𝑗 (𝑡𝑡)}/𝑑𝑑𝑑𝑑をスプライン関数によってモデリングしている. また𝑆𝑆0𝑗𝑗 (. )をワイブル分布でモデリングする場合も Wu et al.(2020)および Shinohara et al.(2020)で提案されている. このようなクラスター化されたデータの解析においては, 𝑍𝑍 = 𝑧𝑧の値は未観測であるものの, クラスターを観測することにより, その値を類推することが出来る. 従って, 興味のある相関
構造は𝑍𝑍 = 𝑧𝑧を所与としたときの条件付きケンドール順位相関係数𝜏𝜏𝜃𝜃|𝑍𝑍 = 4 ∬ 𝑆𝑆𝑇𝑇1𝑇𝑇2 |𝑍𝑍 (𝑡𝑡1 , 𝑡𝑡2 |𝑧𝑧)𝑑𝑑𝑆𝑆𝑇𝑇1𝑇𝑇2|𝑍𝑍 (𝑡𝑡1 , 𝑡𝑡2 |𝑧𝑧) − 1 = 4 ∬𝐶𝐶𝜃𝜃 (𝑣𝑣, 𝑤𝑤)𝑑𝑑𝑑𝑑𝜃𝜃 (𝑣𝑣, 𝑤𝑤) − 1[0,1]×[0,1]である.
しかしながら, 我々の本稿における目的はクラスター化されたデータの解析ではない. つまり
我々は, 𝑍𝑍は未観測であるとし, これをその分布に関する積分で除外し

22

22

∞

𝑃𝑃(𝑇𝑇1 > 𝑡𝑡1 , . . . , 𝑇𝑇𝑘𝑘 > 𝑡𝑡𝑘𝑘 ) ≡ 𝑆𝑆𝑇𝑇1 ... 𝑇𝑇𝑘𝑘 (𝑡𝑡1 , . . . , 𝑡𝑡𝑘𝑘 ) = ∫ 𝑆𝑆𝑇𝑇1…𝑇𝑇𝑘𝑘|𝑍𝑍 (𝑡𝑡1 , . . . , 𝑡𝑡𝑘𝑘 |𝑧𝑧)𝑓𝑓𝑍𝑍 (𝑧𝑧) 𝑑𝑑𝑑𝑑0を得る. 上式の積分内の条件付き生存関数を式(1)で指定し次の定義を与える.

定義(多変量フレイルティ-コピュラモデル): 次式の多変量生存関数の族を定義する.
∞

𝑃𝑃(𝑇𝑇1 > 𝑡𝑡1 , . . . , 𝑇𝑇𝑘𝑘 > 𝑡𝑡𝑘𝑘 ) ≡ 𝑆𝑆𝑇𝑇1 ... 𝑇𝑇𝑘𝑘 (𝑡𝑡1 , . . . , 𝑡𝑡𝑘𝑘 ) = ∫ 𝐶𝐶𝜃𝜃 { 𝑆𝑆01 (𝑡𝑡1 ) 𝑧𝑧 , … , 𝑆𝑆0𝑘𝑘 (𝑡𝑡𝑘𝑘 )𝑧𝑧 }𝑓𝑓𝑍𝑍 (𝑧𝑧) 𝑑𝑑𝑑𝑑

（２）0

ここで𝑓𝑓𝑍𝑍 (. )は正値確率変数の密度関数, 𝐶𝐶𝜃𝜃 (. )は𝑘𝑘変量コピュラ, および𝑆𝑆0𝑗𝑗 (. ), 𝑗𝑗 ∈ {1, . . . , 𝑘𝑘}は任意の生存関数とする.

提案するモデル族(2)は, コピュラによって与えられる相関構造を共有フレイルティの効果でより一般化した族と見做すことができる. 式(2)は Marshall and Olkin(1988)が与えたモデル族の部分族であるが, 彼らは共有フレイルティ分布とコピュラの選択について十分に議論していない.
Charpentier et al.(2014)は, 極値コピュラ族の下で式(2)と類似の族を定義したが, クレイトン型コピ
ュラや FGM 型コピュラなど, 極値コピュラ族に属さない多くの有用なコピュラが存在する.
式(2)のフレイルティ-コピュラモデルは広い分布族であり, 様々なモデルを含む. この名称「フレイルティ-コピュラ」に従い, この族に属する個々のモデルの名称を系統的に与えることが出来る.
たとえば, フレイルティがガンマ分布, コピュラがクレイトン型コピュラで決まるモデルは「ガンマ-ク
レイトンモデル」と呼ぶことにする.
式（２）の積分をどのように計算するかは最も重要な問題であるが, Marshall and Olkin(1988)では議論が不十分であった. 以下の我々の議論では, その問題の解答のいくつかを与える.
2.3 提案モデルの一般的性質
数学的利便性の観点から, フレイルティ-コピュラモデルの式(2)を次のように書き換える
𝑆𝑆𝑇𝑇1 ... 𝑇𝑇𝑘𝑘 (𝑡𝑡1 , . . . , 𝑡𝑡𝑘𝑘 ) = 𝐻𝐻𝜃𝜃,𝜼𝜼 {𝛬𝛬01 (𝑡𝑡1 ), . . . , 𝛬𝛬0𝑘𝑘 (𝑡𝑡𝑘𝑘 )}.

(3)

ここで𝛬𝛬0𝑗𝑗 (𝑡𝑡𝑗𝑗 ) = − log 𝑆𝑆0𝑗𝑗 (𝑡𝑡𝑗𝑗 )は累積ハザード関数, 関数𝐻𝐻𝜃𝜃,𝜼𝜼 は後に例を通して議論する. 式(2)はモデルの構築法を明示する一方, 式(3)は統計的モデリングに便利である. たとえば, 累積ハザード関数にワイブル回帰モデルを設定する場合, 𝛬𝛬0𝑗𝑗 (𝑡𝑡𝑗𝑗 ) = 𝜆𝜆𝑗𝑗 𝑡𝑡𝑗𝑗𝜈𝜈 exp⁡(𝛽𝛽𝑗𝑗 𝑥𝑥𝑗𝑗 )とし, ここで尺度パラメタ

𝜆𝜆𝑗𝑗 > 0, 形状パラメタ𝜈𝜈 > 0, 共変量𝑥𝑥𝑗𝑗 をおいた. 共変量が無い場合, 𝛬𝛬0𝑗𝑗 (𝑡𝑡𝑗𝑗 ) = 𝜆𝜆𝑗𝑗 𝑡𝑡𝑗𝑗𝜈𝜈 とおく.

我々の興味は, 式(2)の積分を計算する方法や, 𝐻𝐻𝜃𝜃,𝜼𝜼 の式を求める方法である. コピュラやフレ

イルティ分布を上手く選び, 式(2)または式(3)が簡潔な形で表されれば, そのモデルは実用的であ23

23

ろう. この問題に対し Charpentier et al.(2014)のアイデアを利用し, 式(2)を∞𝑆𝑆𝑇𝑇1 … 𝑇𝑇𝑘𝑘 (𝑡𝑡1 , … , 𝑡𝑡𝑘𝑘 ) = ∫ exp[−ℓ𝜃𝜃 {𝑧𝑧𝑧𝑧01 (𝑡𝑡1 ), . . . , z𝛬𝛬0𝑘𝑘 (𝑡𝑡𝑘𝑘 )}]𝑓𝑓𝑍𝑍 (𝑧𝑧) 𝑑𝑑𝑑𝑑,0としてみる. ここでℓ𝜃𝜃 = − log(𝐶𝐶𝜃𝜃 )は依存関数（dependence function）と呼ばれる. この依存関数が条件ℓ𝜃𝜃 (𝑧𝑧𝑧𝑧1 , … , 𝑧𝑧𝑧𝑧𝑘𝑘 ) = 𝑧𝑧ℓ𝜃𝜃 (𝑡𝑡1 , … , 𝑡𝑡𝑘𝑘 )を満足すると, 次のように積分が計算できる
𝑆𝑆𝑇𝑇1 … 𝑇𝑇𝑘𝑘 (𝑡𝑡1 , … , 𝑡𝑡𝑘𝑘 ) = 𝐿𝐿𝜼𝜼 [ℓ𝜃𝜃 {𝛬𝛬01 (𝑡𝑡1 ), … , 𝛬𝛬0𝑘𝑘 (𝑡𝑡𝑘𝑘 )}]= 𝐿𝐿𝜼𝜼 [ − log 𝐶𝐶𝜃𝜃 {exp⁡(−𝛬𝛬01 (𝑡𝑡1 )), … , exp⁡(−𝛬𝛬0𝑘𝑘 (𝑡𝑡𝑘𝑘 ))}],

ここで∞𝐿𝐿𝜼𝜼 (𝑠𝑠)⁡= ∫ exp(−𝑧𝑧𝑧𝑧)𝑓𝑓𝑍𝑍 (𝑧𝑧|𝜼𝜼) 𝑑𝑑𝑑𝑑
0はラプラス変換である. 残念ながら, この条件は常に満足されるものでなく, また積分が計算出来るための必要条件でもない. この積分の計算問題は, ３節～５節でさらに個別に議論していく.
提案するモデル（２）の周辺分布を見てみよう. 周辺生存関数は∞

𝑆𝑆𝑇𝑇𝑗𝑗 (𝑡𝑡𝑗𝑗 ) = ∫ 𝑆𝑆0𝑗𝑗 (𝑡𝑡𝑗𝑗 ) 𝑧𝑧 𝑓𝑓𝑍𝑍 (𝑧𝑧|𝜼𝜼) 𝑑𝑑𝑑𝑑 = 𝐿𝐿𝜼𝜼 [ − log⁡𝑆𝑆0𝑗𝑗 (𝑡𝑡𝑗𝑗 ) ]ǡ0

𝑗𝑗 = 1, 2, . . . , 𝑘𝑘,(4)となるので, これはフレイルティ分布に依存する. 式(4)を利用して, 式(2)を次のように書き換える
𝑆𝑆𝑇𝑇1 ... 𝑇𝑇𝑘𝑘 (𝑡𝑡1 , . . . , 𝑡𝑡𝑘𝑘 ) = 𝐶𝐶𝜃𝜃,𝜼𝜼 {𝑆𝑆𝑇𝑇1 (𝑡𝑡1 ), . . . , 𝑆𝑆𝑇𝑇𝑘𝑘 (𝑡𝑡𝑘𝑘 )},ここで𝐶𝐶𝜃𝜃,𝜼𝜼 はコピュラである. 明らかに𝐶𝐶𝜃𝜃,𝜼𝜼 ≠ 𝐶𝐶𝜃𝜃 であり, 実際∞−1𝐶𝐶𝜃𝜃,𝜼𝜼 (𝑢𝑢1 , . . . , 𝑢𝑢𝑘𝑘 ) = ∫ 𝐶𝐶𝜃𝜃 [exp { − 𝑧𝑧𝐿𝐿−1𝜼𝜼 (𝑢𝑢1 )},…, exp { − 𝑧𝑧𝐿𝐿𝜼𝜼 (𝑢𝑢𝑘𝑘 )}]𝑓𝑓𝑍𝑍 (𝑧𝑧|𝜼𝜼) 𝑑𝑑𝑑𝑑0である. 上式から次の基本的性質が導かれる.
(i)
(ii)
(iii)
(iv)

𝐶𝐶𝜃𝜃,𝜼𝜼 は𝑆𝑆0𝑗𝑗 (⋅)または𝛬𝛬0𝑗𝑗 (⋅)に依存しない.

−1
𝐶𝐶𝜃𝜃 (𝑢𝑢1 , . . . , 𝑢𝑢𝑘𝑘 ) = ∏𝑘𝑘𝑗𝑗=1 𝑢𝑢𝑗𝑗 の場合, 𝐶𝐶𝜃𝜃,𝜼𝜼 = 𝐿𝐿𝜼𝜼 [𝐿𝐿−1𝜼𝜼 (𝑢𝑢1 ) + ⋯ + 𝐿𝐿𝜼𝜼 (𝑢𝑢𝑘𝑘 )]である.

𝐶𝐶𝜃𝜃 = min(𝑢𝑢𝑖𝑖 ) の場合, 𝐶𝐶𝜃𝜃,𝜼𝜼 = min(𝑢𝑢𝑖𝑖 )である.
𝑖𝑖

𝑖𝑖

𝑓𝑓𝑍𝑍 (𝑧𝑧|𝜼𝜼)が点𝑧𝑧 = 1に退化するとき, 𝐶𝐶𝜃𝜃,𝜼𝜼 = 𝐶𝐶𝜃𝜃 である.
24

24

(5)

性質(i)は, 𝐶𝐶𝜃𝜃,𝜼𝜼 と𝑆𝑆0𝑗𝑗 (⋅)が別々にモデリングできることを示す. 性質(ii)は, 𝐶𝐶𝜃𝜃,𝜼𝜼 が純粋なフレイル

ティモデルに帰着することを示す. 性質(iii)は, コピュラが上限に達した場合, フレイルティ分布は不定（任意である）ことを示す. 性質(iv) は, 𝐶𝐶𝜃𝜃,𝜼𝜼 が純粋なコピュラ𝐶𝐶𝜃𝜃 に帰着することを示す. 𝐶𝐶𝜃𝜃,𝜼𝜼 のより具体的な性質は, コピュラとフレイルティ分布を実際に与えた３節～５節で調べる.
2.4 ケンドール順位相関係数

故障時間(𝑇𝑇1 , . . . , 𝑇𝑇𝑘𝑘 )の間の相関の強さを表す指標であるケンドール順位相関係数を考える. それ に は , フ レイ ル テ ィ - コ ピ ュ ラ モ デ ル (2) の 下 で 適 当 な ペ ア (𝑇𝑇𝑖𝑖 , 𝑇𝑇𝑗𝑗 ) を 選び , ２変 量 コ ピ ュ ラ𝐶𝐶𝜃𝜃,𝜼𝜼 (1, . . . , 𝑢𝑢𝑖𝑖 , . . . , 𝑢𝑢𝑗𝑗 , . . . , 1)に関する順位相関係数を考えれば良い. 従って, 𝑘𝑘 = 2の場合を議論すれば十分である. Nelsen(2006)に習い, 𝑘𝑘 = 2の場合のケンドール順位相関係数は

𝜏𝜏𝜃𝜃,𝜼𝜼 = 4 ∬ 𝑆𝑆𝑇𝑇1𝑇𝑇2 (𝑡𝑡1 , 𝑡𝑡2 )𝑑𝑑𝑆𝑆𝑇𝑇1𝑇𝑇2 (𝑡𝑡1 , 𝑡𝑡2 ) − 1 = 4 ∬[0,1]2 𝐶𝐶𝜃𝜃,𝜼𝜼 (𝑣𝑣, 𝑤𝑤)𝑑𝑑𝑑𝑑𝜃𝜃,𝜼𝜼 (𝑣𝑣, 𝑤𝑤) − 1,(6)である. ここで𝐶𝐶𝜃𝜃,𝜼𝜼 (𝑣𝑣, 𝑤𝑤)は0 ≤ 𝑣𝑣, 𝑤𝑤 ≤ 1で定義される二変量コピュラである. 𝐶𝐶𝜃𝜃,𝜼𝜼 が特定のコピュラになる場合, 積分を計算し, 𝜏𝜏𝜃𝜃,𝜼𝜼 が簡潔な式で表示できる. しかしながら, 𝜏𝜏𝜃𝜃,𝜼𝜼 の式は一般には求めることが難しく, 数値計算を要する. Marshall and Olkin(1988)に従い, 我々も次のシミュレーショ
ンアルゴリズムを採用して, 𝜏𝜏𝜃𝜃,𝜼𝜼 の値を近似することを推奨する:

ケンドール順位相関τθ,𝛈𝛈を求めるアルゴリズム:

ステップ 1: 𝑍𝑍~𝑓𝑓𝑍𝑍 (𝑧𝑧|𝜼𝜼)を生成.

ステップ 2: (𝑈𝑈1 , 𝑈𝑈2 )~𝐶𝐶𝜃𝜃 (𝑣𝑣, 𝑤𝑤)を生成.
⁄

⁄

−1
−1
ステップ 3: 𝑇𝑇1 = 𝑆𝑆01
(𝑈𝑈11 𝑍𝑍 );⁡𝑇𝑇2 = 𝑆𝑆02
(𝑈𝑈21 𝑍𝑍 )とおく.

ステップ 4 : ステップ 1～3 を𝑅𝑅回繰り返し, 次を得る
𝑅𝑅 −1
𝑗𝑗
𝑗𝑗
𝜏𝜏̂𝜃𝜃,𝜼𝜼 = ( ) ∑sgn{(𝑇𝑇1𝑖𝑖 − 𝑇𝑇1 )(𝑇𝑇2𝑖𝑖 − 𝑇𝑇2 )}.
2
1≤𝑖𝑖<𝑗𝑗≤𝑅𝑅

上記𝜏𝜏̂𝜃𝜃,𝜼𝜼 は𝜏𝜏𝜃𝜃,𝜼𝜼 の不偏推定量である.

3. ガンマ-クレイトンモデル
クレイトン型コピュラ(Clayton 1978)は, 次のように定義される
𝐶𝐶𝜃𝜃 (𝑢𝑢1 , . . . , 𝑢𝑢𝑘𝑘 ) = (𝑢𝑢1−𝜃𝜃 + ⋯ + 𝑢𝑢𝑘𝑘−𝜃𝜃 − 𝑘𝑘 + 1)−1/𝜃𝜃 ,𝜃𝜃 > 0.

クレイトン型コピュラは, 生存時間解析において最も広く用いられているコピュラである(Rivest andWells 2001; Rotolo et al. 2013; Emura et al. 2014; Emura and Chen 2016; Moradian et al. 2019;
Kwon et al. 2021; Emura et al. 2020; Emura, Sofeu and Rondeau 2021). 統計的工程管理にも用いられる（Long et al. 2014; Sun et al. 2020; Huang and Emura 2021; Emura, Lai and Sun 2021; Kim
25

25

et al. 2021; 江村 2021）.
クレイトン型コピュラの下で, 提案するモデル(2)は次のようになる∞−

1

𝑆𝑆𝑇𝑇1 ... 𝑇𝑇𝑘𝑘 (𝑡𝑡1 , . . . , 𝑡𝑡𝑘𝑘 ) = ∫ {𝑆𝑆01 (𝑡𝑡1 )−𝜃𝜃𝜃𝜃 + ⋯ + 𝑆𝑆0𝑘𝑘 (𝑡𝑡𝑘𝑘 )−𝜃𝜃𝜃𝜃 − 𝑘𝑘 + 1} 𝜃𝜃 𝑓𝑓𝑍𝑍 (𝑧𝑧) 𝑑𝑑𝑑𝑑.
0

上記積分はどのようなフレイルティ分布でも陽に書けない. ガンマフレイルティ分布の場合を下で解説する. 対数正規フレイルティ分布の場合の説明は Wang and Emura(2021)に譲る.
𝑍𝑍~𝐺𝐺𝐺𝐺𝑚𝑚𝑚𝑚𝑎𝑎(𝛼𝛼, 𝛽𝛽)の下で生成されるガンマ-クレイトンモデルは次式となる∞−

1

𝑆𝑆𝑇𝑇1…𝑇𝑇𝑘𝑘 (𝑡𝑡1 , . . . , 𝑡𝑡𝑘𝑘 ) = ∫{𝑆𝑆01 (𝑡𝑡1 )−𝜃𝜃𝜃𝜃 + ⋯ + 𝑆𝑆0𝑘𝑘 (𝑡𝑡𝑘𝑘 )−𝜃𝜃𝜃𝜃 − 𝑘𝑘 + 1} 𝜃𝜃 ⁡01𝑧𝑧𝑧𝑧 𝛼𝛼−1 exp (− ) 𝑑𝑑𝑑𝑑.
𝛼𝛼𝛤𝛤(𝛼𝛼)𝛽𝛽
𝛽𝛽

ここで𝜃𝜃 > 0, 𝛼𝛼 > 0,⁡ 𝛽𝛽 > 0である. 上記の同時生存関数を次のように書き換える
𝑆𝑆𝑇𝑇1…𝑇𝑇𝑘𝑘 (𝑡𝑡1 , . . . , 𝑡𝑡𝑘𝑘 ) = 𝐻𝐻𝜃𝜃,𝛼𝛼,𝛽𝛽 {𝛬𝛬01 (𝑡𝑡1 ), … , 𝛬𝛬0𝑘𝑘 (𝑡𝑡𝑘𝑘 )},

ここで

1

1

− ( −log𝑤𝑤 )𝛼𝛼−1𝐻𝐻𝜃𝜃,𝛼𝛼,𝛽𝛽 (𝑎𝑎1 , . . . , 𝑎𝑎𝑘𝑘 ) = ∫ (𝑤𝑤 −𝑎𝑎1 𝜃𝜃𝜃𝜃 + ⋯ + 𝑤𝑤 −𝑎𝑎𝑘𝑘𝜃𝜃𝜃𝜃 − 𝑘𝑘 + 1) 𝜃𝜃0𝛤𝛤(𝛼𝛼)(7)

𝑑𝑑𝑑𝑑 ,⁡⁡ 𝑎𝑎𝑗𝑗 > 0である. この積分は陽に書けないが, 𝐻𝐻𝜃𝜃,𝛼𝛼,𝛽𝛽 (. )を数値計算することは容易である. モデル(7)は, ガ

ンマフレイルティモデルの一般化と見なすことができる. すなわち, 𝜃𝜃 ↓ 0で退化したモデルは
𝐻𝐻0,𝛼𝛼,𝛽𝛽 (𝑎𝑎1 , . . . , 𝑎𝑎𝑘𝑘 ) = (1 + 𝛽𝛽 ∑

𝑘𝑘

𝑗𝑗=1

−𝛼𝛼

𝑎𝑎𝑗𝑗 )となる. 式(7)の周辺生存関数は𝑆𝑆𝑇𝑇𝑗𝑗 (𝑡𝑡𝑗𝑗 ) = {1 + 𝛽𝛽𝛽𝛽0𝑗𝑗 (𝑡𝑡𝑗𝑗 )}−𝛼𝛼 , 𝑗𝑗 = 1, . . . , 𝑘𝑘, となる. したがって
𝑆𝑆𝑇𝑇1...𝑇𝑇𝑘𝑘 (𝑡𝑡1 , . . . , 𝑡𝑡𝑘𝑘 ) = 𝐶𝐶𝜃𝜃,𝛼𝛼 {𝑆𝑆𝑇𝑇1 (𝑡𝑡1 ), . . . , 𝑆𝑆𝑇𝑇𝑘𝑘 (𝑡𝑡𝑘𝑘 )},と書け, ここで

1−1/𝛼𝛼𝐶𝐶𝜃𝜃,𝛼𝛼 (𝑢𝑢1 , . . . , 𝑢𝑢𝑘𝑘 ) = ∫ {𝑤𝑤 (1−𝑢𝑢10は𝛽𝛽に依存しないコピュラである.

)𝜃𝜃

−1/𝛼𝛼

+ ⋯ + 𝑤𝑤 (1−𝑢𝑢𝑘𝑘)𝜃𝜃− 𝑘𝑘 + 1}−1/𝜃𝜃(−log𝑤𝑤)𝛼𝛼−1𝛤𝛤(𝛼𝛼)

𝑑𝑑𝑑𝑑 

二変量の場合, 同時生存関数はであり, ここで

𝑆𝑆𝑇𝑇1𝑇𝑇2 (𝑡𝑡1 , 𝑡𝑡2 ) = 𝐻𝐻𝜃𝜃,𝛼𝛼,𝛽𝛽 {𝛬𝛬01 (𝑡𝑡1 ), 𝛬𝛬02 (𝑡𝑡2 )};1−

1

𝐻𝐻𝜃𝜃,𝛼𝛼,𝛽𝛽 (𝑎𝑎1 , 𝑎𝑎2 ) = ∫ (𝑤𝑤 −𝑎𝑎1𝜃𝜃𝜃𝜃 + 𝑤𝑤 −𝑎𝑎2 𝜃𝜃𝜃𝜃 − 1) 𝜃𝜃 ⁡0

𝜃𝜃 > 0, 𝛼𝛼 > 0, 𝛽𝛽 > 0( −log𝑤𝑤 )𝛤𝛤(𝛼𝛼)

𝛼𝛼−1

𝑑𝑑𝑑𝑑 , ⁡𝑎𝑎1 > 0⁡,⁡⁡⁡ 𝑎𝑎2 > 0である. この関数は, 次の R コードを使用して簡単に計算できる.
26

26

H=function(w){(w^(-a1*theta*beta)+w^(-a2*theta*beta)-2+1)^(-1/theta)*(-log(w))^(alpha-1)/gamma(alpha)}integrate(H,0,1)$value
ガンマ分布のパラメタを𝛼𝛼 = 1/𝜂𝜂, 𝛽𝛽 = 𝜂𝜂 > 0とおくと, 平均値𝐸𝐸(𝑍𝑍) = 1, 分散𝑉𝑉𝑎𝑎𝑟𝑟(𝑍𝑍) = 𝜂𝜂となる.

このとき, フレイルティ分布は分散パラメタ𝜂𝜂によって特徴付けられる.

図１は𝑆𝑆𝑇𝑇1𝑇𝑇2 (𝑡𝑡1 , 𝑡𝑡2 )の等高線𝑆𝑆𝑇𝑇1𝑇𝑇2 (𝑡𝑡1 , 𝑡𝑡2 ) = 𝑝𝑝, 𝑝𝑝 ∈ {0.1, 0.2, … ,0.9}を表す. 𝜃𝜃が 1 から 8 に増えると等高線が垂直および水平になるのがわかる. また, 𝜂𝜂が 0.5 から 1 に増えると, 低い等高線が
縦軸・横軸に近づくことがわかる. よって, コピュラのパラメタとフレイルティ分布のパラメタが異なる働きをしている. 𝜃𝜃と𝜂𝜂が増えると, 完全正相関の等高線min⁡{𝑆𝑆𝑇𝑇1 (𝑡𝑡1 ), ⁡𝑆𝑆𝑇𝑇2 (𝑡𝑡2 )} = 𝑝𝑝に収束することも理解できる. 𝜃𝜃と𝜂𝜂が増えるとケンドール順位相関係数の値も増加する（図１）.

図 1:ガンマ-クレイトンモデルの生存関数の等高線.

27

27

4. グンベルコピュラ
グンベル型コピュラ(Gumbel 1960)は, 次のように定義される
𝐶𝐶𝜃𝜃 (𝑢𝑢1 , . . . , 𝑢𝑢𝑘𝑘 ) = exp[−{(− log 𝑢𝑢1 )𝜃𝜃+1 + ⋯ + (− log 𝑢𝑢𝑘𝑘 )𝜃𝜃+1 }1/( 𝜃𝜃+1 ) ],𝜃𝜃 ≥ 0.

グンベル型コピュラも生存時間解析で多用される(Peng et al. 2018; Wang et al. 2020; Kawakami etal. 2021; Sofeu et al. 2021). グンベル型コピュラの下で, 提案するモデル(2)は∞となり, ここで

𝑆𝑆𝑇𝑇1...𝑇𝑇𝑘𝑘 (𝑡𝑡1 , . . . , 𝑡𝑡𝑘𝑘 ) = ∫ exp(−𝑧𝑧𝐴𝐴𝜃𝜃 )𝑓𝑓𝑍𝑍 (𝑧𝑧) 𝑑𝑑𝑑𝑑
0𝑘𝑘

𝐴𝐴𝜃𝜃 = 𝐴𝐴(𝑡𝑡1 , . . . , 𝑡𝑡𝑘𝑘 ; 𝜃𝜃) = [ ∑{−log⁡𝑆𝑆0𝑗𝑗 (𝑡𝑡𝑗𝑗 )}𝜃𝜃+1 ]

1
𝜃𝜃+1

𝑘𝑘

= [ ∑ Λ𝜃𝜃+1
0𝑗𝑗 (𝑡𝑡𝑗𝑗 ) ]

1
𝜃𝜃+1

𝑗𝑗=1

𝑗𝑗=1

.

上記の積分は, 多くのフレイルティ分布の場合で陽に計算することができる. 積分の計算を以下のフレイルティ分布(ガンマ, 切断正規, 折り畳み正規)の場合で詳説する.
4.1 ガンマ-グンベルモデル:
𝑍𝑍~𝐺𝐺𝐺𝐺𝑚𝑚𝑚𝑚𝑎𝑎(𝛼𝛼, 𝛽𝛽), 𝛼𝛼 > 0,⁡ 𝛽𝛽 > 0とすると, ガンマ-グンベルモデルは次のように計算できる∞

𝑆𝑆𝑇𝑇1...𝑇𝑇𝑘𝑘 (𝑡𝑡1 , . . . , 𝑡𝑡𝑘𝑘 ) = ∫ exp(−𝑧𝑧𝐴𝐴𝜃𝜃 )01𝑧𝑧𝑧𝑧 𝛼𝛼−1 exp (− ) 𝑑𝑑𝑑𝑑⁡ = (1 + 𝛽𝛽𝐴𝐴𝜃𝜃 )−𝛼𝛼 .
𝛼𝛼𝛤𝛤( 𝛼𝛼 )𝛽𝛽𝛽𝛽周辺生存関数は𝑆𝑆𝑇𝑇𝑗𝑗 (𝑡𝑡𝑗𝑗 ) = {1 + 𝛽𝛽Λ 0𝑗𝑗 (𝑡𝑡𝑗𝑗 )}−𝛼𝛼 であるから, 上式は次のように表せる
𝑆𝑆𝑇𝑇1...𝑇𝑇𝑘𝑘 (𝑡𝑡1 , . . . , 𝑡𝑡𝑘𝑘 ) = (1 + [{𝑆𝑆𝑇𝑇1 (𝑡𝑡1 )−1/𝛼𝛼 − 1}𝜃𝜃+1

+ ⋯ + {𝑆𝑆𝑇𝑇𝑘𝑘 (𝑡𝑡𝑘𝑘 )−1/𝛼𝛼 − 1}𝜃𝜃+1 1/(𝜃𝜃+1)

]

−𝛼𝛼

)

.

したがって, このモデルのコピュラは次のようになる
1

𝜃𝜃+1𝜃𝜃+1 𝜃𝜃+1
−1/𝛼𝛼
−1/𝛼𝛼
− 1)
+ ⋯ + (𝑢𝑢𝑘𝑘− 1) }]𝐶𝐶𝜃𝜃,𝛼𝛼 (𝑢𝑢1 , . . . , 𝑢𝑢𝑘𝑘 ) = [1 + {(𝑢𝑢1

−𝛼𝛼

⁡ .

−1
[𝜙𝜙𝜃𝜃,𝛼𝛼 (𝑢𝑢1 ) +このコピュラは𝛽𝛽に依存しない. 上記コピュラはアルキメデス型𝐶𝐶𝜃𝜃,𝛼𝛼 (𝑢𝑢1 , . . . , 𝑢𝑢𝑘𝑘 ) = 𝜙𝜙𝜃𝜃,𝛼𝛼⋯ + 𝜙𝜙𝜃𝜃,𝛼𝛼 (𝑢𝑢𝑘𝑘 )]であり, 生成素𝜙𝜙𝜃𝜃,𝛼𝛼 (𝑡𝑡) = (𝑡𝑡 −1/𝛼𝛼 − 1)𝜃𝜃+1 を持つ(Nelsen 2006).
二変量の場合, 同時生存関数は

1

−𝛼𝛼

𝜃𝜃+1
𝜃𝜃+1
)
𝑆𝑆𝑇𝑇1𝑇𝑇2 (𝑡𝑡1 , 𝑡𝑡2 ) = ( 1 + 𝛽𝛽[ Λ𝜃𝜃+1
01 (𝑡𝑡1 ) + Λ 02 (𝑡𝑡2 ) ]28

28; 𝜃𝜃 ≥ 0,⁡⁡⁡𝛼𝛼 > 0,⁡⁡⁡𝛽𝛽 > 0.

であり, 二変量コピュラ は1

𝜃𝜃+1𝜃𝜃+1 𝜃𝜃+1
−1/𝛼𝛼
−1/𝛼𝛼
− 1)
+ (𝑢𝑢2− 1) }]𝐶𝐶𝜃𝜃,𝛼𝛼 (𝑢𝑢1 , 𝑢𝑢2 ) = [1 + {(𝑢𝑢1−𝛼𝛼.

Genest and MacKay (1986)の公式により, ケンドール順位相関係数は次のように計算される
1

𝜏𝜏𝜃𝜃, 𝛼𝛼 = 1 + 4 ∫0𝜙𝜙𝜃𝜃,𝛼𝛼 (𝑡𝑡)2
𝑑𝑑𝑑𝑑 = 1 −.
1𝑑𝑑𝜙𝜙𝜃𝜃,𝛼𝛼 (𝑡𝑡)/𝑑𝑑𝑑𝑑(𝜃𝜃 + 1) ( + 2)
𝛼𝛼フレイルティ分布を１パラメタに制約した, 𝐺𝐺𝐺𝐺𝑚𝑚𝑚𝑚𝑎𝑎(𝛼𝛼 = 1, 𝛽𝛽) = 𝐸𝐸𝑥𝑥𝑥𝑥(𝛽𝛽)や𝐺𝐺𝐺𝐺𝑚𝑚𝑚𝑚𝑎𝑎(𝛼𝛼 = 𝛾𝛾/2, 𝛽𝛽 =2も考慮できる. 制約𝛼𝛼 = 1/𝜂𝜂と𝛽𝛽 = 𝜂𝜂の下では, 二変量コピュラは
2) = 𝜒𝜒𝑑𝑑𝑑𝑑=𝛾𝛾

−1/𝜂𝜂
1
𝜃𝜃+1
𝜃𝜃+1 𝜃𝜃+1−𝜂𝜂−𝜂𝜂
+ (𝑢𝑢2 − 1) }]𝐶𝐶𝜃𝜃,𝜂𝜂 (𝑢𝑢1 , 𝑢𝑢2 ) = [1 + {(𝑢𝑢1 − 1)である. この𝐶𝐶𝜃𝜃,𝜂𝜂 は, BB1 コピュラ(Joe(1997)の P.150)である. Nelsen(2006)は, 2 パラメタのコピュラの例として挙げている(彼の書籍の例 4.22 と 4.26 を参照). コピュラ𝐶𝐶𝜃𝜃,𝜂𝜂 は𝜃𝜃 = 0のときクレイトンコピ

ュラに, 𝜂𝜂 → 0のときグンベルコピュラに帰着する. 𝜂𝜂がクレイトンコピュラからの乖離を示すため, 頑
健性の研究にも有用である（Emure, Sofeu, and Rondeau(2021)の Supplementary Material）.

BB1 コピュラは, Lu と Bhattacharyya(1990)によって与えられた 2 変量ワイブルモデルからも導かれる. Wang et al.（2020）はガンマ-グンベル分布を競合リスクデータ分析に用いている.
4.2 切断正規-グンベルモデル:
驚くべきことに, 𝑍𝑍~𝑇𝑇𝑇𝑇(𝜇𝜇, 𝜎𝜎 2 )の場合でも, 同時生存関数が次のように陽に書くことが出来る
𝑆𝑆𝑇𝑇1...𝑇𝑇𝑘𝑘 (𝑡𝑡1 , . . . , 𝑡𝑡𝑘𝑘 ) = 2exp (𝜎𝜎 2 2𝐴𝐴 − 𝜇𝜇𝐴𝐴𝜃𝜃 ) Φ(−𝜎𝜎𝐴𝐴𝜃𝜃 )
2 𝜃𝜃

⁡(8)

ここで, Φは𝑁𝑁(0, 1)の累積分布関数である. 導出は付録を参照されたい. 周辺生存関数は
𝜎𝜎 2
𝑆𝑆𝑇𝑇𝑗𝑗 (𝑡𝑡𝑗𝑗 ) = 2exp ( ⁡Λ20𝑗𝑗 (𝑡𝑡𝑗𝑗 ) − 𝜇𝜇⁡Λ 0𝑗𝑗 (𝑡𝑡𝑗𝑗 )) Φ(−𝜎𝜎⁡Λ 0𝑗𝑗 (𝑡𝑡𝑗𝑗 ))2となる. 式(8)は陽に書くことが出来たが, 対応するコピュラとケンドール順位相関係数は陽に書くことが出来ない（下記）. 単純な指数分布のモデルΛ 01 (𝑡𝑡) = 𝑡𝑡を考えると, 周辺生存関数は𝑆𝑆𝑇𝑇1 (𝑡𝑡) =2exp(𝜎𝜎 2 𝑡𝑡 2 /2 − 𝜇𝜇𝜇𝜇)Φ(−𝜎𝜎𝜎𝜎)となる.

制約𝐸𝐸(𝑍𝑍) = 𝜇𝜇 + √2/𝜋𝜋𝜎𝜎 = 1を与え, 基準化された１パラメタのフレイルティ分布を考える.

このフレイルティ分布のラプラス変換を

29

29

∞

𝐿𝐿𝜎𝜎2 (𝑠𝑠) = ∫1−√2/𝜋𝜋𝜎𝜎12exp(−𝑧𝑧𝑧𝑧) √ exp [−𝜎𝜎𝜋𝜋{𝑧𝑧−(1−√2/𝜋𝜋𝜎𝜎)}22𝜎𝜎2] 𝑑𝑑𝑑𝑑とおき, その逆関数を𝐿𝐿−1
𝜎𝜎2 (.)とする. このときコピュラは次のようになる𝐶𝐶𝜃𝜃, 𝜎𝜎2 (𝑢𝑢1 , 𝑢𝑢2 ) = 2exp [2
1
𝜎𝜎 2 −1𝜃𝜃+1 𝜃𝜃+1𝜃𝜃+1𝜃𝜃+1 𝜃𝜃+1− (1 − √2/𝜋𝜋𝜎𝜎){𝐿𝐿−1+ 𝐿𝐿−1]{𝐿𝐿𝜎𝜎2 (𝑢𝑢1 )𝜃𝜃+1 + 𝐿𝐿−1}}
2 (𝑢𝑢2 )
2 (𝑢𝑢1 )
2 (𝑢𝑢2 )
𝜎𝜎
𝜎𝜎
𝜎𝜎
2
1

𝜃𝜃+1
𝜃𝜃+1 𝜃𝜃+1+ 𝐿𝐿−1].
× Φ [ − 𝜎𝜎{𝐿𝐿−1}
𝜎𝜎 2 (𝑢𝑢1 )𝜎𝜎 2 (𝑢𝑢2 )

図２は𝑆𝑆𝑇𝑇1𝑇𝑇2 (𝑡𝑡1 , 𝑡𝑡2 )の等高線𝑆𝑆𝑇𝑇1𝑇𝑇2 (𝑡𝑡1 , 𝑡𝑡2 ) = 𝑝𝑝, 𝑝𝑝 ∈ {0.1, 0.2, … ,0.9}を表す. 𝜃𝜃が 1 から 8 に増えると等高線が垂直および水平になるのがわかる. また, 𝜎𝜎が 0.1 から 1 に増えると, 低い等高線が
縦軸・横軸に近づくことがわかる. よって, コピュラのパラメタとフレイルティ分布のパラメタが異なる働きをしている. 𝜃𝜃と𝜎𝜎が増えると, 完全正相関の等高線min⁡{𝑆𝑆𝑇𝑇1 (𝑡𝑡1 ), ⁡𝑆𝑆𝑇𝑇2 (𝑡𝑡2 )} = 𝑝𝑝に収束することも理解できる. 𝜃𝜃と𝜎𝜎が増えるとケンドール順位相関係数の値も増加する（図 2）.

図 2. 切断正規-グンベルモデルの生存関数の等高線.
30

30

4.3 折り畳み正規-グンベルモデル
さらに驚くべきことに, 折り畳み正規分布𝑍𝑍~𝐹𝐹𝐹𝐹(𝜇𝜇, 𝜎𝜎 2 )のような複雑な密度関数形式であっても,同時生存関数が次のように陽に書くことが出来る

𝜎𝜎 2
𝜇𝜇
𝜎𝜎 2
𝜇𝜇
𝑆𝑆𝑇𝑇1...𝑇𝑇𝑘𝑘 (𝑡𝑡1 , . . . , 𝑡𝑡𝑘𝑘 ) = exp ( 𝐴𝐴2𝜃𝜃 + 𝜇𝜇𝐴𝐴𝜃𝜃 ) Φ (− − 𝜎𝜎𝐴𝐴𝜃𝜃 ) + exp ( 𝐴𝐴2𝜃𝜃 − 𝜇𝜇𝐴𝐴𝜃𝜃 ) Φ ( − 𝜎𝜎𝐴𝐴𝜃𝜃 ).
2
2
𝜎𝜎
𝜎𝜎

導出の詳細については, 付録を参照されたい. ゼロで折り畳んだ(𝜇𝜇 = 0)場合,
𝑆𝑆𝑇𝑇1...𝑇𝑇𝑘𝑘 (𝑡𝑡1 , . . . , 𝑡𝑡𝑘𝑘 ) = 2exp (𝜎𝜎 2 2𝐴𝐴 ) Φ(−𝜎𝜎𝐴𝐴𝜃𝜃 )2 𝜃𝜃

となり, ゼロ切断正規-グンベルモデルの生存関数に等しい(4.2 節). モデルの性質は, 切断正規-グンベルモデルと同じように議論できる. 詳細は Wang and Emura(2021)を参照されたい.

5. ガンマ-FGM モデル
多変量 FGM コピュラ(Kotz et al. 2000)は, 次のように定義される
𝑘𝑘

𝑘𝑘

𝑖𝑖=1

𝑖𝑖=1

𝐶𝐶𝜃𝜃 (𝑢𝑢1 , . . . , 𝑢𝑢𝑘𝑘 ) = ∏ 𝑢𝑢𝑖𝑖 + 𝜃𝜃 ∏ 𝑢𝑢𝑖𝑖 (1 − 𝑢𝑢𝑖𝑖 ) ,⁡⁡⁡⁡⁡⁡⁡𝜃𝜃 ∈ [ − 1,1]⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡(9)

優れた数学的性質から, FGM コピュラは理論・応用上重要である(Domma and Giordano 2013; Shihand Emura 2018; 2020; 2021; Mohtashami-Borzadaran et al. 2019; Shih et al. 2019; de Oliveira etal. 2021; Kawakami et al. 2021; Ota and Kimura 2021; Šeliga et al. 2021). 提案するモデル(2)は次のようになる:∞

𝑆𝑆𝑇𝑇1...𝑇𝑇𝑘𝑘 (𝑡𝑡1 , . . . , 𝑡𝑡𝑘𝑘 ) = ∫ [ ∏𝑘𝑘𝑖𝑖=1 𝑆𝑆0𝑖𝑖 (𝑡𝑡𝑖𝑖 ) 𝑧𝑧 + 𝜃𝜃 ∏𝑘𝑘𝑖𝑖=1 𝑆𝑆0𝑖𝑖 (𝑡𝑡𝑖𝑖 ) 𝑧𝑧 {1 − 𝑆𝑆0𝑖𝑖 (𝑡𝑡𝑖𝑖 ) 𝑧𝑧 } ]𝑓𝑓𝑍𝑍 (𝑧𝑧) 𝑑𝑑𝑑𝑑.
0Marshall and Olkin(1988)は𝑘𝑘 = 2かつ𝑍𝑍~𝐺𝐺𝐺𝐺𝑚𝑚𝑚𝑚𝑎𝑎(𝛼𝛼, 𝛽𝛽)の場合, 積分が陽に書けることを示した.
以下で, 𝑘𝑘 ≥ 2の場合への一般化を考える. 𝑍𝑍~𝐺𝐺𝐺𝐺𝑚𝑚𝑚𝑚𝑎𝑎(𝛼𝛼, 𝛽𝛽)のとき,
∞

𝑘𝑘

𝑘𝑘

𝑆𝑆𝑇𝑇1 …𝑇𝑇𝑘𝑘 (𝑡𝑡1 , . . . , 𝑡𝑡𝑘𝑘 ) = ∫ [ ∏ 𝑆𝑆0𝑖𝑖 (𝑡𝑡𝑖𝑖 ) + 𝜃𝜃 ∏ 𝑆𝑆0𝑖𝑖 (𝑡𝑡𝑖𝑖 ) 𝑧𝑧 {1 − 𝑆𝑆0𝑖𝑖 (𝑡𝑡𝑖𝑖 ) 𝑧𝑧 } ]
0

𝑧𝑧

𝑖𝑖=1

𝑘𝑘

𝑖𝑖=1

= [ 1 + 𝛽𝛽 {∑ Λ 0𝑖𝑖 (𝑡𝑡𝑖𝑖 )} ]
𝑘𝑘

−𝛼𝛼

𝑖𝑖=1

𝑘𝑘

+ 𝜃𝜃 [ 1 + 𝛽𝛽 {∑ Λ 0𝑖𝑖 (𝑡𝑡𝑖𝑖 )} ]
𝑘𝑘

𝑖𝑖=1

−𝜃𝜃 ∑ [ 1 + 𝛽𝛽 { Λ0𝑖𝑖1 (𝑡𝑡𝑖𝑖1 ) + ∑ Λ 0𝑖𝑖 (𝑡𝑡𝑖𝑖 ) } ]𝑖𝑖1 =1
𝑘𝑘

+𝜃𝜃 ∑

𝑖𝑖=1

(−1)𝑗𝑗𝑗𝑗=2∑

𝑧𝑧 𝛼𝛼−1 exp(−𝑧𝑧/𝛽𝛽)𝑑𝑑𝑑𝑑𝛤𝛤(𝛼𝛼)𝛽𝛽 𝛼𝛼1≤𝑖𝑖1 <⋯<𝑖𝑖𝑗𝑗 ≤𝑘𝑘

−𝛼𝛼

−𝛼𝛼

𝑘𝑘

𝑗𝑗

𝑖𝑖=1

𝑚𝑚=1

[ 1 + 𝛽𝛽 { ∑ Λ 0𝑖𝑖 (𝑡𝑡𝑖𝑖 ) + ∑ Λ 0𝑖𝑖𝑚𝑚 (𝑡𝑡𝑖𝑖𝑚𝑚 ) } ]

31

31

−𝛼𝛼

⁡

である. 従って, ガンマ-FGM モデルは, 複雑な形式ではあるが陽に書くことが出来る. 周辺生存関数は𝑆𝑆𝑇𝑇𝑗𝑗 (𝑡𝑡𝑗𝑗 ) = {1 + 𝛽𝛽Λ 0𝑗𝑗 (𝑡𝑡𝑗𝑗 )}−𝛼𝛼 であるから, コピュラは次のようである
−



 k

 k

−1/
C ,  ,  ( u1 , ..., uk ) =)   +   1 −   ( 1 − ui −1/ )  
 1 −   ( 1 − ui=
 i 1=

 i1



k


−   um −1/ −  ( 1 − ui −1/ ) 
m 1=
i 1
=


k

k

−

−

+  ( − 1 ) j

=j 2

−

j

 k

−1/
) +  ( 1 − uim −1/ )   .
 1 −   ( 1 − uim 1
1i1  i j  k  =
=
 i1


このコピュラは文献には存在していないようである. コピュラが𝛽𝛽に依存しないことに注意する.
𝑘𝑘 = 2かつ𝛼𝛼 = 1/𝜂𝜂の場合, コピュラはC ,  ( u1 , u2 ) = ( u1− + u2 − − 1 )−1/ +  ( u1− + u2 − − 1 )−1/− ( 2u1− + u2 − − 2 )−1/ −  ( u1− + 2u2 − − 2 )−1/+ ( 2u1− + 2u2 − − 3 )−1/ .

上記のコピュラは Cook and Johnson (1986)が提案したものと同じであり, 上記のコピュラの導出はMarshall and Olkin(1988)と同じである. 𝜃𝜃 = 0でクレイトンコピュラに帰着し, 𝜂𝜂 → 0のとき FGM コピ

ュラに帰着する. このコピュラは, 他の一般化 FGM コピュラ(Bairamov and Bayramoglu 2013;
Huang et al. 2013; Domma and Giordano 2013; Hürlimann 2017; Saminger-Platz et al. 2020)と同様に, FGM コピュラの狭い相関の範囲を拡張することが出来る.
𝑘𝑘 = 2のとき, ガンマ-FGM 分布の同時生存関数は次のようである

𝑆𝑆𝑇𝑇1𝑇𝑇2 (𝑡𝑡1 , 𝑡𝑡2 ) = 𝐵𝐵𝛽𝛽 (𝑡𝑡1 , 𝑡𝑡2 )−𝛼𝛼 + 𝜃𝜃𝐵𝐵𝛽𝛽 (𝑡𝑡1 , 𝑡𝑡2 )−𝛼𝛼 − 𝜃𝜃{βΛ 01 (𝑡𝑡1 ) + 𝐵𝐵𝛽𝛽 (𝑡𝑡1 , 𝑡𝑡2 )}−𝛼𝛼−𝜃𝜃{βΛ02 (𝑡𝑡2 ) + 𝐵𝐵𝛽𝛽 (𝑡𝑡1 , 𝑡𝑡2 )}−𝛼𝛼 + 𝜃𝜃{2𝐵𝐵𝛽𝛽 (𝑡𝑡1 , 𝑡𝑡2 ) − 1}−𝛼𝛼 .

ここで𝐵𝐵𝛽𝛽 (𝑡𝑡1 , 𝑡𝑡2 ) = 1 + 𝛽𝛽{Λ 01 (𝑡𝑡1 ) + Λ 02 (𝑡𝑡2 )}とおいた. 生存関数𝑆𝑆𝑇𝑇1𝑇𝑇2 (𝑡𝑡1 , 𝑡𝑡2 )の等高線は Wangand Emura(2021)を参照されたい. さらに𝜃𝜃 = 0, 𝛼𝛼 = 1/𝜂𝜂, ⁡𝛽𝛽 = 𝜂𝜂とおくと,
𝑆𝑆𝑇𝑇1𝑇𝑇2 (𝑡𝑡1 , 𝑡𝑡2 ) = [1 + 𝜂𝜂{Λ 01 (𝑡𝑡1 ) + Λ 02 (𝑡𝑡2 )}]となり, これは共有ガンマフレイルティモデルに帰着する.

32

32

−1/𝜂𝜂

６.結論
本稿では, 多変量故障時間モデルの構築に関する２つの主要なモデル族である, コピュラモデル
とフレイルティモデルの統合を目指した. そこで, コピュラモデルとフレイルティモデルの両方を含
む, さらに広い多変量故障時間モデルの族である「多変量フレイルティ-コピュラモデル」を提案した. この提案したモデル族は, Marshall and Olkin(1988)によって提案されたものの部分族である.
本研究は, Marshall and Olkin(1988)が行った不完全な研究を補完・拡張するものと見做せる. 実
際, コピュラとフレイルティ分布の組み合わせにより, 様々な多変量分布を導き出すことに成功し,そのような新しい分布を生み出す手法を示した. さらに, 提案したモデル族に属する分布の名称
を系統的に与える手法も提案した. とりわけ魅力的な多変量分布は, ガンマ-クレイトンモデル, ガ
ンマ-グンベルモデル, 切断正規-グンベルモデル, 折り畳み正規-グンベルモデル, ガンマFGM モデルと命名されたものであった.
本稿で提案したいくつかの分布の実データへの当てはめや, これら分布下での統計的推測理論の構築は, 今後の重要な課題である. 生存時間データを扱うためにこれらの分布を用いることが想定されるが, 生存時間データは打ち切り・切断・競合リスクなどにより完全データを得ることが困難である（Klein and Moeschberger 2003; Dörre and Emura 2019; Dörre 2021; Emura and Ha
2021; Li et al. 2021; Wu et al. 2021）. とりわけ, 従属打ち切り（dependent censoring）の相関構造は同定が困難であり, コピュラを既知と扱う場合がある（Rivest and Wells 2001; Emura and Chen
2018; Emura and Hsu 2020）. 近年この従属打ち切りの相関構造が推定可能なモデルが提案され始めている（Deresa and Van Keilegom 2020). 同様に, 従属切断（dependent truncation）の相関構
造モデリング, とりわけパラメトリックモデリング（Emura and Konno 2012a,b; Emura and Pan 2020）にも提案したモデルが応用できる可能性がある.
謝辞
阿部俊弘先生には，日本統計研究所, 研究集会「統計的モデリングと統計的推測理論」において
本研究の口頭発表の機会および所報へ掲載の機会を頂いたことに心より感謝申し上げます．本研
究は, 台湾科学技術省(MOST 107-2118-M-008-003-MY3)からの助成金によって行われた.

33

33

付録: 数学的な導出
切断正規-グンベルモデルの導出:
切断正規分布⁡𝑍𝑍~𝑇𝑇𝑇𝑇(𝜇𝜇, 𝜎𝜎 2 )の下で, 以下の計算を行うと,

ST1 , ...,Tk ( t1 , ..., t k ) =  ST1 , ...,Tk |Z ( t1 , ..., t k | Z ) f Z ( z )dz−=  exp( − zA )




=


1



1



 ( z −  )2 
exp  −
dz
2 2 


2

2

1


exp  − 2 { 2 2 A z + ( z −  ) 2 } dz

 2


 (  −  2 A ) 2 −  2  1
 { z − (  −  2 A ) }2 
= 2 exp 
exp

−
dz ,

2 2
2 2

   2となる. これに変数変換𝑤𝑤 = {𝑧𝑧 − (𝜇𝜇 − 𝜎𝜎 2 𝐴𝐴𝜃𝜃 )}/𝜎𝜎を適用すると, 求める結果を得る
𝜎𝜎 2
𝑆𝑆𝑇𝑇1...𝑇𝑇𝑘𝑘 (𝑡𝑡1 , . . . , 𝑡𝑡𝑘𝑘 ) = 2exp ( 𝐴𝐴2𝜃𝜃 − 𝜇𝜇𝐴𝐴𝜃𝜃 ) Φ(−𝜎𝜎𝐴𝐴𝜃𝜃 ).
2

折り畳み正規-グンベルモデルの導出:
折り畳み正規分布𝑍𝑍~𝐹𝐹𝐹𝐹(𝜇𝜇, 𝜎𝜎 2 )の下で, 次式を計算する必要がある∞

𝑆𝑆𝑇𝑇1...𝑇𝑇𝑘𝑘 (𝑡𝑡1 , . . . , 𝑡𝑡𝑘𝑘 ) = ∫ exp(−𝑧𝑧𝐴𝐴𝜃𝜃 )01𝜎𝜎√2𝜋𝜋[ exp {−

(𝑧𝑧 + 𝜇𝜇)2
(𝑧𝑧 − 𝜇𝜇)2
}
+
exp
{−
} ] 𝑑𝑑𝑑𝑑.
2𝜎𝜎 2
2𝜎𝜎 2

先の切断正規分布の場合と同様の計算を行い, 次の積分を得る∞∫ exp(−𝑧𝑧𝐴𝐴𝜃𝜃 )01𝜎𝜎√2𝜋𝜋exp {−(𝑧𝑧 + 𝜇𝜇)2
𝜎𝜎 2 2
𝜇𝜇
}
𝑑𝑑𝑑𝑑
=
exp
(
𝐴𝐴 + 𝜇𝜇𝐴𝐴𝜃𝜃 ) Φ (− − 𝜎𝜎𝐴𝐴𝜃𝜃 ).
2𝜎𝜎 2
2 𝜃𝜃
𝜎𝜎

𝜇𝜇を−𝜇𝜇で置き換えると, もう一つの必要な積分を得る. よって

𝜎𝜎 2
𝜇𝜇
𝜎𝜎 2
𝜇𝜇
𝑆𝑆𝑇𝑇1...𝑇𝑇𝑘𝑘 (𝑡𝑡1 , . . . , 𝑡𝑡𝑘𝑘 ) = exp ( 𝐴𝐴2𝜃𝜃 + 𝜇𝜇𝐴𝐴𝜃𝜃 ) Φ (− − 𝜎𝜎𝐴𝐴𝜃𝜃 ) + exp ( 𝐴𝐴2𝜃𝜃 − 𝜇𝜇𝐴𝐴𝜃𝜃 ) Φ ( − 𝜎𝜎𝐴𝐴𝜃𝜃 ).
2
2
𝜎𝜎
𝜎𝜎

34

34

参考文献
Aalen OO (1994). Effects of frailty in survival analysis. Statistical Methods in Medical Research 3(3): 227-243.
Bairamov I, Bayramoglu K (2013). From the Huang-Kotz FGM distribution to Baker’s bivariate distribution,
Journal of Multivariate Analysis, 113: 106-115.
Balakrishnan N, Lai CD. (2009). Continuous Bivariate Distributions. Springer.
Campos E, Braekers R, de Souza DJ, Chaves LM. (2021). Factor copula models for right-censored clusteredsurvival data. Lifetime Data Analysis, 27: 499–535.
Charpentier A, Fougères AL, Genest C, Nešlehová JG. (2014). Multivariate archimax copulas. Journal ofMultivariate Analysis, 126, 118-136.
Clayton DG (1978). A model for association in bivariate life tables and its application to epidemiological studies
of familial tendency in chronic disease incidence. Biometrika. 65, 141-51
Cohen, A (1951). On estimating the mean and variance of singly truncated normal frequency distributions fromthe first three sample moments. Annals of the Institute of Statistical Mathematics 3, 37–44.
Cohen, A (1961). Tables for maximum likelihood estimates: Singly truncated and singly censored samples.
Technometrics 3, 535–541.
Cohen, A (1991). Truncated and Censored Samples: Theory and Applications, CRC Press.
Cook RD, Johnson ME. (1986). Generalized Burr-Pareto-logistic distributions with applications to a uranium
exploration data set. Technometrics, 28(2), 123-131.
Crowder MJ (2012). Multivariate Survival Analysis and Competing Risks, CRC Pressde Oliveira RP, de Oliveira Peres MV, dos Santos MR, Martinez EZ, Achcar JA (2021). A Bayesian inference
approach for bivariate Weibull distributions derived from Roy and Morgenstern methods. Statistics,
Optimization & Information Computing, 9(3), 529-554.
Deresa NW, Van Keilegom I (2020). A multivariate normal regression model for survival data subject todifferent types of dependent censoring. Computational Statistics & Data Analysis, 144, 106879.
Domma F, Giordano S (2013). A copula-based approach to account for dependence in stress-strength models.
Statistical Papers, 54(3), 807-826.
Dörre, A. (2021). Semiparametric likelihood inference for heterogeneous survival data under double truncation
based on a Poisson birth process. Japanese Journal of Statistics and Data Science, doi:10.1007/s42081-02100128-w.
Dörre A, Emura T (2019), Analysis of Doubly Truncated Data, An Introduction, JSS Research Series in
Statistics, Springer.
Duchateau L, Janssen P (2008). The Frailty Model, Springer.
Duchateau L, Janssen P, Lindsey P, Legrand C, Nguti R, Sylvester R (2002). The shared frailty model and thepower for heterogeneity tests in multicenter trials. Computational Statistics & Data Analysis 40(3), 603-620.
Emura T, Chen YH (2016). Gene selection for survival data under dependent censoring, a copula-based
approach, Statistical Methods in Medical Research 25(6): 2840–57.
Emura T, Chen YH (2018). Analysis of Survival Data with Dependent Censoring, Copula-Based Approaches,
JSS Research Series in Statistics, Springer.
Emura T, Ha ID (2021). Special feature: Recent statistical methods for survival analysis. Japanese Journal of35

35Statistics and Data Science, doi:10.1007/s42081-021-00140-0
Emura T, Hsu J H. (2020). Estimation of the Mann–Whitney effect in the two-sample problem under dependent
censoring. Computational Statistics & Data Analysis, 150, 106990.
Emura T, Kao FH, Michimae H (2014). An improved nonparametric estimator of sub-distribution function forbivariate competing risk models, Journal of Multivariate Analysis 132: 229-41.
Emura T, Konno Y (2012a). Multivariate normal distribution approaches for dependently truncated data,
Statistical Papers 53(1): 133-49
Emura T, Konno Y (2012b). A goodness-of-fit tests for parametric models based on dependently truncated data,
Computational Statistics & Data Analysis 56: 2237-50.
Emura T, Lai CC, Sun LH (2021). Change point estimation under a copula-based Markov chain model forbinomial time series, Econometrics and Statistics, doi:10.1016/j.ecosta.2021.07.007.
Emura T, Matsui S, Rondeau V (2019). Survival Analysis with Correlated Endpoints, Joint Frailty-Copula
Models, JSS Research Series in Statistics, Springer
Emura T, Michimae H (2017). A copula-based inference to piecewise exponential models under dependent
censoring, with application to time to metamorphosis of salamander larvae, Environmental and Ecological
Statistics 24(1): 151–73.
Emura T, Nakatochi M, Murotani K, Rondeau V (2017). A joint frailty-copula model between tumour
progression and death for meta-analysis, Statistical Methods in Medical Research 26: 2649-66.
Emura T, Pan CH (2020) Parametric maximum likelihood inference and goodness-of-fit tests for dependently
left-truncated data, a copula-based approach, Statistical Papers 61:479-501.
Emura T, Sofeu C, Rondeau V (2021). Conditional copula models for correlated survival endpoints: individual
patient data meta-analysis of randomized controlled trials, Statistical Methods in Medical Research,doi:10.1177/09622802211046390.
Emura T, Shih JH, Ha ID, Wilke RA (2020). Comparison of the marginal hazard model and the sub-distribution
hazard model for competing risks under an assumed copula, Statistical Methods in Medical Research 29(8):2307–27
江村剛志, 道前洋史 (2020). コピュラを用いた生存時間解析-相関のあるエンドポイントとメタ分析の活用-, 統計数理 第 68 巻(第 1 号): 147-174 (in Japanese).
江村剛志 (2021). コピュラ・マルコフ連鎖モデルによる定常時系列解析-パラメトリック推定と統計的工程管理-, 日本統計学会誌 第 51 巻(第 1 号): 41-73 (in Japanese).
Genest C, MacKay RJ (1986). Copules archimédiennes et families de lois bidimensionnelles dont les marges
sont données. Canadian Journal of Statistics 14(2): 145-159.
Gumbel EJ (1960) Distributions de valeurs extremes en plusieurs dimensions. PubL Inst Stat. Pari 9: 171-3.
Ha ID, Jeong JH, Lee Y. (2017). Statistical Modelling of Survival Data with Random Effects: H-likelihood
Approach, Springer, Singapore.
Ha, ID, Lee Y. (2021), A review of h-likelihood for survival analysis, Japanese Journal of Statistics and Data
Science, doi:10.1007/s42081-021-00125-z.
Hougaard P (1995). Frailty models for survival data. Lifetime Data Analysis, 1(3), 255-273.
Huang JS, Dou X, Kuriki S, Lin GD (2013). Dependence structure of bivariate order statistics with applications

36

36

to Bayramoglu’s distributions, Journal of Multivariate Analysis, 114: 201-208.
Huang X, Emura T (2021). Model diagnostic procedures for copula-based Markov chain models for statisticalprocess control, Communications in Statistics-Simulation and Computation 50(8): 2345-67
Hürlimann W (2017). A comprehensive extension of the FGM copula, Statistical Papers, 58: 373-392.
Joe H (1993). Parametric families of multivariate distributions with given margins. Journal of Multivariate
Analysis 46(2), 262-82.
Joe H (1997). Multivariate Models and Dependence Concepts. Chapman & Hall/CRC, New York.
Kawakami R, Michimae H, Lin YH (2021). Assessing the numerical integration of dynamic prediction formulas
using the exact expressions under the joint frailty-copula model, Japanese Journal of Statistics and Data
Science, doi:10.1007/s42081-021-00133-z.
Kim JM, Baik J, Reller M (2021). Control charts of mean and variance using copula Markov SPC andconditional distribution by copula. Communications in Statistics-Simulation and Computation 50(1): 85-102.
Klein JP, Moeschberger ML (2003). Survival Analysis: Techniques for Censored and Truncated Data, New York:
Springer.
Kotz S, Balakrishnan N, Johnson NL (2000). Continuous Multivariate Distributions, Volume 1: Models, New
York: Springer.
Kwon S, Ha ID, Shih JH, Emura, T. (2021). Flexible parametric copula modeling approaches for clusteredsurvival data. Pharmaceutical Statistics, doi:10.1002/pst.2153.
Leone FC, Nelson LS, Nottingham RB (1961). The folded normal distribution. Technometrics 3(4),543-550.
Li D, Hu XJ, Wang R (2021). Evaluating association between two event times with observations subject toinformative censoring. Journal of the American Statistical Association, doi:10.1080/01621459.2021.1990766.
Lipowski C, Lo SM, Shi S, Ralf A. Wilke RA (2021). Competing risks regression with dependent multiple
spells: Monte Carlo evidence and an application to maternity leave, Japanese Journal of Statistics and Data
Science, doi.org/10.1007/s42081-021-00110-6.
Lodhi C, Mani Tripathi Y, Kumar Rastogi M (2021). Estimating the parameters of a truncated normal
distribution under progressive type II censoring. Communications in Statistics-Simulation and Computation,
50(9): 2757-2781
Lo SM, Stephan G, Wilke RA (2017). Competing risks copula models for unemployment duration: An
application to a German Hartz reform. Journal of Econometric Methods 6(1).
Long TH, Emura T (2014). A control chart using copula-based Markov chain models, Journal of the Chinese
Statistical Association 52(4): 466-96.
Lu JC, Bhattacharyya GK (1990). Some new constructions of bivariate Weibull models. Annals of the Instituteof Statistical Mathematics, 42(3), 543-559.
Marshall AW, Olkin I. (1988). Families of multivariate distributions. Journal of the American Statistical
Association, 83(403), 834-841.
Moradian H, Larocque D, Bellavance F (2019). Survival forests for data with dependent censoring. Statistical
Methods in Medical Research, 28(2), 445-461.
Mohtashami-Borzadaran V, Amini M, Ahmadi J (2019). A generalized bivariate lifetime distribution based onparallel-series structures. Kybernetika, 55(3), 435-454.

37

37

Nelsen RB (2006) An Introduction to Copulas, (2nd ed.). Springer.
Oakes D (1989). Bivariate survival models induced by frailties. Journal of the American Statistical Association
84:487-493.
Ota S, Kimura M (2021). Effective estimation algorithm for parameters of multivariate Farlie–Gumbel–
Morgenstern copula. Japanese Journal of Statistics and Data Science, doi:10.1007/s42081-021-00118-y.
Peng M, Xiang L, Wang S (2018) Semiparametric regression analysis of clustered survival data with semicompeting risks, Computational Statistics & Data Analysis 124: 53-70.
Piancastelli LS, Barreto-Souza W, Mayrink VD (2020). Generalized inverse-Gaussian frailty models withapplication to TARGET neuroblastoma data. Annals of the Institute of Statistical Mathematics
doi.org/10.1007/s10463-020-00774-z
Prenen L, Braekers R, Duchateau L (2018). Investigating the correlation structure of quadrivariate udder
infection times through hierarchical Archimedean copulas, Lifetime Data Analysis. 24:719-742
Rivest LP, Wells MT (2001). A martingale approach to the copula-graphic estimator for the survival function
under dependent censoring. Journal of Multivariate Analysis 79(1): 138-155.
Rondeau V, Pignon JP, Michiels S, MACH-NC Collaborative Group. (2015). A joint model for the dependence
between clustered times to tumour progression and deaths: A meta-analysis of chemotherapy in head andneck cancer. Statistical Methods in Medical Research 24(6): 711-729.
Rotolo F, Legrand C, Van Keilegom I (2013). A simulation procedure based on copulas to generate clusteredmulti-state survival data. Computer Methods and Programs in Biomedicine, 109(3), 305-312.
Šeliga A, Kauers M, Saminger-Platz S, Mesiar R, Kolesárová A, Klement EP (2021). Polynomial bivariate
copulas of degree five: characterization and some particular inequalities. Dependence Modeling, 9(1), 13-42.
Saminger-Platz S, Kolesárová A, Šeliga A, Mesiar R, Klement EP (2021). The impact on the properties of the
EFGM copulas when extending this family. Fuzzy Sets and Systems 415:1-26
Schneider S, Demarqui FN, et al. (2020). An approach to model clustered survival data with dependent
censoring. Biometrical Journal 62(1):157-174.
Shih JH (2014). Copula models, in Handbook of Survival Analysis, Chapter 24, pp.489-510, CRC Press.
Shih JH, Chang YT, Konno Y, Emura T (2019). Estimation of a common mean vector in bivariate meta-analysis
under the FGM copula, Statistics 53(3): 673-95.
Shih JH, Emura T (2018). Likelihood-based inference for bivariate latent failure time models with competingrisks under the generalized FGM copula, Computational Statistics 33(3): 1293-23
Shih JH, Emura T (2019). Bivariate dependence measures and bivariate competing risks models under thegeneralized FGM copula, Statistical Papers 60(4): 1101-18
Shih JH, Emura T (2021). On the copula correlation ratio and its generalization, Journal of Multivariate Analysis
182: 104708
Shinohara S, Lin YH, Michimae H, Emura T (2020) Dynamic lifetime prediction using a Weibull-based
bivariate failure time model: a meta-analysis of individual-patient data, Communications in StatisticsSimulation and Computation, DOI:10.1080/03610918.2020.1855449
Sklar A (1959). Fonctions de répartition à n dimensions et leurs marges", Publ Inst Statist Univ Paris 8: 229-31.
Sofeu C, Emura T, Rondeau V (2021) A joint frailty-copula model for meta-analytic validation of failure time

38

38

surrogate endpoints in clinical trials, Biometrical Journal 63(2): 423-46
Sun LH, Huang XW, Alqawba, MS, Kim JM, Emura T (2020). Copula-based Markov Models for Time SeriesParametric Inference and Process Control, JSS Research Series in Statistics, Springer
Sundberg R (1974). On estimation and testing for the folded normal distribution. Communications in StatisticsTheory and Methods, 3(1): 55-72.
塚原英敦(2021). リスク解析における接合関数, 日本統計学会誌 第 51 巻(第 1 号): 101-121 (inJapanese).
Vaupel JW, Manton KG, Stallard E (1979). The impact of heterogeneity in individual frailty on the dynamics
of mortality. Demography, 16(3), 439-454.
Wang YC, Emura T (2021). Multivariate failure time distributions derived from shared frailty and copulas,
Japanese Journal of Statistics and Data Science, doi:0.1007/s42081-021-00123-1
Wang YC, Emura T, Fan TH, Lo SM, Wilke RA (2020). Likelihood-based inference for a frailty-copula model
based on competing risks failure time data, Quality and Reliability Engineering International 36(5):1622-38
Whitmore GA, Lee MLT. (1991). A multivariate survival distribution generated by an inverse Gaussian mixtureof exponentials. Technometrics 33(1), 39-50.
Wu BH, Michimae H, Emura T (2020). Meta-analysis of individual patient data with semi-competing risks
under the Weibull joint frailty-copula model. Computational Statistics 35:1525-52.
Wu K, Wang L, Yan L, Lio Y (2021). Statistical inference of left truncated and right censored data fromMarshall–Olkin bivariate Rayleigh distribution. Mathematics, 9(21):2703
Zeng X, Gui W (2021). Statistical inference of truncated normal distribution based on the generalized
progressive hybrid censoring. Entropy, 23(2), 186.

39

39

研究所報(最近刊行分)
号数

タイトル

刊行年月日

30 国連ミレニアム開発目標と統計
31 Workshops on “the Population Censuses” and “the Use of
Census Micro Data”
32 ミクロデータとその利用
33 International Symposia on Population Census andMicro Data Archives
34 政府統計の二次的利用

2003.10.20
2003.12.20
2004.04.20
2005.01.10
2005.04.20

35

ジェンダー（男女共同参画）統計

2007.02.20

36

人口センサスの現状と新展開

2007.04.01

37

統計における官学連携2007.04.20

38

ジェンダー（男女共同参画）統計

39

社会生活基本調査とその利用

2010.01.15

40

地方統計の現状と課題

2010.09.15

41

Exploring Potential of Individual Statistical Records

2011.11.05

42

観光統計

2013.02.05

43

国民経済計算関連統計の新たなる展開

2014.01.30

44

タウンページデータによる事業所立地分析

2014.02.15

45

フィンランドのビジネス・レジスター

2015.03.20

Ⅱ

2009.02.10

46 19 世紀ドイツ営業統計史研究

2015.07.20

47

地方統計と統計 GIS

2016.01.25

48

首都圏の人口移動

2017.03.10

49

宿泊業及び飲食業の実証分析

2018.08.01

50

サービス分野の生産物分類

2019.01.31

51

全市区町村産業連関表（平成 23 年表）の推計

2019.10.15

52

商業統計調査2021.01.31

53

産業連関表から供給・使用表へ

2021.03.31

研

究

所

報

No.54

2021 年 11 月 30 日
発行所 法政大学 日本統計研究所〒194-0298 東京都町田市相原町 4342
Tel 042-783-2325,6
Fax 042-783-2332
jsri@adm.hosei.ac.jp
発行人 菅 幹雄

40

BULLETIN
OF
JAPAN STATISTICS RESEARCH INSTITUTE
No.54

November 2021

Statistical Modelling
CONTENTS
Foreword
Pair circulas modeling in higher-order Markov process on the circle
Hiroaki OGATA
On a construction of cylinder distribution
Tomoaki IMOTO
Mixture transition distribution modeling forhigher order circular Markov processes
Takayuki SHIOHAMA
Extended stochastic complexity for sparse estimation
via non-convex regularized regression
Yoichi MIYATA
Multivariate failure time distributions generated from a copula and mixing effects
Takeshi EMURA

Edited byJAPAN STATISTICS RESEARCH INSTITUTE
HOSEI UNIVERSITYTOKYO, JAPAN