问题描述
来自 R文档 ...
在1到400(等于np.arange(1,401)
)的数组上,nrd0将返回31.39367.当我尝试在python中实现类似的东西时...
On an array from 1 to 400 (equivalent to np.arange(1,401)
), nrd0 will return 31.39367. When I try to implement something similar in python...
def nrd0_python(x):
X = min(np.std(x), np.percentile(x,25))
top = 0.9*X
bottom = 1.34*len(x)**(-0.2)
return top/bottom
我超过200(准确地说是224.28217762858455)
I get upwards of 200 (to be exact, 224.28217762858455)
Silverman的经验法则是否有已知的python函数?
Are there any known python functions for silverman's rule of thumb?
推荐答案
您对IQR的计算不正确:
You don't have the right calculation for IQR:
import numpy
def bw_nrd0(x):
if len(x) < 2:
raise(Exception("need at least 2 data points"))
hi = numpy.std(x, ddof=1)
q75, q25 = numpy.percentile(x, [75 ,25])
iqr = q75 - q25
lo = min(hi, iqr/1.34)
if not ((lo == hi) or (lo == abs(x[0])) or (lo == 1)):
lo = 1
return 0.9 * lo *len(x)**-0.2
这与R函数的返回结果相同,表示为:
This returns the same as the R function, which is given as:
> 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)
}
<bytecode: 0x0000000010c688b0>
<environment: namespace:stats>
我对R
代码中第二个if
语句的原始阅读不正确,请参见.它旨在捕获lo == 0
(由于sd(x) == 0
)和矢量全零的情况.如果我的理解是正确的,则应显示为:
if not lo:
if hi:
lo = hi
elif abs(x[0]):
lo = abs(x[0])
else:
lo = 1
Like the R
code, you can shorten this to :
lo = lo or hi or abs(x[0]) or 1
It is important to include these checks, along with the check for the length of the vector.
这篇关于R的nrd0是否有scipy/numpy替代项?的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持!