有没有人用过kernel density estimate?

likm1110

Kernel density estiamte 的程序谁会？
我参考的WIKI

第二张图是kernael estimate.这个图究竟是怎么画的？
第二张图的说明：
‘For the kernel density estimate, we place a normal kernel with variance 2.25 (indicated by the red dashed lines) on each of the data points xi. The kernels are summed to make the kernel density estimate (solid blue curve). The smoothness of the kernel density estimate is evident compared to the discreteness of the histogram, as kernel density estimates converge faster to the true underlying density for continuous random variables’

Wiki连接：http://en.wikipedia.org/wiki/Kernel_density_estimation

我画出来的图形状差不多，但是最大值不一样。

代码：

[Fortran] syntaxhighlighter_viewsource syntaxhighlighter_copycode

!===============================================================
! Description : The Histogram and KDE
! Input : 1. n     number of data points
!         2. x(n)  the data points
! File : 1. Hist.txt    file of Histogram datas
!        2. KDE.txt     file of KDE datas
!===============================================================
Module EstimatePack
Contains
  !---------------------------------------------------
  ! Purpose : The data of histogram
  !---------------------------------------------------
  Subroutine Hist(n,x)
    Implicit None
    Integer :: n
    Real :: x(n)
    Real, Allocatable :: f_hist(:)
    Integer, Allocatable :: num(:)
    Real :: xmin, xmax
    Real :: dbin
    Integer :: xlen, xrange, xranget
    Integer :: i, ip
    dbin = 2.
    xmin = x(1)
    xmax = x(1)
    Do i=1,n
      If (xmin>x(i)) Then
        xmin = x(i)
      Endif
      If (xmax<x(i)) Then
        xmax = x(i)
      Endif
    Enddo
    xmin = xmin-dbin
    xmax = xmax+dbin
    !Write(*,*) 'xmin =', xmin
    !Write(*,*) 'xmax =', xmax
    xlen = int((xmax-xmin)/dbin)
    Allocate(num(xlen),f_hist(xlen))
    f_hist = 0.
    num = 0
    xranget = int(xmin)
    Open(100,file='t_Hist.txt')
    Do i = 1,xlen
      xrange = xranget+int(dbin)
      num(i) = count (xranget<x(:).and.x(:)<xrange)
      f_hist(i) = num(i)*(1./12.)    
      Write(100,*) real(xranget+xrange)/2.,f_hist(i)
      xranget = xranget+int(dbin)
    Enddo
    Close(100)
    Deallocate(num,f_hist)
  End Subroutine Hist
  !---------------------------------------------------
  ! Purpose : function of Gauss Distribution
  ! Input : 1. x        the independent variable
  !         2. mui      the mean(in fact the data point selected)
  !         3. segma    the standard deviation
  !---------------------------------------------------
  Function Fai(x,mui,segma)
    Implicit None
    Real, Parameter :: pi = 3.1415926
    Real, Intent(in) :: x
    Real :: x_t
    Real, Intent(in) :: mui, segma
    Real :: Fai
    x_t = (x-mui)/segma
    Fai = (1./segma)*exp(-x_t**2/2.)/sqrt(2*pi)
  End Function Fai
  !---------------------------------------------------
  ! Purpose : The data of KDE_component
  !---------------------------------------------------
  Subroutine KDE(n,x)
    Implicit None
    Integer :: n
    Real :: x(n)
    Real, Parameter :: xvar = 3.
    Real :: segma = sqrt(2.25)
    Integer :: i, j, step
    Real :: x_in
    Real :: dbin
    Real :: h
    Character(1) :: str 
    Real, Allocatable :: x_kde_c(:,:),f_kde_c(:,:)
    dbin = 0.1
    h = 1.06*segma*(n)**(-0.2)
    step = int(2*xvar/dbin)+1
    Allocate(x_kde_c(n,step),f_kde_c(n,step))
    Do i=1,n
      Write(str,'(I1)') i
      Open(101,file='t_kde_'//str//'.txt')
      x_in = x(i)-xvar
      Do j=1,step   
        Write(101,*) x_in, Fai(x_in,x(i),segma)/(h*real(n))   
        x_kde_c(i,j) = x_in
        f_kde_c(i,j) = Fai(x_in,x(i),segma)/(h*real(n))
        x_in = x_in+dbin
      Enddo
      Close(101)
    Enddo
    !Write(*,*) x_kde_c(2,10),f_kde_c(2,10) 
    Open(102,file='t_kde_x.txt')
    Open(103,file='t_kde_f.txt')
    Do i=1,n
      Do j=1,step
        Write(102,*) i,j,x_kde_c(i,j)
        Write(103,*) i,j,f_kde_c(i,j)
      Enddo
    Enddo
    Close(102)
    Close(103)
    Deallocate(x_kde_c,f_kde_c)
  End Subroutine KDE
End Module EstimatePack
!===============================================================
! Description : The main program
! Post Script : 1. input n(data points dimension) and x(data points)
!===============================================================
Program main
Use EstimatePack
  Implicit None
  Integer, Parameter :: n=6
  Integer :: step, numx, numf, num
  Integer :: ip, ipi, ipj, i, j, iptr
  Real :: temp, tempf
  Real :: x(n)
  Integer :: ios
  Real, Allocatable :: x_kde_temp(:), f_kde_temp(:)
  Real :: x_kde, f_kde
  x = (/-2.1,-1.3,-0.4,1.9,5.1,6.2/)
  numx = 0
  numf = 0
  Call Hist(n,x)
  Call KDE(n,x)
  Open(11,file='t_kde_x.txt')
  Do 
    Read(11,*,iostat=ios) 
    numx = numx+1
    If (ios/=0) Exit
  Enddo
  Close(11)
  Open(12,file='t_kde_f.txt')
  Do 
    Read(12,*,iostat=ios) 
    numf = numf+1
    If (ios/=0) Exit
  Enddo
  Close(12)
  !Write(*,*) numx, numf
  If (numf/=numx) Then
    Write(*,*) 'data pointa have error'
    Else
    num = numx-1
  Endif
  Allocate(x_kde_temp(num), f_kde_temp(num))
  Open(13,file='t_kde_x.txt')
  Do i=1,num
    Read(13,*) ipi, ipj, x_kde_temp(i)
  Enddo
  Close(13)
  Open(14,file='t_kde_f.txt')
  Do i=1,num
    Read(14,*) ipi, ipj, f_kde_temp(i)
  Enddo
  Close(14)
  outer: Do i = 1, num-1
    iptr = i
    inner: Do j =i+1, num
      minval: If(x_kde_temp(j) < x_kde_temp(iptr)) Then
        iptr = j
      Endif minval
    Enddo inner
    swap: If(i /= iptr) Then
      temp = x_kde_temp(i)
      x_kde_temp(i) = x_kde_temp(iptr)
      x_kde_temp(iptr) = temp
      tempf = f_kde_temp(i)
      f_kde_temp(i) = f_kde_temp(iptr)
      f_kde_temp(iptr) = tempf
    Endif swap
  Enddo outer  
  Open(15,file='t_kde_xf.txt')
  Do i=1,num
    Write(15,*) i, x_kde_temp(i), f_kde_temp(i) 
  Enddo
  Close(15)
  Open(16,file='t_kde.txt')
  ip = 0
  Do i=1,num
    If (i<ip.or.i==ip) Cycle
    x_kde = x_kde_temp(i)
    f_kde = f_kde_temp(i)
    Do j=i+1,num
      If (abs(x_kde-x_kde_temp(j))<0.0001) Then
        x_kde = x_kde_temp(i)
        f_kde = f_kde + f_kde_temp(j)
        ip = j
      Endif
    Enddo
    Write(16,*) x_kde, f_kde
  Enddo
  Close(16)
  Deallocate(x_kde_temp, f_kde_temp)
End Program main

画图用的是't_kde.txt'的数据，得到蓝色实线，我的最大值明显比原图小0.2-0.3（图贴不上来....）