Expression for the skew normal distribution specifying up to third-order moments independently

—Analysis of skewness of brand value—

要　旨

本研究の目的は，3 次までの積率を明示的にかつ独立に特定できる確率分布を構成することである．非対称正規分布の母数に変換を施して，平均・分散・歪度を直接パラメタライズできるように新たな確率密度関数を構成する．3 次までの積率を独立に特定することで，(1) 統計モデルの一部として組み込んだ際に，直接歪度を推定すること，(2) 潜在変数（因子）の歪度を直接推定すること，(3) 群間で歪度を比較することが可能となる．非対称正規分布とχ² 分布に関するシミュレーション研究により，母数の推定における妥当性が確認された．また，歪度が観察されやすいブランド価値データに提案手法を適用した結果，集団間での3 次までの積率の違いを細かく比較できることが示された．母数推定には，マルコフ連鎖モンテカルロ法によるベイズ推定を用い，サンプリング手法にはハミルトニアンモンテカルロ法を利用した．

1. 目的

標準正規分布

に従う変数x を，2 つの定数μ₁ とμ₂ を用いて，

とアフィン変換すると，μ₁ は原点回りの1 次の積率（母平均）となり，μ₂ は平均回りの2 次の積率（母分散）となることは，よく知られている．アフィン変換の定数が，そのまま母数になる確率分布は非常に珍しく，正規分布

が，統計解析にもっとも利用される確率分布である主要な理由の1つは，2 次までの積率を明示的に（陽に），かつ，独立に母数そのものとして特定でき，この性質がたいへん便利なためである．

分布の位置（1 次の積率，平均）と分布の散らばり（平均回りの2 次の積率，分散）に加え，統計解析において分布の歪みを考察することは重要である．分布の歪みを表現した確率分布の1 つに，非対称正規分布（skew normal distribution）がある．しかし，従来の非対称正規分布の表現では，潜在変数（因子）の歪度を直接推定したり，群間で歪度の違いを検討したりすることは難しい．本論文の目的は，母歪度μ₃（標準得点の3 次の積率）に着目し

のように，非対称正規分布において，3 次までの積率を明示的に（陽に），かつ，独立に母数そのものとして表現できる確率分布を構成することである．3 次までの積率を明示的に表現することの利点として，(1) 統計モデルの一部として組み込んだ際に，直接歪度を推定することができる，(2) 潜在変数（因子）の歪度を直接推定できる，(3) 群間で歪度の比較を行うことができる，などが挙げられる．以下の章では，3 次までの積率を明示的に特定するための方法について説明し，次いで，シミュレーション研究の結果を示す．そして，適用例において，ブランド価値に関する分析を行い，提案手法によって，潜在変数の歪度の推定や，群間（2 国間）での潜在変数の平均・分散・歪度の比較が可能となることを示す．

2. 方法

確率変数X が，密度関数

に従うとき，X は母数c の非対称正規分布と呼ばれる(Azzalini, 1985)．ここで

である．非対称正規分布の積率母関数は

であり，これを利用すると

が導かれる(Azzalini, 1985)．ただし

である．

2つの定数a, b を用いて，非対称正規分布に従う確率変数X を

とアフィン変換すると，確率分布は，変換のヤコビアンより

となる(Azzalini & Capitano, 1999; Capitano, Azzalini, & Stanghellini, 2003; 蓑谷, 2010)．

ここで独立変数μ₃, μ₂, μ₁ を導入し，以下の3 つの関数g₃, g₂, g₁

を構成する．ただし

である（d は実数の解とする）．

そして確率変数はそのままにして

のように母数のみを変換した新たな確率密度関数に注目する．(7) 式と(8) 式を(5) 式に代入して整理すると，

となる．つまり，関数g₂, g₁は，確率変数X を非対称正規分布の母平均E[X] と母分散V [X] を用いて標準化した後に，再びアフィン変換を施しているのと同値だから，

である．またg₃ は，(4) 式の左辺をμ₃ とおいてc に関して解いた関数である．μ₃ とc は1 対1対応しており，アフィン変換によって母歪度は変化しないから，

である．ただし母歪度μ₃ には制限がある．非対称正規分布のpdf はc→±∞のとき半正規分布のpdf に収束するから，半正規分布の歪度以上の歪みは表現できない．具体的には(4) 式の極限をとると

であることが分かる(Azzalini & Valle, 1996)．したがって約±1.0 を超える歪度は扱うことができない．

