離散イベント時系列からの潜在的な競争の推定
Estimation of Latent Competition from Discrete Event Time Series∗1

宮澤脩一 ∗1∗2

持橋大地 ∗3

Shuichi Miyazawa

Daichi Mochihashi

総合研究大学院大学∗2∗3

AGC 株式会社

The Graduate University for Advanced Studies

AGC Inc.

統計数理研究所

The Institute of Statistical Mathematics

To investigate the industrial market, we analyze patent information to estimate the latent competitive relationships between industrial ﬁrms. In existing studies, text analysis methods are often used for patent analysis
to detect correlations between technologies and companies. However, existing methods cannot capture the competitive relationships between companies that might provide deeper insights for technology marketing. To solve
this problem, we propose a model named Ordinary Diﬀerential Equation integrated Log Gaussian Cox Process(ODE-LGCP) that combines the ODE systems, for modeling competition in this study, with generative processesof event data. In the experiment, we use the real patent data and show that our proposed model can estimate
competitive relationships within companies.

1.

はじめにObserved events
Competitor 1

特許は産業界の企業間競争において企業が自社技術を保護する最も基本的な方法であり，特許データの解析は産業界の技術トレンドを捉えたり，重要な技術を保有する企業を発見するのに有用である. 特許情報を技術マーケティングに活用するための解析は多くの文献があり，例えば，[Tang 12] は文書の潜在トピックを推定する確率的トピックモデルにより，産業界の技術トレンドや各技術テーマに関連する企業の分析を行っている．しかし，既存手法は技術-企業の相関関係は推定できるものの，産業界の構造に対するより深い洞察を与えるであろう企業間の競争関係を推定することができない．そのため本稿では，イベントデータの生成過程に常微分方程式 (ODE: OrdinaryDiﬀerential Equation) によってモデル化された競争システムを仮定する ODE 対数ガウス Cox 過程 (ODE-LGCP) を提案する．提案モデルによって，ODE システムの逆問題，すなわち観測データから ODE システムパラメータを推定する問題を近似的に解くことができ，本稿においては，企業の特許出願イ
ベント時系列データから，企業間の競争関係をモデル化する
ODE システムのパラメータを定量的に推定する．
図 1 に提案モデルの概要を示す．ODE システムの逆問題は，通常は離散時系列上で連続的に観測される数値データを扱うが，本研究では，データは連続時系列上で離散的にしか観測されず，なおかつ数値のないイベントデータである．そこで，ガウス過程により補完を行い，観測のない時刻も含めて潜在的な関数を推定し，その関数を支配する ODE を求める．提案モデルは，対数ガウス Cox 過程を用いてイベントの生成過程をモデル化し，イベントの生成されやすさを表す強度関数に ODEに支配された潜在時系列を用いることで，離散観測イベントデータから ODE パラメータを推論することを可能にする．
ODE システムの逆問題の解法については，ODE パラメータを所与として数値積分でデータをシミュレートし，観測データを再現するようなパラメータを探索するのが基本的なアプローチであるが，このアプローチは計算量が多く，大規模問題へ適用できない．計算量の問題への対応策として，勾配マッチングと呼ばれるベイズ的手法が近年成功を収めており，本稿にCompetitor 2Competitor 3

𝑃𝑜𝑖𝑠 Λ 𝑡

Latent space
Latent functions 𝝀 ∈ ℝ"

TimeODE systems
for modeling competition𝑑𝑥!
= 𝑓! (𝒙, 𝜽)𝑑𝑡
𝜆 = 𝜆" 𝑒 #Time

Latent functions 𝒙 ∈ ℝ

Gaussian Processes
for approximating ODE systems

𝒙𝒌 ~ 𝐺𝑃(𝟎, 𝑲𝝓! )Time

