axjack's blog

### axjack is said to be an abbreviation for An eXistent JApanese Cool Klutz ###

2019年8月をワールドワイドに考える

メモ プログラミング 調べ物

はじめに

UTCを考慮すると途端にわけわからなくなったので整理する。

UTCとは

wikiで協定世界時 - Wikipediaをみませう。

UTCってプラマイいくつまであるの？

時間帯 (標準時) - WikipediaによるとからUTC-12からUTC+14まである。日本はUTC+9

2019年8月をUTCで列挙する

以下の時間帯で考える。

すると次の表のようになる

UTC-12 UTC0 UTC+9 UTC+14
7/30 22:00 7/31 10:00 7/31 19:00 8/01 00:00
7/31 03:00 7/31 15:00 8/01 00:00 8/01 05:00
7/31 03:00 8/01 00:00 8/01 09:00 8/01 14:00
8/01 00:00 8/01 12:00 8/01 21:00 8/02 02:00
… … … …
8/30 22:00 8/31 10:00 8/31 19:00 9/01 00:00
8/31 03:00 8/31 15:00 9/01 00:00 9/01 05:00
8/31 03:00 9/01 00:00 9/01 09:00 9/01 14:00
9/01 00:00 9/01 12:00 9/01 21:00 9/02 02:00

結局2019年8月とは

  • UTC+14が世界で最も早く8/1を迎えるので、UTC0で考えれば7/31 10:00(⇔日本時間 7/31 19:00)
  • UTC-12が世界で最も遅く9/1を迎えるので、UTC0で考えれば9/1 12:00(⇔日本時間 9/1 21:00)

の時間帯を抽出すればとりあえず取りこぼすことはなさそうだ。

参考

タイムゾーンを考慮した日時の扱いのベストプラクティス - エムスリーテックブログ

多変量解析法入門pp.52-53のテコ比の式変形

統計 勉強

やること

多変量解析法入門p.52の式(4.35)あたりの式変形の行間を埋めてみる。

式変形

式(4.35)にて
\begin{eqnarray} \hat{y_k} &=& \hat{\beta_0} + \hat{\beta_1}x_k \\ &=& \bar{y} + \hat{\beta_1}(x_k - \bar{x} ) \\ &=& \frac{\Sigma y_i}{n} + \frac{ (x_k - \bar{x}) \Sigma(x_i-\bar{x})y_i }{ S_{xx} } \end{eqnarray}


ここで、上式の第一項は

\displaystyle \frac{\Sigma y_i}{n} = \frac{1}{n}y_1 + \frac{1}{n}y_2 + \cdots + \frac{1}{n}y_k + \cdots + \frac{1}{n}y_n \\ \displaystyle = \begin{bmatrix} y_1 & y_2 & \cdots & y_k & \cdots & y_n \end{bmatrix} \begin{bmatrix} \frac{1}{n} & \frac{1}{n} & \cdots & \frac{1}{n} \end{bmatrix}'


第二項は

\displaystyle \frac{ (x_k - \bar{x}) \Sigma(x_i-\bar{x})y_i }{ S_{xx} } = \Bigl( \frac{(x_k - \bar{x} )}{S_{xx}} \Bigr) \Bigl( (x_1-\bar{x})y_1 + (x_2-\bar{x})y_2 + \cdots + (x_n-\bar{x})y_n \Bigr) \\ \displaystyle = \begin{bmatrix} y_1 & y_2 & \cdots & y_k & \cdots & y_n \end{bmatrix} \begin{bmatrix} \Bigl( \frac{(x_k-\bar{x})(x_1 - \bar{x})}{S_{xx}} \Bigr) & \Bigl( \frac{(x_k-\bar{x})(x_2 - \bar{x})}{S_{xx}} \Bigr) & \cdots & \Bigl( \frac{(x_k-\bar{x})(x_k - \bar{x})}{S_{xx}} \Bigr) & \cdots & \Bigl( \frac{(x_k-\bar{x})(x_n - \bar{x})}{S_{xx}} \Bigr) \end{bmatrix}'