以上のように，扱える歪みの程度に制約はあるものの，本表現を利用すれば，3 次までの積率を明示的に（陽に）かつ独立に母数そのものとして特定できる．具体的には，任意のμ₁, μ₂, μ₃ を(6), (7), (8) 式を利用してa, b, c に変換し，p(y|a, b, c) に代入することで，平均μ₁，分散μ₂，歪度μ₃ の非対称正規分布を構成・表現することができる．

3. 母数の推定

母数の推定には，マルコフ連鎖モンテカルロ法（Markov Chain Monte Carlo; MCMC 法）によるベイズ推定を利用する．MCMC 法の中でも，ハミルトニアンモンテカルロ法（HamiltonianMonte Carlo method, Duane, Kennedy, Pendleton, & Roweth, 1987; HMC 法）というサンプリング手法を利用する．HMC 法は，ハミルトニアン力学系による確定的な状態遷移と，メトロポリス法による確率的な状態遷移を融合した方法であり(松本, 2004)，母数の対数事後分布の1次微分が容易に求められる場合に適用可能である．HMC 法は，ハイブリッド法とも呼ばれ，仮想的な運動量を導入して，位置変数と運動量変数の結合空間である位相空間上でのエルゴード的なマルコフ連鎖を考える．

ギブスサンプリングやMH アルゴリズムなどの従来のMCMC 法と比べ，HMC 法のメリットとして(1) 共役事前分布を設定しなくても採択率が高い，(2) 条件付き分布によって1 つずつ更新するのではなくすべての母数を同時に更新できる，(3) 時系列相関が低いといった点が挙げられる．以下では，事前分布と事後分布の導出について説明し，ハミルトニアンモンテカルロ法については，付録で説明を行う．

3.1. 母数の推定

ベイズ推定では，母数に関する事前情報を事前分布として表現する．提案した非対称正規分布p(y|μ₁, μ₂, μ₃) の3 つの母数μ₁, μ₂, μ₃ に関して，以下の事前分布を設定する．

平均μ₁ には平均α₁，分散β₁ の正規分布（N）を，分散μ₂ の逆数には形状母数α₂，尺度母数β₂のガンマ分布（G）を，歪度μ₃ には区間[α₃, β₃] の一様分布（U）を仮定する¹．

3.2. 事後分布

確率変数y_i（i = 1, · · · ,N）が独立に非対称正規分布p(y|μ₁, μ₂, μ₃) に従うとすると，尤度関数L(y|μ₁, μ₂, μ₃) と事前分布p(μ₁, μ₂, μ₃) より，母数の事後分布p(μ₁, μ₂, μ₃|y) は

となる．ただし，母数の事前分布は互いに独立であると仮定する．また，母数の対数事後分布は以下となる．

4. シミュレーション

提案手法において，適切に母数の推定が行えるかを検証するために，シミュレーション研究を行った．ここでは，2 つの分布から発生させたデータに対して提案手法を適用する．1 つは，推定の妥当性を確かめるために非対称正規分布から発生させ，もう1 つは，発生モデルが異なる場合の推定値の振る舞いを確かめるためにχ² 分布から発生させた．サンプルサイズN は100, 500,1000, 5000 の4 通りとする．また，非対称正規分布については，歪度μ₃ を0.3, 0.6, 0.9 の3 通りとし，χ² 分布については，自由度df を89, 22, 10 の3 通りとした．χ² 分布の歪度は，となるため，df =89, 22, 10 のχ² 分布の歪度はそれぞれ0.299, 0.603, 0.894 となる．

サンプルサイズと歪度の組合せ（12 通り）に関して，それぞれで50 回分析を行い，母数μ₃ の推定値²の平均，平均平方誤差（Root Mean Squared Error; RMSE），母数μ₃ の事後分布の標準偏差（SD）の平均を次式で算出した．

μ_3,true は真値を，添え字k はシミュレーションの繰り返しを，σ_μ3k は各回の母数μ₃ の事後分布のSD を表す．母数の推定には，NUTS 法を実装しているソフトウェアStan Version 1.3.0 (StanDevelopment Team, 2013) と統計ソフトR を利用した³．5000 回サンプリングを行い，バーンイン期間として最初の1000 個の標本を破棄して，残りの4000 個の標本で事後統計量を算出した．また，シミュレーションデータの発生には，R の乱数発生関数を利用した．