図 1: 提案モデルの概要．ODE ダイナミクスを近似する潜在時系列を用いた離散イベントの生成モデルであり，未知の ODEパラメータを離散イベントから推論することを可能とする．
おいても，勾配マッチングの一種である適応的勾配マッチング(AGM: Adaptive Gradient Matching) [Dondelinger 13] を用
い ODE パラメータの推定を行う．
実験では，米国特許庁 (USPTO) の公開する特許データを用いて，産業界の特定の技術分野における企業間競争関係を推定する．特許を分類するため国際的に統一された分類である
国際特許分類 (IPC) を用い，特許母集団を作成する．さらに，発明の要約テキストに対してトピックモデルを用い，母集団内の技術的なトピックを推定後，特定のトピックに関連する出願イベント時系列データに対し提案モデルを適用することで，細分化された技術テーマごとに企業間競争関係を推定する．推定された ODE パラメータによって，競争の観点で注目すべき企業間関係が発見されることを示す．

2.

予備知識

競争をモデル化するための ODE である Lotka-Volterra 競争方程式並びに提案モデルの前提となる適応的勾配マッチング

連絡先: 宮澤脩一，総合研究大学院大学，miyazawa@ism.ac.jp1

(AGM) および対数ガウス Cox 過程 (LGCP) を紹介する．

2.1

Lotka-Volterra 競争方程式

ODE はシステムの振る舞いをモデル化する有用なツールで幅広い学問分野で応用されている．本稿では，競争をモデル化する ODE として，数理生物学において代表的な，生物種間の競争関係を記述する Lotka-Volterra 競争方程式 [Smale 76] を用いる．Lotka-Volterra 競争方程式は，K 種の生物種の競争下における i 種の個体数 pi の時間変化を下式で表す．
(dpi= ri pidt

∑Kj=1 ai,j pj1−)

ηi.

(1)
図 2: ODE 対数ガウス Cox 過程推論モデルのグラフィカルモ
デル．

ここで，ri , ηi は種 i の内的自然増加率，環境収容力と呼ばれる正の定数，ai,j は，種 j の個体数が種 i の個体数増加率に与える影響の大きさを示す競争係数と呼ばれる正の定数である．
特に，ai,i = 1 ∀i = 1, . . . , K であり，ai,j は種 i の内的競争の影響の大きさに対する種 j の影響の大きさの比を示す．
pi は種 i の個体数を表す非負の変数であるが，本稿では xi ∈R とし，pi = exi と書き換えて，ODE を実数変数を用いて定義しなおす．すると，xi についての導関数は，θ = (r, η, a)と ODE パラメータをまとめて以下の通り表せる．

(ẋi = fi (x, θ) = ri

2.2

∑K1−j=1 ai,j eの確率モデルを，二つ目の確率項は ODE の勾配 ẋODE の確率モデルを表しており，勾配を評価する時刻において，回帰曲線が ODE の導関数 f (X, θ) の制約に従った勾配を持つようX, θ が学習される．なお，g k , Ψk は k 個目のガウス過程回帰曲線の勾配の平均および共分散であり，下式により得られる．
−1′′′⊤−1 ′g k = K′⊤
ϕ k K ϕ k xk , Ψk = K ϕ k − K ϕ k K ϕ k K ϕ k .

)
xjηi.

ここで，Kϕk は xk の共分散，K′ϕk は xk およびその導関数の交差共分散，K′′
ϕk は xk の導関数の自己共分散を示す．

(2)

2.3

適応的勾配マッチング

(3)

3.

ODE 対数ガウス Cox 過程

(4)

3.1

生成モデル

(5)

ODE 対数ガウス Cox 過程は，ODE ダイナミクスを近似する潜在時系列関数を用いた対数ガウス Cox 過程，すなわちODE に従った時系列上のイベントデータの生成モデルである．
θ をパラメータとする K 個の式から構成される ODE システ
ムの下で，k 個目の導関数に従った時系列関数を Fk (t)，区間(k)(k)内平均発生率 λ0 として，強度関数を λk (t) = λ0 Fk (t) と定義する Cox 過程から k 種目のイベントデータを生成する．

p(Ẋ, X, θ, ϕ, γ) = p(X|ϕ)p(θ)p(ϕ)p(γ)×K∏p (x˙k |X, θ, ϕ, γ) .

