R语言学习系列(画向量的密度直方图)

利用R语言画密度直方图比较方便，但为了理解密度函数的意义和如何计算密度值，于是用plot,lines两个画图函数来自己实现密度直

方图的画图程序脚本如下：

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[a<X<=b]，一般直方图x轴采用相等分段方式，前面的bw值就是x轴上直方图的长，高代表落于该直方图的向量个数m.显然有密度d满足：d*m*bw=m/n (P[a<X<=b]) -> d=1/(n*w).从这也可以看出核密度直方图对应的概率密度函数就是一种分段函数

d(x)=mi/n * x(a+(i-1)*bw <x<=a+i * bw, mi为落在i段内的向量个数,n为向量总数,bw是每个直方图想x轴上的长度).
可以很轻易的得出随机变量X的分布函数F(x)(当然也需要分段)，这里不好写表达式，大家自己可以想一下.从上面可以看出，如果采样样本值足够多，分段bw足够小的时候，d(x)就可以不断逼近真实的概率密度函数。
 以上程序在R.2.15上测试通过，并与hist画图对比，图形类似.