上述の分析に加えて，同じシミュレーションデータに対して，標本歪度γ を計算し，提案手法との比較を行った．標本歪度γ とその標準誤差se_γ は，Joanes and Gill (1998) を参考に次式で算出した．

そして，各条件において，標本歪度γ の平均γ_Mean とRMSE_γRMSE を計算した．

4.1. 非対称正規分布に関するシミュレーション

2 つの確率変数Z₁ とZ₂ が独立で，それぞれ標準正規分布に従うとき，

は(1) 式で示した母数c の非対称正規分布に従う(蓑谷, 2012)．そこで，(6) 式でμ₃ をc に変換し，標準正規乱数から(14) 式を利用して非対称正規分布に従うN 個の乱数を発生させた．

提案手法によって得られた推定値 に関する結果（μ_3Mean, μ_3RMSE, σ_μ3Mean）を表1 に，また，標本歪度γ に関する結果（γ_Mean, γ_RMSE, se_γ）を表2 にまとめる．なお，標本歪度の標準誤差se_γ は歪度の値に関係なく，サンプルサイズによって決まるため，表2 には各サンプルサイズに対応した標準誤差を示す．

表1 非対称正規分布に関するシミュレーション結果

表2 非対称正規分布に関するシミュレーション結果（標本歪度）

提案手法によって得られた推定値 に関して，各回で得られた推定値の平均は概ね真値と近い値を示している．また，RMSE は，歪度0.3, 0.6, 0.9 の3 条件の中では，歪度0.9 の場合に小さく，またサンプルサイズが増えると小さくなる傾向にある．事後分布のSD の平均もRMSE と同様の傾向を示した．

推定値 と標本歪度γ を比べると，推定値の平均および標本歪度の平均にはそれほど大きな違いはないが，RMSE は提案手法における推定値 の方が小さいことが分かる．また，母数μ₃の事後分布のSD の平均σ_μ3Mean と標本歪度の標準誤差se_γ の比較においても，提案手法の方が値は小さい．特に，サンプルサイズが小さいときに両者の差は大きく，サンプルサイズが増えるにつれて，μ_3RMSE とγ_RMSE の差，および，σ_μ3Mean とse_γ の差は小さくなることが確認できる．

図1，図2，図3 には，ある回のシミュレーションデータのヒストグラムと，真値および推定値 から求められる確率密度をプロットした図を示す．実線が真値の確率密度，ダッシュ線が推定値の確率密度を表している．図1 が歪度0.3，図2 が歪度0.6，図3 が歪度0.9 のプロットであり，左からサンプルサイズ100, 500, 1000, 5000 である．サンプルサイズが少ない場合は，実線とダッシュ線は一致しないが，多くなるとほぼ一致する傾向にあることが分かる．

図1 非対称正規分布のシミュレーションデータのヒストグラムと密度関数（μ₃ = 0.3）

図2 非対称正規分布のシミュレーションデータのヒストグラムと密度関数（μ₃ = 0.6）k-means 法

図3 非対称正規分布のシミュレーションデータのヒストグラムと密度関数（μ₃ = 0.9）

4.2. χ² 分布に関するシミュレーション

続いて，χ² 分布から乱数を発生させ，シミュレーションデータを作成した．先ほどと同様に，提案手法による推定値 に関する結果（μ_3Mean, μ_3RMSE, σ_μ3Mean）を表3 に，また，標本歪度γに関する結果（γ_Mean, γ_RMSE, se_γ）を表4 にまとめる．

表3 χ² 分布に関するシミュレーション結果

表4 χ² 分布に関するシミュレーション結果（標本歪度）

提案手法による推定値 に関して，各回の推定値の平均は，歪度が小さい場合には，サンプルサイズが増えると真値に近づくが，歪度が大きい場合には真値と多少のずれが生じていることが分かる．RMSE は，サンプルサイズが増えるほど小さくなる傾向にあるが，非対称正規分の場合と異なり，歪度が小さい方がRMSE は小さいことが確認できる．一方，事後分布のSD の平均は，歪度が大きいほど小さく，サンプルサイズが多いほど小さい傾向にある．

推定値 と標本歪度γ の各回の平均を比べると，標本歪度γ の方が真値に近い値を示している．また，標本歪度γ のRMSE は，推定値 のRMSE と同様に，歪度が小さい方が値が小さくなることが確認できる．