対数ガウス Cox 過程

時系列上のイベントの発生時刻について，イベントの発生率が時刻によって変動する場合の確率過程として Cox 過程がよく知られている．時系列上の区間 A におけるイベント発生数∫は，時刻 t における発生率を λ(t) として，Λ(A) = A λ(t)dtをパラメータとするポアソン分布 P o(Λ(A)) に従うと考える．
ここで，時間によって値が変動する未知の強度関数 λ(t) をデータから推定するために関数生成器としてガウス過程を仮定した
確率過程を対数ガウス Cox 過程 [Møller 98] という．非負の強度関数を実数空間上で定義される関数 f (t) によって表すため，指数変換を用いて λ(t) = λ0 · ef (t) (≥ 0) とする．ここで，λ0は区間内の平均的な発生率とする．

勾配マッチングは，ODE システムの逆問題のベイズ的アプローチによる近似解法であり，計算量の多い数値積分を行うことなく ODE パラメータ θ を推定する．時系列観測データ Y に対し，潜在的な回帰曲線 X をフィットすることで ODE ダイナミクスを近似し，回帰曲線の勾配および ODE の導関数の値が近い値を取るような確率モデルを定義する．[Dondelinger 13]
は，回帰にガウス過程を用いる適応的勾配マッチング (AGM:Adaptive Gradient Matching) を提案した．
AGM では，K 個の式から ODE システムが構成されるとすると，モデルの同時分布を回帰曲線の勾配 Ẋ を周辺化したものとして以下のように定義する．

p(Y , X, θ, ϕ, σ, γ) = p(Y |X, σ)p(σ)p(X, θ, ϕ, γ),∫p(X, θ, ϕ, γ) = p(Ẋ, X, θ, ϕ, γ)dẊ,(7)k=1

ここで，ϕ はガウス過程のカーネル関数 k(t, t′ ) のパラメータ，
γ は ODE 勾配の確率モデルのノイズの分散を表す．さらに，(5) 式に含まれる勾配の条件付き確率分布を，ガウス過程回帰曲線の勾配の確率モデルと ODE に従う勾配の確率モデルをProduct-of-Experts[Hinton 02] の枠組みで統合した確率モデルとして以下のように定義する．

3.2

推論モデル

図 2 に提案手法の推論モデルに関するグラフィカルモデルを示す．ϕk をパラメータとするガウス過程から生成する潜在時系列 xk ∼ N (0, Kϕk ) を用いて Fk (t) ≈ exk と近似することを考える．そのために，(4)– (6) 式の AGM の枠組みを用い，xk の回帰曲線の勾配 ẋk が，(2) 式のように実数変数の導関数に変換した ODE の勾配 fk (x, θ) とマッチするような尤度を導入し，潜在時系列 exk が ODE に従うように誘導する．
通常の AGM は観測時系列データ y に対し，ガウス過程で回帰関数をフィッティングするが，提案モデルでは観測時系列ではp (x˙k |X, θ, ϕ, γ) ∝ p (x˙k |xk , ϕ) p (x˙k |X, θ, γ)= N (x˙k |g k , Ψk ) N (x˙k |f k (X, θ), γI) .
(6)
ここで，右辺一つ目の確率項はガウス過程回帰曲線の勾配 ẋGP

2

なく潜在時系列 y に対するフィッティングを行う．xk , y k は，そW
れぞれ有限の時点数 T, W の時刻 τ = {τt }T
t=1 , ω = {ωw }w=1 ,においてのみ評価するのものとし，潜在時系列 y k は潜在時系
列 xk のガウス過程の予測分布に標準偏差 σk の出力ノイズが加えて以下のように生成する．

表 1: トピックの例とそれを代表する確率上位の単語.
トピック

12345y k ∼ N (µ̂k , Σ̂k ),−1µ̂k = K⊤
xy Kxx xk

上位単語

content item ﬁle page web page displayword text sentence candidate model term
message query entity natural language node question
user interface template provide ﬁeld form
text display image unit input character

(8)

−12
Σ̂k = Kyy − K⊤xy Kxx Kxy + σk I

