2014/07/21
データいろいろ > 階級幅だけ決めてヒストグラムを描く
データいろいろ > 階級幅だけ決めてヒストグラムを描く
R の hist 関数はヒストグラムを描けて便利な一方、グラフィックパラメータを上書きしてしまうようでレイアウト調整がしにくい。また引数 breaks で全階級の境界点を渡す必要があり、主に seq 関数で始点・終点・階級幅の三つを渡す。だからデータの値範囲が変わると breaks の始点・終点を変える必要がある。同じ系統のデータで裾がちょっと伸びたような時、面倒。
そこで階級幅だけ決めて、後はデータ範囲から自動的に breaks に渡す中身を生成する短いコードを書いたが、ついでに階級ごとの計数（そこに含まれるデータ個数）を算出して plot + rect で描画するようにした。また plot の前までと同じことを PostgreSQL のクエリでも書いた。実行環境は 2014/04/17 を参照。
↓ テンプレートと実行結果。サンプルデータには rnorm 関数による正規乱数を使った。設定するのは基本的にデータと階級幅（変数 wid）だけ。
↓ コードの中で作られるヒストグラム用データ。blw と upp が各階級の上下限。cnt が計数。このデータから rect 関数で細長い矩形を描画している。
↓ 二つデータ・階級幅での例に、hist 関数のデフォルト描画を併記。最初のデータは上と同じ rnorm(10000)。次に全体を 100 倍し、合わせて階級幅を変えた。
同じことの描画の前まで、つまり各階級の上下限と計数を出すところまでを PostgreSQL でも行う。まず正規乱数のサンプルデータとして、2013/11/23（ボックス=ミュラー法を SQL で）と同様の式で ↓ 一時テーブルを作成。
このテーブルに対し、階級幅だけ設定して ↓ を実行するとデータ範囲から自動的に全階級を出してカウントする。実行結果も合わせて示す。
↓ 同じクエリで、元データと階級幅ともに 100 倍して実行。クエリの変更箇所はブロック a と b だけ。実行の結果、階級の上下限だけが変わる。
↓ 確認のため、上の元テーブルを R で読み込み、先ほどと同様にヒストグラムを描画。各階級のデータは SQL の結果と同じ。
そこで階級幅だけ決めて、後はデータ範囲から自動的に breaks に渡す中身を生成する短いコードを書いたが、ついでに階級ごとの計数（そこに含まれるデータ個数）を算出して plot + rect で描画するようにした。また plot の前までと同じことを PostgreSQL のクエリでも書いた。実行環境は 2014/04/17 を参照。
↓ テンプレートと実行結果。サンプルデータには rnorm 関数による正規乱数を使った。設定するのは基本的にデータと階級幅（変数 wid）だけ。
# サンプルデータ（正規乱数）
dat = rnorm(10000)
# 階級幅
wid = 0.25
# 各階級の計数を算出
tmp = floor(range(dat) / wid ) * wid
seq = seq(tmp[1], tmp[2], by=wid)
mat = cbind(blw=seq, upp=seq+wid, cnt=NA)
mat[,'cnt'] = apply(mat, 1, function(x){
length(which(x['blw'] <= dat & dat < x['upp'])) })
# 低水準作図関数で描画
windows(height=344/96, width=436/96)
par(mfrow=c(1,1), plt=c(15, 95, 15, 95)/100)
plot(0, type='n', xlim=c(min(mat[,'blw']), max(mat[,'upp'])), xlab=''
, ylim=c(0,max(mat[,'cnt']) * 1.05), yaxs='i', ylab='')
rect(mat[,'blw'],0,mat[,'upp'],mat[,'cnt'], border=1, col='gray')
↓ コードの中で作られるヒストグラム用データ。blw と upp が各階級の上下限。cnt が計数。このデータから rect 関数で細長い矩形を描画している。
↓ 二つデータ・階級幅での例に、hist 関数のデフォルト描画を併記。最初のデータは上と同じ rnorm(10000)。次に全体を 100 倍し、合わせて階級幅を変えた。
windows(height=1000/96, width=445/96)
par(mfrow=c(4,1), omd=c(10, 95, 5, 100)/100, plt=c(0, 100, 5, 85)/100)
cof = c(1, 100) # 元データに掛ける数
wid = c(0.25, 20) # 階級幅
for(i in 1:2){
dat = rnorm(10000) * cof[i]
hist(dat, main='', xlab='', ylab='') # 通常のヒストグラム
# 区分値と計数を作成
tmp = floor(range(dat) / wid[i] ) * wid[i]
seq = seq(tmp[1], tmp[2], by=wid[i])
mat = cbind(blw=seq, upp=seq+wid[i], cnt=NA)
mat[,'cnt'] = apply(mat, 1, function(x){
length(which(x['blw'] <= dat & dat < x['upp'])) })
# 低水準作図関数で描画
plot(0, type='n', xlim=c(min(mat[,'blw']), max(mat[,'upp'])), xlab=''
, ylim=c(0,max(mat[,'cnt']) * 1.05), yaxs='i', ylab='')
rect(mat[,'blw'],0,mat[,'upp'],mat[,'cnt'], border=1, col='gray')
}
同じことの描画の前まで、つまり各階級の上下限と計数を出すところまでを PostgreSQL でも行う。まず正規乱数のサンプルデータとして、2013/11/23（ボックス=ミュラー法を SQL で）と同様の式で ↓ 一時テーブルを作成。
CREATE TABLE "201407"."21_sample_norm" AS
SELECT sqrt(-2 * ln(random())) * sin(2 * pi() * random()) val
FROM generate_series(1, (10^4) :: int) ;
このテーブルに対し、階級幅だけ設定して ↓ を実行するとデータ範囲から自動的に全階級を出してカウントする。実行結果も合わせて示す。
WITH a AS (
-- サンプルデータ
SELECT val FROM "201407"."21_sample_norm"
), b AS (
-- 階級幅を設定
SELECT 0.25 :: float wid
), c AS (
-- データ範囲（最小値・最大値を縦に積む）
SELECT unnest(ARRAY[min(val), max(val)]) rng FROM a
), d AS (
-- 階級幅で下側に丸めた最小値 ※画像にある「最小最大値」は誤記
SELECT floor(rng / wid) * wid tmp FROM b, c
), e AS (
-- 各階級準備
SELECT wid, min(tmp) + wid * generate_series(0
, ((max(tmp) - min(tmp)) / wid) :: int) blw
FROM b, d GROUP BY wid
), f AS (
-- 各階級
SELECT blw, blw + wid upp FROM e
)
SELECT blw, upp, count(*) FROM a JOIN f ON blw <= val AND val < upp
GROUP BY blw, upp ORDER BY 1 ;
↓ 同じクエリで、元データと階級幅ともに 100 倍して実行。クエリの変更箇所はブロック a と b だけ。実行の結果、階級の上下限だけが変わる。
↓ 確認のため、上の元テーブルを R で読み込み、先ほどと同様にヒストグラムを描画。各階級のデータは SQL の結果と同じ。
library(RPostgreSQL)
con = dbConnect(PostgreSQL(), dbname='daily', port=5434, user='xxx', password='xxx')
dat = dbGetQuery(con, 'SELECT val FROM "201407"."21_sample_norm"')[,1]
dbDisconnect(con)
wid = 0.25
tmp = floor(range(dat) / wid ) * wid
seq = seq(tmp[1], tmp[2], by=wid)
mat = cbind(blw=seq, upp=seq+wid, cnt=NA)
mat[,'cnt'] = apply(mat, 1, function(x){
length(which(x['blw'] <= dat & dat < x['upp'])) })
print(mat)
windows(height=344/96, width=445/96)
par(mfrow=c(1,1), plt=c(15, 95, 15, 95)/100)
plot(0, type='n', xlim=c(min(mat[,'blw']), max(mat[,'upp'])), xlab=''
, ylim=c(0,max(mat[,'cnt']) * 1.05), yaxs='i', ylab='')
rect(mat[,'blw'],0,mat[,'upp'],mat[,'cnt'], border=1, col='gray')