计算R中的人口标准偏差

Question

寻找一种方法来计算 R 中的总体标准偏差 - 使用超过 10 个样本。无法在R中提取C源代码来查找计算方法。

# Sample Standard Deviation 
# Note: All the below match with 10 or less samples
n <- 10 # 10 or greater it shifts calculation
set.seed(1)
x <- rnorm(n, 10)

# Sample Standard Deviation
sd(x)
# [1] 0.780586
sqrt(sum((x - mean(x))^2)/(n - 1))
# [1] 0.780586
sqrt(sum(x^2 - 2*mean(x)*x + mean(x)^2)/(n - 1)) # # Would like the Population Standard Deviation equivalent using this.
# [1] 0.780586
sqrt( (n/(n-1)) * ( ( (sum(x^2)/(n)) ) - (sum(x)/n) ^2 ) )
# [1] 0.780586

现在，人口标准偏差需要匹配 sd(x) 和 100 个计数。

# Population Standard Deviation 
n <- 100 
set.seed(1)
x <- rnorm(x, 10)

sd(x)
# [1] 0.780586

sqrt(sum((x - mean(x))^2)/(n))
# [1] 0.2341758

sqrt(sum(x^2 - 2*mean(x)*x + mean(x)^2)/(n)) 
# [1] 0.2341758

# Got this to work above using (eventual goal, to fix the below):
# https://en.wikipedia.org/wiki/Algebraic_formula_for_the_variance
sqrt( (n/(n-1)) * ( ( (sum(x^2)/(n)) ) - (sum(x)/n) ^2 ) )  # Would like the Population Standard Deviation equivalent using this.
# [1] 3.064027

Answer 1

我认为最简单的方法是从sd 快速定义它：

sd.p=function(x){sd(x)*sqrt((length(x)-1)/length(x))}

Answer 2

我刚刚花了相当多的时间寻找一个包含人口标准偏差的现成函数的包。结果如下：

1) radiant.data::sdpop 应该是个不错的函数（见documentation）

2) multicon::popsd 也很有效，但请检查 documentation 以了解第二个参数是什么

3) muStat::stdev 与 unbiased=FALSE 不能正常工作。在 github 页面上，似乎在 2012 年有人将其设置为 sd(x)*(1-1/length(x)) 而不是 sd(x)*sqrt(1-1/length(x))...

4) 如果没有 ml.data.frame (MarkLogic Server)，rfml::sd.pop 将无法工作

我希望这会有所帮助。

Answer 3

请检查问题。 rnorm 的第一个参数应该是 n。

总体和样本标准差为：

sqrt((n-1)/n) * sd(x) # pop
## [1] 0.8936971

sd(x) # sample
## [1] 0.8981994

它们也可以这样计算：

library(sqldf)
library(RH2)

sqldf("select stddev_pop(x), stddev_samp(x) from X")
##   STDDEV_POP("x") STDDEV_SAMP("x")
## 1       0.8936971        0.8981994

注意：我们使用了这个测试数据：

set.seed(1)
n <- 100
x <- rnorm(n)
X <- data.frame(x)

Answer 4

## Sample Standard Deviation  
n <- 10 # Sample count
set.seed(1)
x <- rnorm(n, 10)

sd(x) # Correct 
# [1] 0.780586
sqrt(sum((x - mean(x))^2)/(n - 1)) # Correct 
# [1] 0.780586
sqrt(sum(x^2 - 2*mean(x)*x + mean(x)^2)/(n - 1)) # Correct 
# [1] 0.780586
sqrt( (n/(n-1)) * ( ( (sum(x^2)/(n)) ) - (sum(x)/n) ^2 ) ) # Correct 
# [1] 0.780586
sqrt((sum(x^2) - (sum(x)^2/n))/(n-1)) # Correct 
# [1] 0.780586
sqrt( (n/(n - 1)) * ( (sum(x^2)/(n))  - (sum(x)/n) ^2 ) ) # Correct 
# [1] 0.780586


## Population Standard Deviation  
n <- 100 # Note: 10 or greater biases var() and sd()
set.seed(1)
x <- rnorm(n, 10)

sd(x) # Incorrect Population Standard Deviation!!
# [1] 0.8981994
sqrt(sum((x - mean(x))^2)/(n)) # Correct
# [1] 0.8936971
sqrt(sum(x^2 - 2*mean(x)*x + mean(x)^2)/(n)) # Correct 
# [1] 0.8936971
sqrt((sum(x^2) - (sum(x)^2/n))/(n)) # Correct
# [1] 0.8936971
sqrt( (n/(n)) * ( (sum(x^2)/(n))  - (sum(x)/n) ^2 ) ) # Correct 
# [1] 0.8936971