ここで，Kxx , Kxy , Kyy は，それぞれカーネル関数値 k(τt , τt′ ),k(τt , ωw ), k(ωw , ωw′ ) を要素とする行列である．こうして生成された実数空間上の潜在時系列 y k を指数変換すると，Fk (t) ≈eyk と近似でき，生成モデルで用いる対数ガウス Cox 過程の強度関数が得られる．
対数ガウス Cox 過程の学習のため，連続時間を時間間隔 ∆の W 個の時間窓に離散化し，それぞれの窓でポアソン分布の尤度を評価する．w 番目の時間窓の中間時刻 ωw における k 種の潜在状態を yk,w ，同窓内における k 種のイベント頻度をmk,w∑(k)とし，各時間窓における平均イベント頻度 λ0w mk,w=∆

図 3: トピック 2 に関連する特許出願イベントの時系列．
分類の粒度が粗い国際特許分類を用いて集めた特許集合に対し，より細分化された技術的なトピック毎に分析を行うため，対象特許集合の要約テキストに潜在的に含まれるトピックのモデリングを行った．LDA[Blei 03] に代表されるトピックモデルは，文書をトピックの多項分布とし，さらにトピックを単語の多項分布とする階層モデルである．特に，Correlated Topic
Models (CTM)[Blei 06] は文書のトピック事前分布を多変量正規分布とするモデルであり，グローバルなトピックの出現確率やトピック間の相関関係を推定できる．本稿では，対数尤度を指標 [Griﬃths 04] として用いトピック数を 9 と設定してCTM によるトピックの推論を行った．表 1 に一部のトピックのトピック–単語分布における確率上位 6 単語を示した．これ
らの各トピックについて，各文書の文書–トピック分布における当該トピックの確率が経験的に設定した閾値 15%を超える特許のみを各トピックの関連特許として抽出した．本稿では，
図 3 で示されるトピック 2 に関する特許出願件数上位 5 企業
A, B, C, D, E のイベントデータのみに実験対象を絞り．5 社間の競争モデルを仮定した実験結果を示す．

を(k)用いて Pois(mk,w ; λ0 eyk,w ) の尤度を評価する．
まとめると，提案モデルの同時分布は下式の通り表される．

p(m, y, x, σ, θ, ϕ, γ) = p(m|y)p(y|x, σ, ϕ)p(σ)∫× p(ẋ, x, θ, ϕ, γ)dẋ,p(m|y) =∏∏w

kp(y|x, σ, ϕ) =∏(k)Pois(mk,w ; λ0 eyk,w ),N (y k ; µ̂k , Σ̂k ),k∫p(ẋ, x, θ, ϕ, γ)dẋ ∝ p(θ)p(ϕ)p(γ)∏N (xk ; 0, Kxx )k1⊤× e− 2 (f k −gk ) (Ψk +γI)

−1

(f k −g k )

.

4.2

(9)

3.3

MCMC 法による推論(9) 式は解析的に計算できないため，MCMC による近似推論を行う．ODE パラメータ θ に対してはロジット変換，その他のパラメータ σ, ϕ, γ に対しては対数変換を行って得られる実数変数について，それぞれ事前分布を正規分布とすることで，ハミルトニアンモンテカルロ法 (HMC)[Neal 11] 等の実数空間上の潜在状態に対する推論アルゴリズムが利用可能になる．また，[Dondelinger 13] で指摘されている通り，AGMの推論において，非線形な ODE システムの影響で，事後分布はごつごつとした多峰形状となる可能性があるため，レプリカ交換モンテカルロ法 [Hukushima 96] によりサンプラーが局所的な溝に陥ることを防ぐ．

4.

実験

4.1

実験データの準備

実験設定