より、バラバラにした項を足して元に戻せば、

\displaystyle \hat{y_k} = \begin{bmatrix} y_1 & \cdots & y_k & \cdots & y_n \end{bmatrix} \begin{bmatrix} \Bigl( \frac{1}{n} + \frac{(x_k-\bar{x})(x_1 - \bar{x})}{S_{xx}} \Bigr) & \cdots & \Bigl( \frac{1}{n} + \frac{(x_k-\bar{x})(x_k - \bar{x})}{S_{xx}} \Bigr) & \cdots & \Bigl( \frac{1}{n} + \frac{(x_k-\bar{x})(x_n - \bar{x})}{S_{xx}} \Bigr) \end{bmatrix}' \\ \displaystyle \iff \hat{y_k} = \begin{bmatrix} y_1 & \cdots & y_k &\cdots& y_n\end{bmatrix} \begin{bmatrix} h_{k1} & \cdots & h_{kk} & \cdots & h_{kn} \end{bmatrix}'


となる。

ということで、y_kの係数

\displaystyle h_{kk} = \frac{1}{n} + \frac{ (x_k - \bar{x})^2 }{ S_{xx} }

をテコ比と呼ぶ。式(4.36)

とりあえずアソシエーション分析する〜その３〜

調べ物 統計

手を動かして形は理解したっぽいので、今度こそ(？)仕組みを理解する編です。いきなり読むと分からなくなるので、何となくアソシエーション分析のパッケージ{arules}は動かせたけどイマイチ意味がわからないんだよなーという状態になった時に読むのをお勧めします。


前回↓
axjack.hatenablog.jp

もくじ

アソシエーション分析の内部で使っているデータはどんな形式なのか？

アソシエーション分析の内部で使っているデータとは、{arules}のread.transactions関数で読み込んだものを指しています。この読み込んだデータをSQLで喩えてみるべく、TransというSQLのテーブルにアソシエーション分析の内部で使っているデータが入っているとします。関係スキーマ表記とcreate table 表記を用いてTransテーブルを表すと、以下のようになります。なおここで、アイテムが全部でp種類あるとしましょう。

関係スキーマ表記

Trans( tranID, item_1, item_2, ..., item_p )
となります。下線部は主キーです。

create table 表記

create table Trans(
    tranID nvarchar(max) not null primary key
  , item_1  bit not null
  , item_2 bit not null
  , ...
  , item_p bit not null
);

となります。

記号Mと演算子σを定義する

ここで唐突ですが、記号Mと演算子σをSQLを用いて定義します。
なお、記号「:=」は「左辺を右辺で定義する」、記号「=」は「左辺と右辺は等しい・上の式と右の式は等しい」と解釈してください。

記号M は Transの全行数

記号M は Transの全行数を表します。select文で書くと、

M := select count(*) from Trans;

となります。

演算子σ は Transの抽出条件つきの行数

幾つか例示します。ここでXやYはTransテーブルにあるtranIDではない任意の列(すなわちitem_1〜item_pのどれか)とします。

カッコの中が１つの場合

σ( X ) := select count(*) from Trans where X = 1;

カッコの中が２つの場合

σ( X∪Y ) := select count(*) from Trans where X = 1 and Y = 1;

σ( Y∪X ) := select count(*) from Trans where Y = 1 and X = 1;
= select count(*) from Trans where X = 1 and Y = 1;
= σ( X∪Y )

Xが複合アイテム( X = {X1, X2} )のような場合は、

σ( X ) = σ( {X1, X2} ) := select count(*) from Trans where X1 = 1 and X2 = 1;

となります。*1

相関ルールの３指標は何者なのか？

相関ルールの３指標とは、

  • Support(サポート; 支持度)
  • Confidence(コンフィデンス; 確信度)
  • Lift(リフト; リフト値)

を指します。

さて以下にて、σはSupport・Confidence・ Liftを用いて表せることを示します。

Support(サポート; 支持度)

Suppと略します。

