有边界区间上的核密度估计
我现在就碰到一个这样的问题。根据住院费用的一定赔付数据求出住院赔付的核密度估计。这里的费用明显是非负的。进而需要求解赔付的期望和方差。这里需要用到积分。
但我遇到的问题是,因为要求根据病人男女不同,医院等级不同等,需要进一步求出每一个细分赔付的核密度估计,并求期望和方差。 根据文章中给出的核密度函数 den.est=function(u,ui,h) 。只需要求 den.est(u,ui,h)*u 在0-Inf上的积分 就可求解期望。 但是这样作用在三个数上的自定义函数,怎样运用tapply 和by 之类的函数实现快捷地对dat在其他条件限制下的核密度函数?一般运用tapply都只是求出了一个值,但是这里需要求出一个函数。如果每次人工将dat按条件挑出来,费时且当限制条件种类很多(如医院等级很多),容易出错。 不知道有没有什么好的方法?
边学R边用,对R了解不是很深。
BTW: The ULR for this article says "unbounded" region instead of "bounded" region: http://cos.name/2010/04/kernel-density-estimation-with-unbounded-region/#comment-2564
Cloud2016 @Ihavenothing 我觉得原标题是对的,只是 URL 中的词写错了。这篇文章讨论的是 x 取值有边界限制的情况,所以 URL 里的 unbounded 应该改成 bounded。如果改了 URL 记得在 _redirects 文件中做重定向。
- 已编辑
@邱怡轩 Cloud2016 大神们,我有一个问题想请教一下。主要还是bandwidth(下面简称h
)的估计的问题。
Wiki
我看了一下wikipedia里面给出的近似的计算:
R
还有我查看了R中的bw.nrd0
function (x)
{
if (length(x) < 2L)
stop("need at least 2 data points")
hi <- sd(x)
if (!(lo <- min(hi, IQR(x)/1.34)))
(lo <- hi) || (lo <- abs(x[1L])) || (lo <- 1)
0.9 * lo * length(x)^(-0.2)
}
bw.nrd
function (x)
{
if (length(x) < 2L)
stop("need at least 2 data points")
r <- quantile(x, c(0.25, 0.75))
h <- (r[2L] - r[1L])/1.34
1.06 * min(sqrt(var(x)), h) * length(x)^(-1/5)
}
结果对比
数据:原贴中的数据。
JMP(因为公司限制,不好贴图)计算结果是:0.22706547075559
R中:
- bw.nrd: 0.2225342
- bw.nrd0: 0.1889441
我比较疑惑,到底应该怎么算以及JMP中使用的是什么方法(有时候模拟的数据JMP和bw.nrd0是一致的)?
Jonie_Y 大神,主要也是bandwidth(下面简称h)的问题,我用的也是Wiki的近似计算公式,但我用的是python编写代码,意思应该也是差不多的,我的x,y数据指的是经纬度,能否帮我看看代码中的公式是对的吗?因为我的数据是二维的,您使用的是一维数据。还有r <- quantile(x, c(0.25, 0.75)) h <- (r[2L] - r[1L])/1.34这两个行代码这样写的原理是什么呢?希望能得到您的解答。
x_mean=np.mean(x)
y_mean=np.mean(y)
n = len(x)
for i in range(len(x)):
dis.append(math.pow((math.pow(x-x_mean,2))+(math.pow(y-y_mean,2)),0.5))
dis_array = np.array(dis)
q1 = np.percentile(dis_array, 25)
q3 = np.percentile(dis_array, 75)
dis_std=np.std(dis_array)
delta=min(dis_std,(q3-q1)/1.34)
h=1.06deltapow(len(x),-0.2)