有没有人用过kernel density estimate?

likm1110 · 发表于 2015-5-15 10:29:17

本帖最后由 likm1110 于 2015-5-16 11:52 编辑

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] 纯文本查看 复制代码

001!===============================================================
002! Description : The Histogram and KDE
003! Input : 1. n     number of data points
004!         2. x(n)  the data points
005! File : 1. Hist.txt    file of Histogram datas
006!        2. KDE.txt     file of KDE datas
007!===============================================================
008Module EstimatePack
009Contains
010  !---------------------------------------------------
011  ! Purpose : The data of histogram
012  !---------------------------------------------------
013  Subroutine Hist(n,x)
014    Implicit None
015    Integer :: n
016    Real :: x(n)
017    Real, Allocatable :: f_hist(:)
018    Integer, Allocatable :: num(:)
019    Real :: xmin, xmax
020    Real :: dbin
021    Integer :: xlen, xrange, xranget
022    Integer :: i, ip
023    dbin = 2.
024    xmin = x(1)
025    xmax = x(1)
026    Do i=1,n
027      If (xmin>x(i)) Then
028        xmin = x(i)
029      Endif
030      If (xmax<x(i)) Then
031        xmax = x(i)
032      Endif
033    Enddo
034    xmin = xmin-dbin
035    xmax = xmax+dbin
036    !Write(*,*) 'xmin =', xmin
037    !Write(*,*) 'xmax =', xmax
038    xlen = int((xmax-xmin)/dbin)
039    Allocate(num(xlen),f_hist(xlen))
040    f_hist = 0.
041    num = 0
042    xranget = int(xmin)
043    Open(100,file='t_Hist.txt')
044    Do i = 1,xlen
045      xrange = xranget+int(dbin)
046      num(i) = count (xranget<x(:).and.x(:)<xrange)
047      f_hist(i) = num(i)*(1./12.)   
048      Write(100,*) real(xranget+xrange)/2.,f_hist(i)
049      xranget = xranget+int(dbin)
050    Enddo
051    Close(100)
052    Deallocate(num,f_hist)
053  End Subroutine Hist
054  !---------------------------------------------------
055  ! Purpose : function of Gauss Distribution
056  ! Input : 1. x        the independent variable
057  !         2. mui      the mean(in fact the data point selected)
058  !         3. segma    the standard deviation
059  !---------------------------------------------------
060  Function Fai(x,mui,segma)
061    Implicit None
062    Real, Parameter :: pi = 3.1415926
063    Real, Intent(in) :: x
064    Real :: x_t
065    Real, Intent(in) :: mui, segma
066    Real :: Fai
067    x_t = (x-mui)/segma
068    Fai = (1./segma)*exp(-x_t**2/2.)/sqrt(2*pi)
069  End Function Fai
070  !---------------------------------------------------
071  ! Purpose : The data of KDE_component
072  !---------------------------------------------------
073  Subroutine KDE(n,x)
074    Implicit None
075    Integer :: n
076    Real :: x(n)
077    Real, Parameter :: xvar = 3.
078    Real :: segma = sqrt(2.25)
079    Integer :: i, j, step
080    Real :: x_in
081    Real :: dbin
082    Real :: h
083    Character(1) :: str 
084    Real, Allocatable :: x_kde_c(:,:),f_kde_c(:,:)
085    dbin = 0.1
086    h = 1.06*segma*(n)**(-0.2)
087    step = int(2*xvar/dbin)+1
088    Allocate(x_kde_c(n,step),f_kde_c(n,step))
089    Do i=1,n
090      Write(str,'(I1)') i
091      Open(101,file='t_kde_'//str//'.txt')
092      x_in = x(i)-xvar
093      Do j=1,step   
094        Write(101,*) x_in, Fai(x_in,x(i),segma)/(h*real(n))  
095        x_kde_c(i,j) = x_in
096        f_kde_c(i,j) = Fai(x_in,x(i),segma)/(h*real(n))
097        x_in = x_in+dbin
098      Enddo
099      Close(101)
100    Enddo
101    !Write(*,*) x_kde_c(2,10),f_kde_c(2,10) 
102    Open(102,file='t_kde_x.txt')
103    Open(103,file='t_kde_f.txt')
104    Do i=1,n
105      Do j=1,step
106        Write(102,*) i,j,x_kde_c(i,j)
107        Write(103,*) i,j,f_kde_c(i,j)
108      Enddo
109    Enddo
110    Close(102)
111    Close(103)
112    Deallocate(x_kde_c,f_kde_c)
113  End Subroutine KDE
114End Module EstimatePack
115!===============================================================
116! Description : The main program
117! Post Script : 1. input n(data points dimension) and x(data points)
118!===============================================================
119Program main
120Use EstimatePack
121  Implicit None
122  Integer, Parameter :: n=6
123  Integer :: step, numx, numf, num
124  Integer :: ip, ipi, ipj, i, j, iptr
125  Real :: temp, tempf
126  Real :: x(n)
127  Integer :: ios
128  Real, Allocatable :: x_kde_temp(:), f_kde_temp(:)
129  Real :: x_kde, f_kde
130  x = (/-2.1,-1.3,-0.4,1.9,5.1,6.2/)
131  numx = 0
132  numf = 0
133  Call Hist(n,x)
134  Call KDE(n,x)
135  Open(11,file='t_kde_x.txt')
136  Do
137    Read(11,*,iostat=ios)
138    numx = numx+1
139    If (ios/=0) Exit
140  Enddo
141  Close(11)
142  Open(12,file='t_kde_f.txt')
143  Do
144    Read(12,*,iostat=ios)
145    numf = numf+1
146    If (ios/=0) Exit
147  Enddo
148  Close(12)
149  !Write(*,*) numx, numf
150  If (numf/=numx) Then
151    Write(*,*) 'data pointa have error'
152    Else
153    num = numx-1
154  Endif
155  Allocate(x_kde_temp(num), f_kde_temp(num))
156  Open(13,file='t_kde_x.txt')
157  Do i=1,num
158    Read(13,*) ipi, ipj, x_kde_temp(i)
159  Enddo
160  Close(13)
161  Open(14,file='t_kde_f.txt')
162  Do i=1,num
163    Read(14,*) ipi, ipj, f_kde_temp(i)
164  Enddo
165  Close(14)
166  outer: Do i = 1, num-1
167    iptr = i
168    inner: Do j =i+1, num
169      minval: If(x_kde_temp(j) < x_kde_temp(iptr)) Then
170        iptr = j
171      Endif minval
172    Enddo inner
173    swap: If(i /= iptr) Then
174      temp = x_kde_temp(i)
175      x_kde_temp(i) = x_kde_temp(iptr)
176      x_kde_temp(iptr) = temp
177      tempf = f_kde_temp(i)
178      f_kde_temp(i) = f_kde_temp(iptr)
179      f_kde_temp(iptr) = tempf
180    Endif swap
181  Enddo outer  
182  Open(15,file='t_kde_xf.txt')
183  Do i=1,num
184    Write(15,*) i, x_kde_temp(i), f_kde_temp(i)
185  Enddo
186  Close(15)
187  Open(16,file='t_kde.txt')
188  ip = 0
189  Do i=1,num
190    If (i<ip.or.i==ip) Cycle
191    x_kde = x_kde_temp(i)
192    f_kde = f_kde_temp(i)
193    Do j=i+1,num
194      If (abs(x_kde-x_kde_temp(j))<0.0001) Then
195        x_kde = x_kde_temp(i)
196        f_kde = f_kde + f_kde_temp(j)
197        ip = j
198      Endif
199    Enddo
200    Write(16,*) x_kde, f_kde
201  Enddo
202  Close(16)
203  Deallocate(x_kde_temp, f_kde_temp)
204End Program main

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

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

规矩勋章