有没有人用过kernel density estimate?
本帖最后由 likm1110 于 2015-5-16 11:52 编辑Kernel density estiamte 的程序谁会？
我参考的WIKIhttp://en.wikipedia.org/wiki/File:Comparison_of_1D_histogram_and_KDE.png#/media/File:Comparison_of_1D_histogram_and_KDE.png
http://bbs.fcode.cn/forum.php?mod=image&aid=727&size=300x300&key=ebd6409e127c6d46&nocache=yes&type=fixnone
第二张图是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
我画出来的图形状差不多，但是最大值不一样。
代码：
!===============================================================
! 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（图贴不上来....）
http://bbs.fcode.cn/forum.php?mod=image&aid=734&size=300x300&key=ee749bb64a2b0340&nocache=yes&type=fixnone
页:
[1]