众所周知,计算相关性非常的简单,因为R
语言中有函数cor.test()
,该函数可以计算多种方法的相关性检验,返回相关性,Pvalue等检验值,但是这个函数在Julia
中并不存在,让Julia作为一门科学计算语言显得并不完美。
首先,我们来看一下R语言计算线性相关的结果是什么样子的?
a<-c(1,3,4,5,7)
b<-c(2,4,6,8,9)
c<-cor.test(a,b)
c
Pearson's product-moment correlation
data: a and b
t = 7.7771, df = 3, p-value = 0.004423
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.6757774 0.9984873
sample estimates:
cor
0.976086
首先,我们生成了两个示例的向量,这两个向量必须等长
然后,我们使用了cor.test
,因为默认就是皮尔森,因此我们不再使用method参数,得到的结果也很丰富。
但是,我们比较Julia呢?
首先,julia标准库并没有计算相关性的函数,这里我们可以调用Statistics包来计算相关性
a = [1,3,4,5,7]
b = [2,4,6,8,9]
using Statistics
c = cor(a,b)
0.9760860118037878
然后,我们参考R 语言的定义计算p-value和95%的置信区间。
using Distributions ##载入包,用来计算T统计和t分布的cdf值
using Statistics ##计算cor
## 根据r和样本量n计算置信区间
function interval(r::Float64, n::Int64)
z = 0.5 * log(exp(1), (1 + r) / (1 - r))
s = 1 / (sqrt(n - 3))
zl = z - 1.96 * s
zh = z + 1.96 * s
low = (exp(2 * zl) - 1) / (exp(2 * zl) + 1)
high = (exp(2 * zh) - 1) / (exp(2 * zh) + 1)
(low, high)
end
## 输入两个向量,返回各个参数
function cor_test(x::AbstractVector, y::AbstractVector)
if length(x) == length(y)
r = cor(x, y) #相关性值
n = length(x)
df = n - 2 ## 自由度
t = r * sqrt((df) / (1 - r^2)) ##T统计量
tdist = cdf(TDist(df), t) ##t分布的累计
if r > 0
pvalue = (1 - tdist) * 2
else
pvalue = (tdist) * 2
end
else
println("请输入等长的两个向量")
end
fin = (r = r, pvalue = pvalue, df = df, t = t, cfi = interval(r, n))
end
测试以下函数的返回是否和R语言一致
cor_test(a,b)
(r = 0.9760860118037878, pvalue = 0.004423314933107214, df = 3, t = 7.777137710478182, cfi = (0.6757635765936858, 0.9984873283114327))
c = cor_test(a,b)
c.r
0.9760860118037878
c.pvalue
0.004423314933107214