Supp( X ) := σ(X) / M
より、σ(X) = Supp( X ) * M

Supp( X ⇒ Y ) := σ( X∪Y ) / M
より、σ( X∪Y ) = Supp( X ⇒ Y ) * M

ということで、SuppにMを掛けると「ルールが含まれるトランザクション数」を取り出せます。

Confidence(コンフィデンス; 確信度)

Confと略します。

Conf( X ⇒ Y ) := Supp( X ⇒ Y ) / Supp( X ) 
より、Supp( X ) = Supp( X ⇒ Y )/Conf( X ⇒ Y )
ところでSupp( X ) = σ( X ) / Mなので
σ( X ) = ( Supp( X ⇒ Y )/Conf( X ⇒ Y ) ) * M

ということで、SuppをConfで割ってからMを掛けると
「ルールの前件(X ⇒ YにおけるXの部分)が含まれるトランザクション数」が取り出せます。

Lift(リフト; リフト値)

Lift( X ⇒ Y ) := Conf( X ⇒ Y ) / Supp(Y) 
より、Supp( Y ) = Conf( X ⇒ Y )/Lift( X ⇒ Y )
ところでSupp( Y ) = σ( Y ) / Mなので
σ( Y ) = ( Conf( X ⇒ Y )/Lift( X ⇒ Y ) ) * M

ということで、ConfをLiftで割ってからMを掛けると
「ルールの後件(X ⇒ YにおけるYの部分)が含まれるトランザクション数」が取り出せます。


以上をまとめると、

σ Supp,Conf,Lift で表すと
σ( X∪Y ) Supp( X ⇒ Y ) * M
σ( X ) ( Supp( X ⇒ Y )/Conf( X ⇒ Y ) ) * M
σ( Y ) ( Conf( X ⇒ Y )/Lift( X ⇒ Y ) ) * M

となります。

これの何が嬉しいかというと、アソシエーション分析で得られたSupp・Conf・Liftの値を割り算掛け算することによって比率から個数に変換することができるわけです。(なぜならば「演算子σ は Transの抽出条件つきの行数」であるから、です。)

実際のデータで確かめてみる

アソシエーション分析(1)から表を２つ引用します。

表：トランザクション
f:id:axjack:20190806211933p:plain

表：アソシエーション分析の結果
f:id:axjack:20190806211949p:plain

いくつかのルールにてσを求めてみます。ここでM = 5とします。

例１

ルール lhs ⇒ rhs supp conf lift
6 ソーセージ ⇒ オムツ 0.2 1 1.66

は

σ( ソーセージ∪オムツ ) = supp * M = 0.2 * 5 = 1
σ( ソーセージ ) = (supp/conf) * M = (0.2/1) * 5 = 1
σ( オムツ ) = (conf/lift) * M = (1/1.666) * 5 = 3

となります。

例２

ルール lhs ⇒ rhs supp conf lift
11 タバコ ⇒ 弁当 0.2 1 2.5

は

σ( タバコ∪弁当 ) = supp * M = 0.2 * 5 = 1
σ( ソーセージ ) = (supp/conf) * M = (0.2/1) * 5 = 1
σ( オムツ ) = (conf/lift) * M = (1/2.5) * 5 = 2

となります。

例３

ルール lhs ⇒ rhs supp conf lift
17 {オムツ, ソーセージ } ⇒ ビール 0.2 1 1.25

は

σ( {オムツ, ソーセージ}∪{ビール} ) = supp * M = 0.2 * 5 = 1
σ( {オムツ, ソーセージ} ) = (supp/conf) * M = (0.2/1) * 5 = 1
σ( {ビール} ) = (conf/lift) * M = (1/1.25) * 5 = 4

となります。

おわりに

アソシエーション分析{arules}をいじって一週間ぐらい経過しました。支持度・確信度・リフト値の意味がわからずムムムとなっていた状態からはだいぶ進歩したような気持ちです。比率を個数に直すと結構「なるほど！」と思えてきました。

