结果不匹配。浮点数错误？

Question

我正在尝试将计算自相关的 R function acf 重写为 C#：

class AC
{
    static void Main(string[] args)
    {
        double[] y = new double[] { 772.9, 909.4, 1080.3, 1276.2, 1380.6, 1354.8, 1096.9, 1066.7, 1108.7, 1109, 1203.7, 1328.2, 1380, 1435.3, 1416.2, 1494.9, 1525.6, 1551.1, 1539.2, 1629.1, 1665.3, 1708.7, 1799.4, 1873.3, 1973.3, 2087.6, 2208.3, 2271.4, 2365.6, 2423.3, 2416.2, 2484.8, 2608.5, 2744.1, 2729.3, 2695, 2826.7, 2958.6, 3115.2, 3192.4, 3187.1, 3248.8, 3166, 3279.1, 3489.9, 3585.2, 3676.5 };
        Console.WriteLine(String.Join("\n", acf(y, 17)));
        Console.Read();
    }

    public static double[] acf(double[] series, int maxlag)
    {
        List<double> acf_values = new List<double>();
        float flen = (float)series.Length;
        float xbar = ((float)series.Sum()) / flen;
        int N = series.Length;

        double variance = 0.0;
        for (int j = 0; j < N; j++)
        {
            variance += (series[j] - xbar)*(series[j] - xbar);
        }
        variance = variance / N;

        for (int lag = 0; lag < maxlag + 1; lag++)
        {
            if (lag == 0)
            {
                acf_values.Add(1.0);
                continue;
            }

            double autocv = 0.0;
            for (int k = 0; k < N - lag; k++)
            {
                autocv += (series[k] - xbar) * (series[lag + k] - xbar);
            }
            autocv = autocv / (N - lag);

            acf_values.Add(autocv / variance);
        }
        return acf_values.ToArray();
    }
}

这段代码有两个问题：

1. 对于大型数组（长度 = 25000），此代码大约需要 1-2 秒，而 R 的 acf 函数在不到 200 毫秒内返回。
2. 输出与 R 的输出不完全匹配。

关于我在哪里搞砸或对代码进行任何优化的任何建议？

        C#          R
    1   1           1
    2   0.945805846 0.925682317
    3   0.89060465  0.85270658
    4   0.840762283 0.787096604
    5   0.806487301 0.737850083
    6   0.780259665 0.697253317
    7   0.7433111   0.648420319
    8   0.690344341 0.587527097
    9   0.625632533 0.519141887
    10  0.556860982 0.450228026
    11  0.488922355 0.38489632
    12  0.425406196 0.325843042
    13  0.367735169 0.273845337
    14  0.299647764 0.216766466
    15  0.22344712  0.156888402
    16  0.14575994  0.099240809
    17  0.072389526 0.047746281
    18  -0.003238526    -0.002067146

Answer 1

您可以尝试更改此行：

autocv = autocv / (N - lag);

到这里：

autocv = autocv / N;

其中任何一个都是期望值的可接受除数，R 显然使用的是第二个。

要在不使用 C# 编译器的情况下查看这一点，我们可以阅读您拥有的表格，并通过将 C# 列中的每个值除以 N/(N - lag) 来调整值，并查看它们是否同意来自 R 的值。

N在这里是47，滞后范围是0到17，所以N - lag是47:30。

将上面的表格复制到我的本地剪贴板后：

cr <- read.table(file='clipboard', comment='', check.names=FALSE)
cr$adj <- cr[[1]]/47*(47:30)
max(abs(cr$R - cr$adj))
## [1] 2.2766e-09

更接近的近似值。

如果将flen 和xbar 定义为double 类型，您可能会做得更好，因为浮点数没有9 位小数精度。

R 速度如此之快的原因是 acf 是作为本机和非托管代码（C 或 FORTRAN）实现的。