gpt4 book ai didi

algorithm - 核密度估计 julia

转载 作者:塔克拉玛干 更新时间:2023-11-03 02:52:56 28 4
gpt4 key购买 nike

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

算法如下:

\$ f(x) = \frac{1}{n*h} * \sum_{i = 1}^n K(\frac{x - X_i}{h}) \$

在哪里

\$ K(u) = 0.5*I(|u| <= 1)\$ with \$ u = \frac{x - X_i}{h}\$

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

这是我的实现:

#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的公式,因此我真的不明白错误在哪里。

最佳答案

我不是 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 内核 :)

关于algorithm - 核密度估计 julia,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/32394393/

28 4 0
Copyright 2021 - 2024 cfsdn All Rights Reserved 蜀ICP备2022000587号
广告合作:1813099741@qq.com 6ren.com