「アソシエーション分析してみたけど、比率だけだと実感が湧かない！」という人は比率を個数に直して考えると理解が深まるかもしれません。また、個数に直してみると「比率はどうなってるの？」とやっぱり比率が恋しく？なりますね笑　

こんなかんじで比率・個数を行き来していくとアソシエーション分析の面白さがより一層分かってくるのかなぁ、と感じました。

参考

en.wikipedia.org

www.cis.doshisha.ac.jp

*1:Association rule learning - Wikipediaによると、"Note that supp(X∪Y) means the support of the union of the items in X and Y. This is somewhat confusing since we normally think in terms of probabilities of events and not sets of items. " と書かれている。一見「∪」の記号からORを連想してしまうが決してそうではないということを注として書いておく。X∪Yに於けるXは集合でありYも集合である。すなわちX = {みかん,ビール}, Y = {ビール, おむつ} ならば X∪Y = {みかん, ビール, おむつ}である。(※ X∩Y は {ビール} なことに注意。) そして結果的にσ( X∪Y )をselect count(*)文に翻訳するとそれぞれのアイテムをbit型な列とした時に「=1」でandを使って結んだ抽出条件と同じになる。この場合はσ( X∪Y ) = select count(*) from Trans where みかん = 1 and ビール = 1 and おむつ = 1である。集合の要素と列で同じ名前を使っていることが混乱の原因なのかもしれないものの、そこは上手くLet's 脳内変換。

とりあえずアソシエーション分析する〜その２〜

前回に引き続きとりアソしてみます。

前回↓ axjack.hatenablog.jp

tranID で group-by した時グループに属するitemが1個のもの、を除外する

ファイルの準備

dat_singleをコピってitemが1つしかないtranIDの行を作成する。

axjack:dat $ diff dat_single.csv dat_single.v2.csv
35a36,41
> 11,ビール
> 12,ビール
> 13,ビール
> 14,ワイン
> 15,ワイン
> 16,ワイン

前処理

アソシエーション分析用のファイルはCSVから読み込むと良い、とのことなのでファイル読み込み→加工→ファイル書き出し→再度読み込みします。

# load library ####
library(dplyr)
library(readr)

# single形式のCSVファイルを読み込む
mytemp <- read_csv(file="dat/dat_single.v2.csv"
                ,col_types = c("cc")
                ,col_names = c("tranID","item")
                ,quote = '\"'
                #,skip = 1 #ヘッダがある場合は1行目をスキップする
                )

# transactionの長さ1ではないものを抽出し
mytemp %>% unique() %>% group_by(tranID) %>% filter( n() != 1 ) -> mytemp.f

# ファイルに書き出す
write_csv(mytemp.f,path="dat/dat_single.f.csv", col_names = FALSE)

末尾の６行が消えたか確認すると、

> tail(mytemp,   n=6)
# A tibble: 6 x 2
  tranID item  
  <chr>  <chr> 
1 11     ビール
2 12     ビール
3 13     ビール
4 14     ワイン
5 15     ワイン
6 16     ワイン
> tail(mytemp.f, n=6)
# A tibble: 6 x 2
# Groups:   tranID [3]
  tranID item      
  <chr>  <chr>     
1 8      チーズ    
2 9      ワイン    
3 9      牛肉      
4 9      じゃがいも
5 10     ワイン    
6 10     こんにゃく

消えました。ではCSVファイルをread.transactionsで読み込みます。

データの確認

###########################
# アソシエーション分析 ####
###########################
# load library #### 
library(arules)

# load CSV ####
# (1)元データ
ds <- read.transactions(
  file = "dat/dat_single.csv"
  , format = "single"
  , sep = ","
  , cols = c(1,2)
  , rm.duplicates = T
#  , skip = 1
)

# (2)レコード追加したデータ
dsv2 <- read.transactions(
  file = "dat/dat_single.v2.csv"
  , format = "single"
  , sep = ","
  , cols = c(1,2)
  , rm.duplicates = T
  #  , skip = 1
)