先程と同様に，ある回のシミュレーションデータのヒストグラム，χ² 分布（実線），推定値から求められる非対称正規分布（ダッシュ線）を図4，図5，図6 に示す．サンプルサイズが少ない場合はズレが大きいが，多くなるにつれて比較的似た形状を示していることが分かる．χ² 分布の場合には，非対称正規分布の場合と比べて少しズレが生じているが，甚だしく分布を表現できていないという訳ではなく，χ² 分布の場合でも，全体的な傾向は再現できている．

2 つの分布を利用したシミュレーション研究より，限られた中ではあるが，提案モデルにおいて適切に母数推定を行うことができると考えられる．ただし，歪んだ分布はχ² 分布の他にも多数あり，総合的な頑健性については今後検証が必要である．

5. 適用例

データに無視できない歪度が観察されるのは，大別して2 つの場合があろう．1 つ目は異質な発生過程を有するデータが含まれて歪む場合である．例としては，100 点満点の試験のデータに入力ミスによって200 点が観察されたデータや，一般人と少数の関取が混じった体重のデータを挙げることができる．これらの場合は，データを削除したり，集団別に分析したり，混合分布をあてはめたりなど，状況に応じて対処することが多い．

図4 χ² 分布のシミュレーションデータのヒストグラムと密度関数（μ₃ = 0.300, df = 89）

図5 χ² 分布のシミュレーションデータのヒストグラムと密度関数（μ₃ = 0.603, df = 22）

図6 χ² 分布のシミュレーションデータのヒストグラムと密度関数（μ₃ = 0.894, df = 10）

2つ目は変数や観測対象自身の性質として歪む場合である．例としては，収入や貯蓄高の分布が有名である．歪度は，格差の1 つの指標として参照することができる．収入や貯蓄高は直接観測される変数である．他方，ブランド価値や企業価値の分布などのある種の構成概念も，少数の有名ブランド（ガリバー的存在のブランド）と多数の発展途上のブランドの価値の差が大きく，測定値が正に歪む場合が多い．その歪度は，ブランド間格差の指標として解釈でき，要約統計量として実質科学的意味を有する．そこで，ブランド価値に関するデータに提案手法を適用し，歪度に関する解釈を行う．

5.1. ブランドアジア

ブランドアジア2013 (日経BP コンサルティング, 2013b) は，アジア12 地域における主要グローバルブランドおよびローカルブランドのブランド力を調査・測定する目的で実施された．ブランドアジア2013 では，ブランドイメージの観点からブランド力を測定するブランドジャパン(日経BP コンサルティング, 2013a) の評価・分析方法を踏襲し，消費者に対して消費行動上のブランド価値を示すコンシューマ指標（以下BtoC とする）が計算される．具体的には，2 次因子分析を行い，「総合力」と4 つの下位指標である「フレンドリー」「イノベーディブ」「アウトスタンディング」「コンビニエント」を算出する．調査における観測変数およびパス図を表5 および図7に示す⁴．

表5 観測変数

図7 ブランドアジア2013 のパス図

ここでは，ブランドアジア2013 のデータの一部を利用して，ブランド価値の歪みの推定を行う．アジア12 地域のうち，インドネシアとマレーシアを取り上げ，2 か国の違いを考察する．インドネシアとマレーシアの調査手法および回収数を表6 に示す．また，調査対象者は20 歳以上の男女であり，調査期間は2012 年11 月から2013 年1 月である．

ブランドアジアでは，全地域に共通する60 のグローバルブランドと各地域で選定された40 のローカルブランドを合わせて，100 ブランドを調査対象としている⁵．共通する60 ブランドおよび各地域固有のブランドを表7 に示す．インドネシアおよびマレーシアのそれぞれのブランドが表5 に示した15 個の観測変数によって評価される．

表6 調査概要

表7 グローバルブランドとローカルブランド

5.2. モデル

ここでは，平均0，分散1 に標準化されたブランドi (i = 1, · · · ,N) の測定z_i = (z_i1, · · · , z_ij , · · · , z_i15)′ を，2 次因子分析モデルを利用して

と表現する．(15) 式のΛ は，サイズ（15 × 4）の1 次因子パタン行列であり，

である．また，f_i は，1 次因子ベクトルであり，次式で表される．

4 つの1 次因子は，表5 に示したフレンドリー，イノベーティブ，アウトスタンディング，コンビニエントの4 つの下位指標である．λ は，2 次因子パタンベクトルであり，

