[其他行业算法] 有没有人用过kernel density estimate?

!===============================================================
! 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