2

我有一个看起来像这样的向量：

> y
 [1] 6.7 5.3 3.3 6.7 3.3 4.7 4.7 6.7 3.3 6.7

我试图h = 3使用以下公式计算时移时的估计自协方差：

在此处输入图像描述

这就是我所做的，我想知道有没有比硬编码更简单的方法：

> 1/10 * (
+     (y[4] - mean(y))*(y[1] - mean(y)) + 
+     (y[5] - mean(y))*(y[2] - mean(y)) + 
+     (y[6] - mean(y))*(y[3] - mean(y)) + 
+     (y[7] - mean(y))*(y[4] - mean(y)) + 
+     (y[8] - mean(y))*(y[5] - mean(y)) + 
+     (y[9] - mean(y))*(y[6] - mean(y)) + 
+     (y[10] - mean(y))*(y[7] - mean(y))
+ )
[1] -0.04848
4

4 回答 4

4

您可以像这样定义自己的函数：

my.cov <- function(y,h){
  y.bar <- mean(y)
  T <- length(y)
  (1/T) * sum((y[(h+1):T] - y.bar) * (y[1:(T-h)] -y.bar))
}
#Checking the function
> y <- c(6.7, 5.3, 3.3, 6.7, 3.3, 4.7, 4.7, 6.7, 3.3, 6.7)
> my.cov(y, h=3)
[1] -0.04848

请注意，由于是矢量化运算符，因此不for涉及循环。*

于 2013-10-26T20:01:02.313 回答
3

我会做这个矢量化而不是使用循环。更快，更清洁：

> y <- c(6.7,5.3,3.3,6.7,3.3,4.7,4.7,6.7,3.3,6.7)
> t <- 1:7
> (1/10)*sum((y[t+3]-mean(y))*(y[t]-mean(y)))
[1] -0.04848

把它变成一个函数是微不足道的......

foo <- function(x,lag) {
    n <- length(x)
    i <- 1:(n-lag)
    (1/n) * sum( (y[ i + lag ] - mean(y) ) *( y[ i ] - mean(y) ) )
    }

foo( y , 3 )
# [1] -0.04848
于 2013-10-26T20:00:34.880 回答
2

这是一种功能化和矢量化的方法：

autocov <- function(x, lag = 3)
  sum((tail(x, -lag) - mean(x)) *
      (head(x, -lag) - mean(x))) / length(x)

x <- c(6.7, 5.3, 3.3, 6.7, 3.3, 4.7, 4.7, 6.7, 3.3, 6.7)
autocov(x)
# [1] -0.04848
autocov(c(1:3, 1:3))
# [1] 0.33333333

最后一个示例是向您展示我如何怀疑您的自协方差公式是错误的。我本来想：

autocov <- function(x, lag = 3) cov(tail(x, -lag), tail(x, -lag))

autocov(c(1:3, 1:3))
# [1] 1

（即，每个信号都应该有自己的平均值，并且使用您的示例，缩放比例应该是 1/7。）

于 2013-10-26T20:00:18.387 回答
0

为什么不使用包acf中的自协方差函数stats？

顺便说一句：当您处理时间序列时，我建议使用时间序列对象而不是普通向量。我更喜欢包zoo。

library(zoo)
test_zoo=zoo(c(6.7,5.3,3.3,6.7,3.3,4.7,4.7,6.7,3.3,6.7),1:10)
acf(test_zoo,type="covariance",lag.max=3,plot=FALSE)[3]

#Autocovariances of series ‘test_zoo’, by lag
#      3 
#-0.0485
于 2013-10-27T00:46:26.803 回答