である．g_i が総合的なブランド価値を表す2 次因子である． ε_i はサイズ15，ζi_i はサイズ4 の誤差変数ベクトルであり，

のように正規分布に従う．Σ_ε とΣ_ζ は対角行列であり，観測変数の分散が1.0 であることから，

という制約が置かれる．

5.3. 事前分布と事後分布

母数θ をθ = (Λ,λ, f_i, g_i, (Σ_ε), (Σ_ζ), μ₁, μ₂, μ₃) と置く⁶と，全データZ に関する尤度L(Z|θ)は

と表現される．また，母数θ の事前分布は，互いに独立であることを仮定し

と表される．よって，母数θ の事後分布は，

となる．また，対数事後分布は，以下となる．

観測変数が標準化されていることから，1 次因子パタンλ_jk の事前分布は，区間[−0.999, +0.999]の一様分布とする．区間を[−1.0, +1.0] としなかったのは，誤差の分散が著しく小さくなるのを防ぐためである．また，1 次因子の分散を1.0 に固定するために，同様の理由から，2 次因子パタンλ_k の事前分布は，区間[−0.999, +0.999] の一様分布とする．

総合的なブランド価値を表す2 次因子g_i の事前分布は，非対称正規分布の表現を利用し，

とする⁷．さらに，μ₃ の事前分布は，区間[−0.99, +0.99] の一様分布とする．

5.4. 分析例1

自社ブランドの価値を創造・強化する際には，既存ブランドの総合力（2 次因子）g_i の分布を確認することが重要となる．かりに，ガリバー的存在のブランドが無ければ，企業努力によって，自社ブランドを上位ブランドとして確立できる可能性がある．一方，ガリバー的存在のブランドがあるならば，総合力ではなく，特定の分野（例えば，フレンドリーやコンビニエントなど）に特化してブランド価値を強化することも考えられる．そこで，インドネシアとマレーシアのデータに対して，上述の総合力g_i に3 次までの積率を特定できる分布を仮定したモデルを適用し，歪度を推定することによって，既存ブランドの総合力の分布について考察を行う．シミュレーション研究と同様に，NUTS 法を実装しているStan とR を利用して，HMC 法で推定を行った．バーンイン期間を1000 回とし，その後4000 個の標本を利用して事後統計量（事後分布の平均値，標準偏差（SD），中央値）を算出した．分析結果を表8 に，また，母数μ₃ のトレースと分布を図8 に示す．

図8 μ₃ のトレースと分布（インドネシア（上）とマレーシア（下））

表8 推定結果

すべての母数に関してGeweke 指標の絶対値が，1.96 以内に収まっており，マルコフ連鎖の収束が示唆される．また，トレース図では，似た値が出続ける傾向が無く，状態空間を行き来ているため，収束していると考えられる．推定された歪度μ₃ から，ガリバー的存在のブランドの総合力g_i がどれほど突出しているかを検討することができる．実際に，両国の歪度μ₃ をみると，インドネシアとマレーシアの歪度は共に正の値であり，ガリバー的存在のブランドの総合力が突出して高いと解釈することができる．

生活の中で携帯電話や通信環境をよく利用するインドネシアでは，IT・家電メーカーやインターネット企業（Nokia, Google, Apple, Telmomsel など）のブランド価値が高く，世界的に有名で，品質が良く信頼性の高い憧れのブランドが上位にランクインしている．一方，マレーシアでは，国内最大手の銀行（Maybank）と自動車メーカー（HONDA, TOYOTA, Mercedes-Benz）が上位を占める．公共交通機関が発達していないマレーシアでは自家用車への依存率が高いため，自動車メーカーのブランド価値が高く評価されている．また，両国ともに，食品や衣類などの日常生活で利用するブランドが，下位に位置する結果となった．

インドネシアとマレーシア共に，総合力g_i の分布は正に歪んでいるが，文化や生活の違いにより，上位を占めるガリバー的存在のブランドは異なる．インドネシアでは，IT・インターネット系のブランド価値が高いために分布が正に歪むのに対して，マレーシアでは，銀行と自動車メーカーが高く評価されるために正に歪むことが示唆された．また，両国の歪度μ₃ を比べると，インドネシアの方が値は大きいが，それほど大きな差ではないことが分かる．歪度μ₃ に差がある場合には，ガリバー的存在のブランドの総合力の突出度合いが異なるため，2 国間でブランド価値の創造・強化に対して異なるアプローチを考える必要がある．しかし，単純に両国間で総合力（2 次因子）の分布の歪度μ₃ を比較することはできず，群間で潜在変数の分布の歪みを検討するためには，測定不変を仮定した上で多母集団解析を行う必要がある．分析例2 では，2 国間で総合力g_iの分布の歪度に違いが見られるかを検討する．