# (3) "(2)"からlength=1のtransactionを削除したデータ
dsf <- read.transactions(
  file = "dat/dat_single.f.csv"
  , format = "single"
  , sep = ","
  , cols = c(1,2)
  , rm.duplicates = T
  #  , skip = 1
)

トランザクションの行数列数を確認すると

> ds
transactions in sparse format with
 10 transactions (rows) and
 19 items (columns)
> dsv2
transactions in sparse format with
 16 transactions (rows) and
 19 items (columns)
> dsf
transactions in sparse format with
 10 transactions (rows) and
 19 items (columns)

たしかにdsv2は6トランザクション消すことに成功しています。summaryをみて見ると、

> summary(ds)
transactions as itemMatrix in sparse format with
 10 rows (elements/itemsets/transactions) and
 19 columns (items) and a density of 0.1842105 

most frequent items:
ソーセージ     ビール     ワイン       牛肉     生野菜    (Other) 
         3          3          3          3          3         20 

element (itemset/transaction) length distribution:
sizes
2 3 4 5 
1 4 4 1 

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    2.0     3.0     3.5     3.5     4.0     5.0 
> summary(dsv2)
transactions as itemMatrix in sparse format with
 16 rows (elements/itemsets/transactions) and
 19 columns (items) and a density of 0.1348684 

most frequent items:
    ビール     ワイン ソーセージ       牛肉     生野菜    (Other) 
         6          6          3          3          3         20 

element (itemset/transaction) length distribution:
sizes
1 2 3 4 5 
6 1 4 4 1 

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.000   1.000   3.000   2.562   4.000   5.000

よい感じです。あとはaprioriに作ったtransactionを食わせればアソシエーション分析の完成かと。ここから先は前回と同じなので省略。

参考

とりあえずアソシエーション分析する

調べ物 統計

タイトル通りです。とりあえず動かしてみたい！という願望でやってみました。理解は後からついてくると信じて、まずは手を動かしたいという人向け？の記事です。

読み込むデータの準備

前処理なしのトランザクションデータを{arules}パッケージで読み込む方法 - 六本木で働くデータサイエンティストのブログ にあるデータに１行追加したものを利用します。gistにアップロードしました。

ところでどうやら、CSVの最終行に空行が無いと読み込みエラーになるっぽい？ので、エラーが出たら空行付け足してみてください。

dat_basket.csvはbasket形式のデータで、dat_single.csvは同データをsingle形式にしたものです。

gist.github.com

データの読み込み

read.transactionsの引数はどうなってるんじゃ？という人は?read.transactionsしてみましょう。冷静に読むと丁寧に説明が書かれております。

# load library #### 
library(arules)

# load CSV ####
# format = "basket" の場合
mytran_basket <- read.transactions(
  file = "dat/dat_basket.csv"
  , format = "basket"
  , sep = ","
  , cols = NULL
  , rm.duplicates = T
)

# format = "single" の場合
mytran_single <- read.transactions(
  file = "dat/dat_single.csv"
  , format = "single"
  , sep = ","
  , cols = c(1,2)
  , rm.duplicates = T
)

生成したtranデータの確認

mytran_basketもmytran_singleも同じ内容なので、mytran_basketを使います。

> mytran_basket
transactions in sparse format with
 10 transactions (rows) and
 19 items (columns)
> summary(mytran_basket)
transactions as itemMatrix in sparse format with
 10 rows (elements/itemsets/transactions) and
 19 columns (items) and a density of 0.1842105 

most frequent items:
    ビール ソーセージ     柿ピー       牛肉     生野菜    (Other) 
         4          3          3          3          3         19 

element (itemset/transaction) length distribution:
sizes
2 3 4 5 
1 4 4 1 

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    2.0     3.0     3.5     3.5     4.0     5.0 

includes extended item information - examples:
      labels
1     アイス
2     おむつ
3 こんにゃく

適当にピックアップすると、

transactions in sparse format with
 10 transactions (rows) and
 19 items (columns)

transactionsは10件あり、品目数は19件と云うことがわかります。dat_basket.csvは10行だしdat_single.csvのtransaction-IDは10個なので確かにあってそうです。

