核密度估计 julia

Question

我正在尝试实现内核密度估计。但是，我的代码没有提供应有的答案。它也是用 julia 编写的，但代码应该是不言自明的。

算法如下：

在哪里

因此，该算法测试 x 与由某个常数因子（binwidth）加权的观测值 X_i 之间的距离是否小于 1。如果是这样，它会将 0.5 / (n * h) 分配给该值，其中 n = #of 观察值。

这是我的实现：

#Kernel density function.
#Purpose: estimate the probability density function (pdf)
#of given observations
#@param data: observations for which the pdf should be estimated
#@return: returns an array with the estimated densities 

function kernelDensity(data)
|   
|   #Uniform kernel function. 
|   #@param x: Current x value
|   #@param X_i: x value of observation i
|   #@param width: binwidth
|   #@return: Returns 1 if the absolute distance from
|   #x(current) to x(observation) weighted by the binwidth
|   #is less then 1. Else it returns 0.
|  
|   function uniformKernel(x, observation, width)
|   |   u = ( x - observation ) / width
|   |   abs ( u ) <= 1 ? 1 : 0
|   end
|
|   #number of observations in the data set 
|   n = length(data)
|
|   #binwidth (set arbitraily to 0.1
|   h = 0.1 
|   
|   #vector that stored the pdf
|   res = zeros( Real, n )
|   
|   #counter variable for the loop 
|   counter = 0
|
|   #lower and upper limit of the x axis
|   start = floor(minimum(data))
|   stop = ceil (maximum(data))
|
|   #main loop
|   #@linspace: divides the space from start to stop in n
|   #equally spaced intervalls
|   for x in linspace(start, stop, n) 
|   |   counter += 1
|   |   for observation in data
|   |   |
|   |   |   #count all observations for which the kernel
|   |   |   #returns 1 and mult by 0.5 because the
|   |   |   #kernel computed the absolute difference which can be
|   |   |   #either positive or negative
|   |   |   res[counter] += 0.5 * uniformKernel(x, observation, h)
|   |   end
|   |   #devide by n times h
|   |   res[counter] /= n * h
|   end
|   #return results
|   res
end
#run function
#@rand: generates 10 uniform random numbers between 0 and 1
kernelDensity(rand(10))

这是被退回的：

> 0.0
> 1.5
> 2.5
> 1.0
> 1.5
> 1.0
> 0.0
> 0.5
> 0.5
> 0.0

其和为：8.5（累积分布函数，应为1。）

所以有两个bug：

1. 值未正确缩放。每个数字应约为其当前值的十分之一。事实上，如果观察次数增加 10^n n = 1, 2, ... 那么 cdf 也会增加 10^n

例如：

> kernelDensity(rand(1000))
> 953.53 

1. 它们的总和不等于 10（如果不是因为缩放错误，则不等于 1）。随着样本量的增加，错误变得更加明显：大约有。 5% 的观察结果未包括在内。

我相信我实现了公式 1:1，因此我真的不明白错误在哪里。

Answer 1

我不是 KDE 方面的专家，所以对所有这些都持保留态度，但您的代码的一个非常相似（但要快得多！）的实现将是：

function kernelDensity{T<:AbstractFloat}(data::Vector{T}, h::T)
  res = similar(data)
  lb = minimum(data); ub = maximum(data)
  for (i,x) in enumerate(linspace(lb, ub, size(data,1)))
    for obs in data
      res[i] += abs((obs-x)/h) <= 1. ? 0.5 : 0.
    end
    res[i] /= (n*h)
 end
 sum(res)
end

如果我没记错的话，密度估计应该积分为 1，也就是说，我们希望 kernelDensity(rand(100), 0.1)/100 至少接近 1。在上面的实现中，我到达那里，给或取 5%，但话又说回来，我们不知道 0.1 是最佳带宽（使用 h=0.135 而不是我到达 0.1% 以内），并且已知统一的内核只有大约 93% 的“效率”。

无论如何，Julia 中有一个非常好的 Kernel Density 包可用 here，所以您可能应该只使用 Pkg.add("KernelDensity") 而不是尝试编写自己的 Epanechnikov 内核 :)

Answer 2

指出错误：您有 n 个大小为 2h 的箱 B_i 覆盖 [0,1]，随机点 X 落在预期数量的 箱中。你除以 2 n h。

对于 n 个点，您的函数的期望值为。

实际上，您有一些大小

编辑：顺便说一句，如果您假设 bin 在 [0,1] 中具有随机位置，则该偏差很容易计算。那么这些垃圾箱平均会丢失其大小的 h/2 = 5%。