5.5. 分析例2

インドネシアとマレーシアの総合力（2 次因子）g_i の分布の歪みの違いについて検討するために，提案した2 次因子分析を多母集団に拡張して，分析を行った．分析では，3 次の積率だけでなく，平均や分散の違いも検討するため，複数のモデルを仮定した⁸．仮定したモデルを表9 に示す．総合力の平均・分散・歪度に関して様々な仮定を置き，モデル比較を行う．表中の∗ は自由母数を表し，数値は，固定母数を表している．M2 は歪度だけが異なるモデルであり，M3 は歪度が等しいという等値制約を置いたモデルである．また，M4 は分散と歪度が異なるモデルであり，M5 は分散が異なり歪度が等しいという等値制約を課したモデルである．なお，これらのモデルは，測定不変が成り立つという仮定の下で分析を行った．分析例1 と同様に5000 回のサンプリングを行って，後半の4000 個の標本を利用して事後統計量を算出した．事後分布の平均値，標準偏差（SD），中央値を表10 にまとめる⁹．

表9 モデル比較

分析の結果，DIC の観点からは，M3 の適合が良いことが分かった．よって，インドネシアとマレーシアの総合力（2 次因子）の分布の歪度には違いがないことが示唆された．両国ともに正に歪んでおり，上位のランキングには違いがみられるが，全体的な分布の特徴として，ブランドの格差にはそれほど大きな違いは無いと考えられる．つまり，ガリバー的存在のブランドの総合力の突出度合いは同程度であるため，両国においては，類似した方法でブランド価値の創造・強化を行うことが企業戦略の1 つとして考えることができる．以上のように，提案手法を利用することで，群間で平均・分散・歪度に差が見られるかどうかを，詳細に検討できることが示された．

表10 多母集団解析の推定結果

6 考察

本研究の目的は，3 次までの積率を明示的にかつ独立に特定できる確率分布を構成することであった．シミュレーション研究より，提案手法において適切に歪度を推定できることが確認できた．ただし，シミュレーションでは2 つの分布のみに注目して推定精度の評価を行っており，その他の様々な歪んだ分布を含めた総合的な頑健性については今後検討が必要である．一方，適用例では，3 次までの積率を独立に特定することで，潜在変数の歪度を直接推定できること，また，2つの国の潜在変数の分布の違いをより詳しく比較できることが示された．適用できる範囲は歪度が約±1 と限定されるものの，提案手法は歪んだ分布の分析を行う際に有益であると考えられる．

本論文で提案した確率分布の表現方法に関する課題を2 点述べる．1 つ目は，扱うことのできる歪度の範囲である．本表現では±1 を超える歪度を表現することはできない．この問題に対処するには，対称的なpdf であるg(x) とそのcdf であるG(x) を用いて

のように非対称分布を構成する方法(Azzalini & Capitano, 2003) を用いることが有望であろう．裾の重たいg(x) を利用すれば歪度の範囲が広がる．

2 つ目は，多変量分布への表現の拡張である．非対称正規分布は，Azzalini and Valle (1996),Arnold and Beaver (2004), Valle (2004) などによって，自然な形式で多変量分布に拡張されている．その多変量非対称正規分布に，本表現法を適用することは可能である．しかし多変量同時分布の歪みを実質科学的に納得できるように解釈することは難しく，現時点では実践的開発意義に乏しいように思われる．

脚 注

1  正規分布の母数に関する推測では，共役事前分布として，平均には正規分布を，分散の逆数にはガンマ分布を設定することが多い．

2  推定値 には，母数μ₃ の事後分布の平均を用いた．

3  事前分布として，平均μ₁ には正規分布N(0, 10⁴) を，分散μ₂ の逆数にはガンマ分布G(10⁻³, 10⁻³) を，歪度μ₃ には一様分布U(0, 0.99) を設定した．なお，NUTS 法については，付録を参照のこと．

4  ブランドジャパン(日経BP コンサルティング, 2013a) の分析においても，同様のパス図および観測変数を利用している．

5  ただし，シンガポールはグローバルブランドのみ，マレーシアは20 のローカルブランドを利用した．

