在R中计算总体标准差

Question

寻找在R中计算样本数量大于10的总体标准差的方法。目前无法提取R中的源代码以查找计算方法。

# 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

现在，人口标准差需要与 100 个观测值的样本标准差 sd(x) 相匹配。

# 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

- eyeOfTheStorm

Answer 1

请检查问题。 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)

- G. Grothendieck

---

感谢您的帮助。看起来可以在不依赖平均数/方差/标准差的情况下进行计算 - 请参见以下内容。 - eyeOfTheStorm

但是你为什么想要这样做呢？ - Ben Bolker

感谢提供RH2示例。 - eyeOfTheStorm

@BenBolker 对于使用“充分统计量” - eyeOfTheStorm

哦，我不确定那是什么意思。 - Ben Bolker

对于固定的n，总体标准差是样本标准差的标量倍数，因此它们都代表从足够统计观点来看相同的样本空间划分，因此是等效的。 - G. Grothendieck

Answer 2

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

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

- Bernhard

Answer 3

## 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 

- eyeOfTheStorm

Answer 4

我刚刚花了相当长的时间寻找一个带有计算总体标准差的函数包。以下是结果：

1) radiant.data::sdpop 应该是一个不错的函数 (查看文档)

2) multicon::popsd 也很好用，但需要查看文档以理解第二个参数的含义

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

4) rfml::sd.pop 将无法正常工作，除非使用ml.data.frame (MarkLogic Server)

希望这能帮到您。

- bik_92