[未知用户] 是的。谢谢三楼。对于离散的分布列(在连续成为密度函数)的值是在[0,1](由于此时的分布列即是概率值);连续的密度函数值不一定在[0,1]之间,因为它只是函数值不是概率(三楼);例如三楼给的dnorm(0,0,000001)=398942.3;均匀分布[0,0.01]密度函数值全是100!
文章提供的方法很有实际意义。现实中的很多分布原本就有区间限制。
我现在就碰到一个这样的问题。根据住院费用的一定赔付数据求出住院赔付的核密度估计。这里的费用明显是非负的。进而需要求解赔付的期望和方差。这里需要用到积分。
但我遇到的问题是,因为要求根据病人男女不同,医院等级不同等,需要进一步求出每一个细分赔付的核密度估计,并求期望和方差。 根据文章中给出的核密度函数 den.est=function(u,ui,h) 。只需要求 den.est(u,ui,h)*u 在0-Inf上的积分 就可求解期望。 但是这样作用在三个数上的自定义函数,怎样运用tapply 和by 之类的函数实现快捷地对dat在其他条件限制下的核密度函数?一般运用tapply都只是求出了一个值,但是这里需要求出一个函数。如果每次人工将dat按条件挑出来,费时且当限制条件种类很多(如医院等级很多),容易出错。 不知道有没有什么好的方法?
边学R边用,对R了解不是很深。
6 个月 后
I like this boundary bias reduction trick, and it reminds me the somewhat remotely related generalized jackknife bias correction.
[未知用户] 你终于来主站了,小的恭候大驾已几年:)
I'm looking for more biostats/genomcis related topics here :D. One of the problems is to estimate the density of $p$-values on the unit interval, for which this trick could apply.

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
    8 年 后

    Cloud2016 @Ihavenothing 我觉得原标题是对的,只是 URL 中的词写错了。这篇文章讨论的是 x 取值有边界限制的情况,所以 URL 里的 unbounded 应该改成 bounded。如果改了 URL 记得在 _redirects 文件中做重定向。

      @邱怡轩 Cloud2016 大神们,我有一个问题想请教一下。主要还是bandwidth(下面简称h)的估计的问题。

      Wiki

      我看了一下wikipedia里面给出的近似的计算:

      h=(4σ^53n)151.06σ^n1/5h=(\frac{4\hat{\sigma}^5}{3n})^{\frac{1}{5}}\approx 1.06 \hat{\sigma} n^{-1/5}

      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是一致的)?

        6 个月 后
        3 年 后

        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)

          5 天 后

          yayalisa 你好,这块我也没搞太深,如果你需要,我可以把Silverman的<Density Estimation for statistics and data analysis>发给你。那里面有。