6  本研究では，誤差分散に制約を置くため，Σ_ε とΣ_ζ は母数推定の対象とはならない．

7  潜在変数の単位は任意に定めることができるため，ここでは，平均μ₁ と分散μ₂ をそれぞれ0 と1 に固定する．

8  自由母数として推定する場合には，事前分布として，平均μ₁ には標準正規分布N(0, 1.0) を，分散の逆数1/μ₂にはガンマ分布G(10⁻³, 10⁻³) を仮定する．

9  表10 の“—” は，固定母数として扱っているために値が算出されない部分を表す．

10  ベイズ推定の枠組みでは，位置変数が母数に相当する．

11  ハミルトン力学には，時間反転性（reversibility），ハミルトンH の不変（Conservation of the Hamiltonian），位相空間の体積保存（Volume Preservation）という重要な3 つの性質がある（Neal (2009) を参照のこと）．

12  U(θ) = −log p(θ|Data) と置いてもよいが，事後分布の分母は定数であり規格化定数としてZ_H に吸収される．

13  詳しくは，Bishop (2012) を参照のこと．

付録

A. ハミルトニアンモンテカルロ法

A.1. ハミルトン方程式

力学の観点から，物体の位置変数θ = (θ₁, · · · , θ_n) は時間t により変化すると考え，各位置変数に対して対応する運動量変数p = (p₁, · · · , p_n) を導入する¹⁰．物体の運動に関して，位置エネルギーU(θ) は位置変数θ によって，運動エネルギーK(p) は運動量変数p によって与えられる．エネルギー保存則より，系の全エネルギーは位置エネルギーと運動エネルギーの和である．

(A.1) 式のH(θ, p) はハミルトン関数と呼ばれる．ハミルトン関数を利用することで，ハミルトンの運動方程式¹¹が

で与えられ，この方程式を解くことで，物体の運動の変化を知ることができる．なお，運動エネルギーは一般に，以下のように定義される．

A.2. 数値積分

ハミルトン方程式を解くためには，数値積分を利用する必要があり，HMC 法ではleap frog 法が利用される．leap frog 法は，ステップサイズεと更新回数L を定め，はじめに運動量変数をステップサイズε/2 で半ステップ更新し，次にステップサイズεで位置変数の更新を行い，そして運動量変数の2 回目の半ステップ更新を行う．leap flog 法は，体積保存の性質を保ったまま数値積分を行うことができる手法であり，更新式は以下となる．

A.3. ハミルトニアン力学を利用したMCMC

位置変数θ と運動量変数p の位相空間上での同時分布を，ハミルトニアン関数を用いて

と定義する．ここで，位置エネルギーU(θ) を

とする．p(θ) は母数の事前分布であり，L(Data|θ) は尤度を表す¹²．同時分布は位置変数θ と運動量変数p の周辺分布の積で表現できる（独立である）ため，同時分布p(θ, p) からサンプリングを行うことができれば，得られた標本は事後分布p(θ|Data) からのサンプルと見なすことができる．提案手法では，U(θ) に(13) 式の定数を除いた部分を指定することになる．

A.4. アルゴリズム

HMC 法のアルゴリズムを以下に示す．

step1 θ の初期値θ⁰，ステップサイズε，更新回数L を決める．

step2 t = 1, 2, · · · , T に対して，次を繰り返す．

(a) 運動量変数p⁰ を正規分布から発生させる(p⁰ ∼ N(0, I_n))．

(b) θ^t−1 とp⁰ を初期値とし，leap frog 法(ε,L) により，候補点 を求める．

(c) 採択確率α を計算する．

(d) u ∼ U(0, 1) を発生させ，θ^t を以下に従って決める．

ハミルトンH の値は常に一定であるため，ハミルトン力学の下での時間発展によって，p(θ, p)からエルゴード的なマルコフ連鎖を構成することはできない．そこで，エルゴート的なマルコフ連鎖を得るために，分布p(θ, p) を不変に保ちつつ，H の値を変えるような位相空間での追加的な移動を導入する．これを行うのが，p の値をθ を条件とする分布から抽出したもので置き換えることであり，step2 の(a) の手順p⁰ ∼ N(0, I_n) に対応する．ただし，θ とp は独立なため，条件付き分布p(p|θ) からのサンプリングは容易である¹³．また，step2 の(c) で計算される採択確率αは詳細釣り合い条件が満たされるために用いられる．HMC 法における詳細釣り合い条件については，Neal (2009) やBishop (2012) などが詳しい．