most frequent items:
    ビール ソーセージ     柿ピー       牛肉     生野菜    (Other) 
         4          3          3          3          3         19

transactionに含まれる品目のうち、最も多いのはビールで次はソーセージで・・・と云うことがわかります。CSVの中を検索すると確かにあっています。

element (itemset/transaction) length distribution:
sizes
2 3 4 5 
1 4 4 1

transactionの長さとはtransactionに含まれる品目の数だそうで、長さ2と長さ5が1件で長さ3と長さ4が4件だと云うことがわかります。目視で確認すると確かにあっています。

というわけで、ここまでは読み込んだtransactionについての情報でした。

いよいよアソシエーション分析

流れは、apriori関数にデータを突っ込んでrulesを生成→inspectで中身を見る、です。

> mytran_basket.apr <- apriori(mytran_basket, parameter = NULL)
Apriori

Parameter specification:
 confidence minval smax arem  aval originalSupport maxtime support minlen maxlen target   ext
        0.8    0.1    1 none FALSE            TRUE       5     0.1      1     10  rules FALSE

Algorithmic control:
 filter tree heap memopt load sort verbose
    0.1 TRUE TRUE  FALSE TRUE    2    TRUE

Absolute minimum support count: 1 

set item appearances ...[0 item(s)] done [0.00s].
set transactions ...[19 item(s), 10 transaction(s)] done [0.00s].
sorting and recoding items ... [19 item(s)] done [0.00s].
creating transaction tree ... done [0.00s].
checking subsets of size 1 2 3 4 5 done [0.00s].
writing ... [136 rule(s)] done [0.00s].
creating S4 object  ... done [0.00s].

[136 rule(s)]ということで136の相関ルールが抽出されました。inspectしてみましょう。

> inspect( head(sort(mytran_basket.apr,by="support"),15) )
     lhs                    rhs            support confidence lift      count
[1]  {柿ピー}            => {ビール}       0.3     1           2.500000 3    
[2]  {じゃがいも}        => {牛肉}         0.2     1           3.333333 2    
[3]  {果物}              => {豆腐}         0.2     1           3.333333 2    
[4]  {果物}              => {生野菜}       0.2     1           3.333333 2    
[5]  {ソーセージ,柿ピー} => {ビール}       0.2     1           2.500000 2    
[6]  {ソーセージ,ビール} => {柿ピー}       0.2     1           3.333333 2    
[7]  {果物,豆腐}         => {生野菜}       0.2     1           3.333333 2    
[8]  {果物,生野菜}       => {豆腐}         0.2     1           3.333333 2    
[9]  {生野菜,豆腐}       => {果物}         0.2     1           5.000000 2    
[10] {チーズ}            => {ワイン}       0.1     1           5.000000 1    
[11] {チーズ}            => {ソーセージ}   0.1     1           3.333333 1    
[12] {おむつ}            => {ドレッシング} 0.1     1          10.000000 1    
[13] {ドレッシング}      => {おむつ}       0.1     1          10.000000 1    
[14] {おむつ}            => {ビール}       0.1     1           2.500000 1    
[15] {おむつ}            => {生野菜}       0.1     1           3.333333 1    
>

とりあえず一行目のこれ↓

     lhs                    rhs            support confidence lift      count
[1]  {柿ピー}            => {ビール}       0.3     1           2.500000 3

を解釈するには リフト値とは何？ Weblio辞書 を読むと良いです。すると、一行目は

  • support ： 全体の30%が柿ピーとビールを一緒に購入している。
  • confidence：柿ピー購入者の100%が柿ピーとビールを一緒に購入している。
  • lift：ビールをいきなり購入するよりも、柿ピーを買ってビールを購入する確率が2.5倍大きい。

と解釈できます。ちなみにcountの3をtransactions数の10で割るとsupportとなっています。

おわりに