本稿では，y, x, θ のみを推論対象としてサンプリングを行い，σ, γ, ϕ については推論対象とはせずハイパーパラメータとして扱った．ハイパーパラメータは経験的に σ = 0.02，γ = 0.1と設定した．また，ガウス過程のカーネル関数にRBF カーネ()′ 2|
ル k(x, x′ ) = ϕ21 exp − |x−x2ϕ2を用い，ϕ1 = 5.0, ϕ2 = 0.22で固定した．また，ODE パラメータの値域は，自然増加率
r について [1, 10]，環境収容力 η について [5, 15]，競争係数a について [0, 10] を用いた．潜在変数 x および y を評価する時刻については，T = 21, W = 100 で， 0 から 1 に正規化した時刻において，τ = [0.00, 0.05, 0.10, . . . , 0.95, 1.00]，
ω = [0.005, 0.015, 0.025, . . . , 0.985, 0.995] とした．
MCMC シ ミュレ ー ション は ，異 な る 10 段 階 の 温 度
1.10 , 1.11 , . . . , 1.19 でチェーンのレプリカを作成し，各温度で HMC により事後分布からのサンプルを得，チェーン間のスワップ試行を毎イテレーション行う設定とした．HMC のパラメータは，リープフロッグステップ数は L = 30 で固定し，ステップサイズは事前のシミュレーションで調整した値を用いた．2,000 回のバーンインシミュレーションを行い，その後の10,000 回分のシミュレーションで得られたサンプルを 10 サン
プル毎に間引いて収集した 1,000 個のサンプルを用いて事後パラメータ推論を行った．

提案手法の有効性を検証するために，米国特許庁 (USPTO)が公開する現実の特許データに対して，特定の技術テーマに関する企業間競争関係を推定する．Handling natural language
data という国際特許分類 (IPC) グループを含む 2005 年から2020 年の間に出願され特許の内，第一出願人に組織名が含まれる特許に限定し，計 17,422 の特許を解析対象として選定した．

3

A社

B社

C社

D社

E社

210−1−2
0.00

0.25

0.50

0.75

1.00

0.00

0.25

0.50

0.75

1.00

0.00

0.25

0.50

0.75

1.00

0.00

0.25

0.50

0.75

1.00

0.00

0.25

0.50

0.75

1.00

図 4: Lotka-Volterra 競争方程式に従うダイナミクスを近似する各社の潜在関数の事後サンプル. X 軸: 時刻，Y 軸: xk の値.
2.0
1.5
1.0
0.50-0.5-1.0-1.5-2.0A
B
CDEA

B

CD(a)

5.

E

本稿では．時系列イベントデータの生成過程に任意の ODEシステムを仮定することを可能とする ODE 対数ガウスを提案した．産業界の企業間の競争関係をモデル化するために，ODEとして Lotka-Volterra 競争方程式を用い，特定の技術に関連する特許出願イベント時系列から ODE パラメータの推定を行った．推定された ODE パラメータの競争係数により，企業間の 1 対 1 の競争関係における影響の度合いを定量化し，注目すべき関係性を発見し，解釈・考察を与えることを可能にした．
本稿では，技術分野別にモデリングしたが，企業活動が複数の異なる技術分野をまたがることも多いため，ある企業内における異なる技術分野間の資源の奪い合いもまた競争モデルで表現できる．企業間の競争と企業内の競争を同時にモデリングすることで，より示唆に富む推定結果が得られると期待される．
また，提案手法は特許テキストの内容に関しては考慮せず，単にイベントデータとして扱うに限っている．本稿では特許母集団を仕分けるために前段でテキスト解析を行い，その結果を提案モデルによるモデリングに用いているが，確率的トピックモデルや文書分散表現といったテキストデータに対する手法とODE システムを直接組み合わせたモデルも考えられる．

DAC
BE(b)

図 5: 競争関係の推定結果. (a) i 行 j 列の要素が 対数競争係数
log aij の事後平均を示す隣接行列．赤が影響の強い関係，青が影響の弱い関係を示す．(b) 競争関係の有効グラフ．エッジの太さは対数競争係数の絶対値に対応し，影響の度合いは，赤く太いエッジ，赤く細いエッジ，内的競争，青く細いエッジ，青く太いエッジの順に強い．

4.3

まとめと今後の展望

結果と考察