A.5. No-U-Turn Sampler (NUTS)

HMC 法を実行する際には，分析者がステップサイズεと更新回数L を定める必要がある．更新回数が短すぎるとランダムウォーク的な振る舞いが現れ，逆に長すぎると無駄に計算時間がかかるため，これらのパラメタを決定するのは難しい．この点を克服した方法が，ノーユーターンサンプラー（No-U-Turn Sampler; NUTS 法）であり，自動でステップサイズεと更新回数L を決定することができる．具体的には，母数空間に関して，現在の位置 と候補点˜θ との距離Q を

と定め，これを更新の停止基準に用いる．更新回数L を増やしても距離Q が大きくならないのであれば，これ以上更新を行うのをやめ，これまでに計算された候補点の中からサンプリングを行う．また，ステップサイズεは，dual averaging 法(Nesterov, 2009) を利用して決定する．

B. 各国の調査票

インドネシアとマレーシアで利用した調査票の質問項目表現を表B.1 に示す．

B.1. 観測変数（質問項目の表現）

References

• Arnold, B. C., & Beaver, R. J. (2004).Elliptical models subject to hidden truncation or selective sampling.In Genton, M. J. (Eds.),Skew-Elliptical Distributions and their Applications (pp.101-112). Chapman & Hall/CRC.
• Azzalini, A. (1985).A class of distributions which includes the normal ones.Scandinavian Journal of Statistics, 12, 171-178.
• Azzalini, A., & Capitanio, A. (1999).Statistical Applications of the Multivariate Skew Normal Distribution.Journal of the Royal Statistical Society. Series B, 61, 579-602.
• Azzalini, A., & Valle, A. D. (1996).The multivariate skew-normal distribution.Biometrika, 83, 715-726.
• Azzalini, A., & Capitano, A. (2003). Distribution, generated by perturbation of symmetry with emphasis on a multivariate skew-t distribution.Journal of the Royal Statistical Society B, 65, 367-389.
• Bishop, C. M. (2012).パターン認識と機械学習 下 —ベイズ理論による統計的予測— (元田浩・栗田多喜夫・樋口知之・松本裕治・村田 昇,監訳).丸善出版. (原本は2011年刊行)
• Capitanio, A., Azzalini, A., & Stanghellini, E. (2003).Graphical Models for Skew-Normal Variates.Scandinavian Journal of Statistics, 30, 129-144.
• Duane, S., Kennedy, A. D., Pendleton, B. J., & Roweth, D. (1987).Hybrid Monte Carlo.Physics Letters B, 195, 216-222.
• Hoffman, M., & Gelman, A. (2012). The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo.Arxiv preprint arXiv:1111.4246
• Joanes, D. N., & Gill, C. A. (1998).Comparing measures of sample skewness and kurtosis.The Statistician, 47, 183-189.
• 松本隆 (2004).非線形ダイナミカルシステムの再構成と予測.石黒真木夫・松本隆・乾敏郎・田辺国士・伊庭幸人・甘利俊一・竹内啓・竹村彰通 (編),階層ベイズモデルとその周辺 —時系列・画像・認知への応用— (pp. 89-170). 岩波書店.
• 蓑谷千凰彦 (2010).統計分布ハンドブック(増補版).朝倉書店.
• 蓑谷千凰彦 (2012).正規分布ハンドブック.朝倉書店.
• Neal, R. M. (2010).MCMC using Hamiltonian dynamics.In Brooks, S., Gelman, A., Jones, G. L., & Meng,X. L. (Eds.),Handbook of Markov Chain Monte Carlo (113-162). Chapman & Hall/CRC.
• Nesterov, Y. (2009).Primal-dual subgradient methods for convex problems.Mathematical programming, 120, 221-259.
• 日経BPコンサルティング (2013a).ブランドジャパン2013総合報告書[解説書].日経BPコンサルティング.
• 日経BPコンサルティング (2013b).BRAND ASIA 2013 DATA BOOK.日経BPコンサルティング.
• Stan Development Team. (2013).Stan Modeling Language User's Guide and Reference Manual, Version 1.3.
• Valle, A. D. (2004).The skew-normal distribution.In Genton, M. J. (Eds.),Skew-Elliptical Distributions and their Applications (pp.1-24). Chapman & Hall/CRC.