とりあえず動かすことができたら、support(支持度)・confidence(確信度)・lift(リフト値)について調べて見ると良いかと思われます。自分もまさに調べているところですので参考URLをご参照あれ。併売率・購買率などの用語で調べてもだいたいアソシエーション分析についての記事が見つかるので併せてどうぞ。

参考

2019年6月 統計検定準一級 問11

復習

問題の超大雑把な概要

詳細はお手元に問題を用意してください。

定数の情報：

\theta,\phi \in [0,1], \phi + \theta \le 1

推移確率行列M：

\begin{equation}M=\begin{pmatrix} 1-\theta & \theta & 0 \\\\ \theta & 1-\theta-\phi & \phi \\\\ 0 & \phi & 1-\phi \end{pmatrix} \end{equation}

問題

[1]

格付けの状態(A,B,C)を番号(1,2,3)に対応させ、 i,j \in \{1,2,3\} とする。ある年の格付けAが100社・格付けBが20社・格付けCが0社。翌年において、
・A→Bに推移した企業が5社
・B→Aに推移した企業が1社
・他に変化はなかった
・\phi = 0.01とする

このとき、\thetaの最尤推定値を、小数点第3位を四捨五入して答えよ。

[2]

Mの固有値をλとし、直交行列(U)_{ij}= u_{ij}を
\begin{equation}U^{T}MU = \begin{pmatrix} \lambda_{1} & 0 & 0 \\\\ 0 & \lambda_{2} & 0 \\\\ 0 & 0 & \lambda_{3} \end{pmatrix} \end{equation}
を満たすように取る。
t年に於いて格付けAの企業がt+nにおいて格付けCである確率を、n・λ・uを用いて表せ。

解答

[1]

ある年→翌年の格付けの変化は、

\begin{pmatrix} 95 & 5 & 0 \\\\ 1 & 19 & 0 \\\\ 0 & 0 & 0 \end{pmatrix}

と表される。求めたいのはθの最尤推定値である。確率は与えられている。確率分布は３項の状態を扱っているので多項分布の３項版とする。


多項分布の式は、\displaystyle Mult( m | \mu, N ) = \frac{N!}{m_{1}! \cdots m_{K}! } \prod_{k=1}^{K}\mu_{k}^{m_k} ただし N=\sum_{k=1}^K m_k

よって、尤度関数Lは

\displaystyle L = \frac{ (95+5+0)! }{95!5!0!} (1-\theta)^{95} \theta^{5} 0^0 \times \frac{ (1+19+0)! }{1!19!0!} \theta^1 (1-\theta-\phi)^{19} \phi^0 \times \frac{ (0+0+0)! } {0!0!0!} 0^0 \phi^0 (1-\phi)^0


\phi = 0.01*1を代入しながら整理すると、

\displaystyle L = \frac{ (100)! }{95!5!0!} (1-\theta)^{95} \theta^{5} 0^0 \times \frac{ (20)! }{1!19!0!} \theta^1 (1-\theta-0.01)^{19} (0.01)^0 \times \frac{ (0+0+0)! } {0!0!0!} 0^0 (0.01)^0 (1-0.01)^0


さらに、 t^0 = 0^0 = 0! = 1 , t \in \mathbb{R}から、

\displaystyle L = \frac{ 100! }{95!5!} (1-\theta)^{95} \theta^{5} \times \frac{ 20! }{1!19!} \theta (0.99-\theta)^{19}

と簡約化できる。


従って、\thetaの最尤推定値を \displaystyle \frac{\partial}{\partial \theta} \log L = 0 な θで導出できると信ずれば*2、

\displaystyle \frac{ \partial } {\partial \theta } \log L = 0 から、

\displaystyle \frac{ \partial } {\partial \theta } \log L = \frac{ \partial } {\partial \theta } \Bigl( \log(100!) - \log(95!) - \log(5!) + 95\log(1-\theta) + 5\log(\theta) + \log(20!) - \log(1!) - \log(19!) + \log(\theta) + 19\log(0.99 - \theta) \Bigr) = 0 \\\\ \displaystyle \iff 95 \frac{-1}{1-\theta} + 6 \frac{1}{\theta} + 19\frac{-1}{0.99 - \theta}= 0 \\\\ \displaystyle \iff 120\theta^2 -124.99\theta + 5.94 = 0 \\\\ \displaystyle \iff \theta = 0.991667 \ or \ 0.0499159

