R言語学習シリーズ(ベクトルの密度ヒストグラムを描く)


R言語で密度ヒストグラムを描くのは便利だが,密度関数の意味と密度値をどのように計算するかを理解するためにplot,linesの2つのグラフ関数を用いて自分で密度直を実現する.
方図のグラフィックプログラムのスクリプトは以下の通りです.
DrawDensity = function(x,bw=5)
{
   if(any(bw<=0))
   {
     bw <- 5
    
   }
   #print(bw)
   n <- length(x)
   if(any(n>0))
   {
     mn <- min(x)
     mx <- max(x)
     rmn <- (mn %/% bw)
     begx <-  rmn * bw
     rmn <- (mx %/% bw)
     endx <- (rmn) * bw
     if ((mx-rmn * bw) > 0)
     {
        endx <- (rmn+1) * bw
     }
     
     segs <- (endx-begx)  %/% bw
     segcounts <- numeric(segs)
     validsegs <- 0
     print(segs)
     for(i in 1:segs)
     {
        slw <- begx+(i-1) * bw
        sup <- begx+(i) * bw
        tmp<- x[x > slw & x<=sup];
        segcounts[i]<- length(tmp)
        if(length(tmp)>0)
        {
           validsegs <- validsegs+1
        }
     }
     #      :               .
     baseIndesity = 1 / (n*bw);
     ymax = baseIndesity * (max(segcounts)+1)
     plot(c(begx,endx),c(0,0),xlim=c(begx,endx),ylim=c(0,ymax))
     lines(c(begx,endx),c(0,0))
     lines(c(0,0),c(0,ymax))
     #print(segs)
     for (i in 1:segs)
     {
       x1 <- begx+(i-1) * bw
       x2 <- begx+(i) * bw
       y<- baseIndesity  * segcounts[i]
       drawHist(x1,x2,y)
       #print(i)
     }

   }
}
drawHist = function(x1,x2,y)
{
  #print(x1);print(x2);print(y)
  lines(c(x1,x1),c(0,y))
  lines(c(x1,x2),c(y,y))
  lines(c(x2,x2),c(0,y))
}

核密度関数を理解する鍵は単位密度を理解する計算式:1/(n*bw)である.ヒストグラムを描くにはセグメントが必要であるため、各ヒストグラムの面積はそのパターン領域に落ちる確率P[ad=1/(n*w)である.これからも核密度ヒストグラムに対応する確率密度関数がセグメント関数であることがわかる.
d(x)=mi/n*x(a+(i-1)*bw