「統計学関連なんでもあり」の過去ログ--- 048

No.21922　AIC による，ヒストグラム（度数分布表）の最適階級分割の探索　　【明石】　2016/03/06(Sun) 20:27

青木先生様，
お忙しいところを失礼いたします，明石と申します。

「AIC による，ヒストグラム（度数分布表）の最適階級分割の探索」
　http://aoki2.si.gunma-u.ac.jp/R/AIC-Histogram.html
について，ご教示をいただきたいと思います。
よろしくお願いいたします。

−−−
例示では，ヒストグラムのビンが19本あります。

(1)
最初の2本⇒併合
次の5本　⇒併合
次の5本　⇒併合
次の7本　⇒併合

結果として，4階級になるという理解でよろしいでしょうか？

(2)
ans$df[1:10,]の読み方が分かりません。

AICの値の小さい順に上位10個を表示していると理解していますが，
c1 r c2 の読み方が，分かりません。

出力結果の見方を，教えていただけないでしょうか？

    c1 r c2      AIC
78   2 5  7 487.5261
60   2 5  2 488.7143
105  3 5  1 490.2904
72   2 6  5 490.6599
130  3 4  8 491.1111
59   2 3  2 491.6356
123  3 5  6 492.0030
116  3 4  4 492.4315
55   2 4  1 492.9979
129  3 2  8 493.2800

(3)
ヒストグラムと ans$df[1:10,] を一緒に見ています。

ヒストグラム（19本のビンを，4階級に併合）の区間を知りたい場合には，
ansのどの要素を参照すればよろしいのでしょうか？

(4)
数量データの離散化については，
値等分割，レコード数等分割（四分位など），平均と分散
などの方法があります。

本方式「AIC による，ヒストグラム（度数分布表）の最適階級分割の探索」
は，数量データの離散化の一つの手法であるという理解でよろしいでしょうか？

ご教示を，よろしくお願いいたします。

No.21923　Re: AIC による，ヒストグラム（度数分布表）の最適階級分割の探索　　【青木繁伸】　2016/03/07(Mon) 09:31

(1) 結果として，4階級になるという理解

その通り

(2) c1 r c2 の読み方が，分かりません。

