階級数を変えてヒストグラムを描く関数

ヒストグラムは図のように、階級数の違いによって、分布の概観がかなり異なってみえる。この事例では階級数が13の場合は正規分布のように見えるが、階級数が33の場合は右側に盛り上がりがあるようにみえる。
どのように階級を設定すべきかについては、決定的なものはないが、一般的に良く使われている階級幅はスタージェス(Sturges)の公式と呼ばれるものから計算したもので、データ数をN、階級数をLしたときに
L= floor(1 + \log_2N)
で計算されるものである。ここで、floor(x)は実数xに対してx以下の最大の整数を与える関数の意味で用いている。Rのデフォルトの階級数もこの設定である。一方、探索的データ解析では、階級数の最大値として、
L=floor(10\log_{10}N)
を推奨している。この方法はデータ数が少ない場合に、より多めの階級数を取ることになる。いずれにしても最適な階級数があるわけでゃなく、階級数を変化させて分布形状を調べるのが良いとされているようである。この探索的データ解析の方法をデフォルトとして、階級数をかえてヒストグラムを描くRの関数を作成した。

#-----------------------------------
# 階級数を変えてヒストグラムを求める関数
#-----------------------------------
exhist <- function(		x, 				#xはデータベクトル
                       	n=NULL,    	# 階級数の指定
                        ...)          	# hist に引き渡す任意の引数
{
	N <- length(x)				#データ数を求める
	scale_n <- floor(10*log10(N))		#階級数を求める
	#階級数が与えられていたら
		if(!is.null(n)){		
			scale_n <- n 
		}
	digits_n <- max(nchar(x-trunc(x)))	#有効桁数を求める
		if (digits_n > 2){ 
			digits_n <- digits_n -2
		}else{
			digits_n <- 0
		}
	#階級幅を求める
	M <- round((max(x)-min(x))/scale_n,digits=digits_n)
	#スタートポイントとエンドポイントを求める
	SP <- round((min(x)-M/2),digits=digits_n)
	EP <- round((max(x)+M/2),digits=digits_n)
	#階級ベクトルを設定する
	brk=seq(SP,EP,length=scale_n+1)
	#ヒストグラムを描く
	hist(x,breaks=brk,...)
	#階級ベクトルを返す
	return(brk)
}