ここで、\phi + \theta \le 1 \iff 0.01 + \theta \le 1 \iff \theta \le 0.99 から、 \theta = 0.0499159 \simeq 0.05

となる。

[2]

チャップマン–コルモゴロフ方程式 より、推移確率行列Mのn乗はn回の推移の確率に等しいので、

  • M^nを計算し
  • 格付けA → 格付けCに対応するM^n上の成分すなわち
  • M^nの1行3列目の成分を求め

れば良い。さて、

\begin{equation}U^{T}MU = \begin{pmatrix} \lambda_{1} & 0 & 0 \\\\ 0 & \lambda_{2} & 0 \\\\ 0 & 0 & \lambda_{3} \end{pmatrix} \end{equation}

の両辺をn乗すると、

\begin{equation}U^{T}M^nU = \begin{pmatrix} \lambda_{1}^n & 0 & 0 \\\\ 0 & \lambda_{2}^n & 0 \\\\ 0 & 0 & \lambda_{3}^n \end{pmatrix} \end{equation}

となるので、

\begin{equation} M^n = U \begin{pmatrix} \lambda_{1}^n & 0 & 0 \\\\ 0 & \lambda_{2}^n & 0 \\\\ 0 & 0 & \lambda_{3}^n \end{pmatrix} \end{equation} U^T

従って、

\begin{equation} M^n = \begin{pmatrix} u_{11} & u_{12} & u_{13} \\\\ u_{21} & u_{22} & u_{23} \\\\ u_{31} & u_{32} & u_{33} \end{pmatrix} \begin{pmatrix} \lambda_{1}^n & 0 & 0 \\\\ 0 & \lambda_{2}^n & 0 \\\\ 0 & 0 & \lambda_{3}^n \end{pmatrix} \end{equation} \begin{pmatrix} u_{11} & u_{12} & u_{13} \\\\ u_{21} & u_{22} & u_{23} \\\\ u_{31} & u_{32} & u_{33} \end{pmatrix} ^T

\iff \begin{equation} M^n = \begin{pmatrix} u_{11} & u_{12} & u_{13} \\\\ u_{21} & u_{22} & u_{23} \\\\ u_{31} & u_{32} & u_{33} \end{pmatrix} \begin{pmatrix} \lambda_{1}^n & 0 & 0 \\\\ 0 & \lambda_{2}^n & 0 \\\\ 0 & 0 & \lambda_{3}^n \end{pmatrix} \end{equation} \begin{pmatrix} u_{11} & u_{21} & u_{31} \\\\ u_{12} & u_{22} & u_{32} \\\\ u_{13} & u_{23} & u_{33} \end{pmatrix}

となるので、M^nの1行3列目の成分は、

\begin{pmatrix} \lambda_{1}^n u_{11} & \lambda_{2}^n u_{12} & \lambda_{3}^n u_{13} \end{pmatrix} \begin{pmatrix} u_{31} & u_{32} & u_{33} \end{pmatrix}^T \\\\ = \displaystyle \lambda_{1}^n u_{11}u_{31} + \lambda_{2}^n u_{12}u_{32} + \lambda_{3}^n u_{13}u_{33} \\\\ = \displaystyle \sum_{k=1}^{3} \lambda_{k}^n u_{1k}u_{3k}


である。

*1:ここで尤度関数の変数が(θ,φ)からθのみになるのだと思われる。慈悲に感謝。

*2:真面目にやる場合はラグランジュの未定なんちゃら法とか極値が極大値になっているかとかを調べる必要があるゆえ。然し乍らどうやら信ずるものは救われるらしく、たとえば多項分布の推定 - MOXBOXを参照。

axjack is said to be an abbreviation for An eXistent JApanese Cool Klutz.