プログラム中の注釈を参照
for (c1 in 1:(c-2)) { # 左端で併合する階級数を探索
for (c2 in 1:(c-2)) { # 右端で併合する階級数を探索
:
for (r in 1:(c-c1-c2)) { # 中央付近で併合する階級数を探索

(3) ヒストグラム（19本のビンを，4階級に併合）の区間を知りたい

df 要素ではなく，breaks 要素を見ます。ans$breaks

(4) 数量データの離散化の一つの手法

その通り

No.21924　【御礼】 Re: AIC による，ヒストグラム（度数分布表）の最適階級分割の探索　　【明石】　2016/03/07(Mon) 09:37

青木先生様，
お忙しいところを失礼いたします，明石と申します。

ご教示いただき，大変に助かりました。

今まで，EMアルゴリズム（混合分布）を用いて，離散化していましたが，
うまくいかないことが多く，苦慮していました。

心から，心より，御礼を申し上げます。

No.21928　Re: AIC による，ヒストグラム（度数分布表）の最適階級分割の探索　　【明石】　2016/03/07(Mon) 21:50

青木先生　様；

お忙しいところを失礼いたします，明石と申します。
お世話になります。

追加の質問がございます。
ご教示を何卒どうぞ，よろしくお願いいたします。

(4)
ans$n

> [1] 16 22 11 6 7 9 11 4 3 2 3 2 1 1 0 0 0 0 2

ビン（19本）の高さ，と理解しました。

(5)

ans$breaks

> [1] 0.015000 5.369737 10.724474 16.079211 21.433947 26.788684 32.143421 37.498158 42.852895
[10] 48.207632 53.562368 58.917105 64.271842 69.626579 74.981316 80.336053 85.690789 91.045526
[19] 96.400263 101.755000

ビン（19本）の区間の左端，右端の値，と理解しました。

(6)
.frame': 416 obs.

⇒　収束計算の繰り返し回数と理解しました。

(7)
併合した結果の4階級の「4」という値を取得するには，どうすればよいのでしょうか？

(8)
4階級を，左から，1,2,3,4 と番号付けしたといたします。
ベクトルxの要素（要素数は100）が，どの階級（1〜4）に属するのかを出力する
Rプログラムを考えていますが，作成できません。
誠に申し訳ございませんが，Rコードを例示していただけると助かります。

何卒どうぞ，ご教示をよろしくお願いいたします。

No.21929　Re: AIC による，ヒストグラム（度数分布表）の最適階級分割の探索　　【青木繁伸】　2016/03/07(Mon) 22:52

併合後のヒストグラム（の長方形）を描くのは，以下の部分で，その長方形の X 座標は，lo+cmg[i]*delta，lo+cmg[i+1]*delta
なので，	sapply(1:(m+2), function(i) rect(lo+cmg[i]*delta, 0,            # 併合後のヒストグラム
		lo+cmg[i+1]*delta, mean(p[(cmg[i]+1):(cmg[i+1])]),
		border="red", col="red", density=15))
それを，cat でコンソールに書き出すとすれば<
	sapply(1:(m+2), function(i) { rect(lo+cmg[i]*delta, 0,            # 併合後のヒストグラム
		lo+cmg[i+1]*delta, mean(p[(cmg[i]+1):(cmg[i+1])]),
		border="red", col="red", density=15)
		cat(lo+cmg[i]*delta, lo+cmg[i+1]*delta, "\n") # これを追加しますね。カテゴリー数は m＋2 かな？
		})

> ans <- AIC.Histogram(x, 0.01)
0.015 10.72447 
10.72447 37.49816 
37.49816 64.27184 
64.27184 101.755

No.21930　【Special Thanks !】 Re: AIC による，ヒストグラム（度数分布表）の最適階級分割の探索　　【明石】　2016/03/08(Tue) 05:11

青木先生様，

お忙しいところを失礼いたします，明石と申します。

ご教示いただき，大変に助かりました。

ご教示いただきました内容を基にして，自分の頭でトレースをいたします。

心から，心より，御礼を申し上げます。

No.21934　Re: AIC による，ヒストグラム（度数分布表）の最適階級分割の探索　　【明石】　2016/03/11(Fri) 16:29

青木先生様，
お忙しいところを失礼いたします，明石と申します。

お世話になります。
ご教示いただきました内容を追試して，理解をしました。

最適階級の区間（左端の値，右端の値）を，コンソールに出力できるようになりましたが，
関数の戻り値に，加えたいと思い，
先生からご教示いただきましたコードに下記のように加筆しましたが，
うまくいきません。

ご教示をいただければ，助かります。
どうぞ，よろしくお願いいたします。
k <- m+2   #　最適階級数

mx <- matrix(0, nrow=k, ncol=2)

sapply(1:(m+2), function(i) {
	rect(lo+cmg[i]*delta, 0,            # 併合後のヒストグラム
	lo+cmg[i+1]*delta, mean(p[(cmg[i]+1):(cmg[i+1])]),
	border="red", col="red", density=15)
	r1 <- lo+cmg[i]*delta
	r2 <- lo+cmg[i+1]*delta
	cat(r1, r2, "\n")

	mx[i,1] <- r1
	mx[i,2] <- r2
	})

invisible(list(df=df, n=n, breaks=brks, k=k, mx=mx))  # 結果を返す
以上，よろしくお願いいたします。

No.21935　Re: AIC による，ヒストグラム（度数分布表）の最適階級分割の探索　　【青木繁伸】　2016/03/11(Fri) 16:58

一番少ない修正量で動かすには

mx[i,1] <<- r1
mx[i,2] <<- r2

しかし，<<- はなるべく使わないほうがよい

正統的には

#       k = m+2                          # 使わない     
#       mx <- matrix(0, nrow=k, ncol=2)  # 使わない
        cmg <- cumsum(c(0, c1, rep(r, m), c2))
        mx <- sapply(1:(m+2), function(i) { rect(lo+cmg[i]*delta, 0, # 結果を mx に代入
                lo+cmg[i+1]*delta, mean(p[(cmg[i]+1):(cmg[i+1])]),
                border="red", col="red", density=15)
                r1 <- lo+cmg[i]*delta
                r2 <- lo+cmg[i+1]*delta
                cat(i, r1, r2, "\n")
                return(c(r1, r2)) # sapply から返される値
        })
        mtext(paste("AIC =", round(AIC, 2)), side=3, line=-1.2, 
                at=lo+0.8*c*delta)
        invisible(list(df=df, n=n, breaks=brks, k=ncol(mx), mx=t(mx)))

No.21936　【御礼】 Re: AIC による，ヒストグラム（度数分布表）の最適階級分割の探索　　【明石】　2016/03/11(Fri) 20:40

青木先生様，

お忙しいところを失礼いたします，明石と申します。

今回も，大変に良い勉強をさせていただき，誠にありがとうございました。

ご教示いただきました内容を基にして，自分の頭でトレースをいたします。

心から，心より，御礼を申し上げます。