図 4, 5 にトピック 2 に関連する特許出願データからの提案手法によるパラメータの推定結果を示した．図 4 は，各社の潜在関数 xk の事後分布からのサンプルであり，関数の浮き沈みは 3 のイベントの出現傾向を反映している．また，図 5 に競争関係の推定結果を可視化した．(a) は，i 行 j 列の要素を対数競争係数 log aij とする隣接行列である．対角成分は各企業の内的な対数競争係数 log aii = 0 であり，これを基準として白色で示し，log aij > 0 を赤，log aij < 0 を青で示している．(b) は (a) の隣接行列をネットワーク状に可視化した結果である．競争係数の値は，関数 x に反映されており，log aijが大きい場合，i の関数の下降と j の関数の上昇が同時期に起こり，log aij が小さい場合，j の関数の i の関数への影響が弱いことを示す．
推定された対数競争係数からは，A, C, D が E から強い影響を受けることが示唆されており，A, C, D の潜在関数の下降と E の潜在関数の上昇が同時期に起きている結果に対応する．
上記の結果は，競争係数が推定できなくとも，例えば対数ガウス Cox 過程で関数 x の形状さえ可視化できれば，目視で発見可能かもしれない．しかしながら，B のように上昇，下降，再上昇とより複雑な動きをする関数については目視で企業間の関係性を捉えることは難しい．また，競争モデルで取り扱う企業数 K を増やした場合，企業間の 1 対 1 の関係の数は O(K 2 )で増加するため，目視での発見は困難になるため，競争係数行列を図 5(a) のように可視化できることは実用上意味がある．
log aij < 0 である青で示される関係性に着目すると，企業j の企業 i への影響の小ささもまた評価できる．C, D の関数は似た形状で時系列としては相関するため，一見すると関係性があるようにも見える．しかしながら，競争係数からは互いに影響しあわない関係であるという推定結果が得られており，これは関数の目視確認では発見することが難しい関係であろう．

参考文献[Blei 03] Blei, D. M., Ng, A. Y., and Jordan, M. I.: Latent dirichlet
allocation, JMLR, Vol. 3, No. Jan, pp. 993–1022 (2003)[Blei 06] Blei, D. and Laﬀerty, J.: Correlated topic models, Advances in Neural Information Processing Systems, Vol. 18, p.
147 (2006)[Dondelinger 13] Dondelinger, F., Husmeier, D., Rogers, S., andFilippone, M.: ODE parameter inference using adaptive gradient matching with Gaussian processes, in AISTATS, pp. 216–228PMLR (2013)[Griﬃths 04] Griﬃths, T. L. and Steyvers, M.: Finding scientiﬁctopics, Proceedings of the National Academy of Sciences, Vol.
101, No. suppl 1, pp. 5228–5235 (2004)[Hinton 02] Hinton, G. E.: Training products of experts by minimizing contrastive divergence, Neural Computation, Vol. 14, No. 8,pp. 1771–1800 (2002)[Hukushima 96] Hukushima, K. and Nemoto, K.: Exchange Monte
Carlo Method and Application to Spin Glass Simulations, Journalof the Physical Society of Japan, Vol. 65, No. 6, p. 1604–1608(1996)[Møller 98] Møller, J., Syversveen, A. R., and Waagepetersen, R. P.:
Log Gaussian Cox Processes, Scandinavian Journal of Statistics,Vol. 25, No. 3, pp. 451–482 (1998)[Neal 11] Neal, R. M., et al.: MCMC using Hamiltonian dynamics,
Handbook of Markov Chain Monte Carlo, Vol. 2, No. 11, p. 2
(2011)[Smale 76] Smale, S.: On the diﬀerential equations of species in competition, Journal of Mathematical Biology, Vol. 3, pp. 5–7 (1976)[Tang 12] Tang, J., Wang, B., Yang, Y., Hu, P., Zhao, Y., Yan, X.,
Gao, B., Huang, M., Xu, P., Li, W., et al.: Patentminer: topicdriven patent analysis and mining, in Proceedings of the 18th
ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pp. 1366–1374 (2012)4