周期性高斯的实现答案

【问题标题】：implementation of periodic gaussian周期性高斯的实现
【发布时间】：2014-09-24 16:57:44
【问题描述】：

我正在尝试在 C、MATLAB 或 Python 中实现周期性高斯。

评估如下定义的周期性高斯函数的正确方法是什么

我目前正在根据下面的公式进行评估，以避免从负到正无穷的总和：

提前致谢。

【问题讨论】：

这个问题是关于特定编程语言或工具的吗？
不，它不是任何语言特有的，我问的是如何以编程方式在任何语言中实现周期性高斯，困难在于如何解决从负无穷到正无穷的求和。如果您可以在 C、Matlab 或 Python 中提供答案，将对我有所帮助。
@user3657953 你最好问你的问题here。
我不明白在使用模运算符时如何避免无穷大的总和。 k的分辨率是多少？ k 是整数吗？实数？如果它是一个实数，那么除非您使用符号数学，例如 Python 中的 Theano 或 MATLAB 中的 Symbolic Math Toolbox，否则您无法做到这一点。
是的，k是整数

标签： python c matlab numeric numerical-methods

【解决方案1】：

好吧，你不应该计算无限和，因为一旦你达到 (x-kL) >> 2sigma，你就会达到浮点精度的极限。

因此，您应该能够从找到 x - kL 的最小值开始（即，只需设置 x = x mod L 和 k=0，这样做是合法的，因为这是一个无限和），然后在 k 处添加项= +/- 1, +/- 2, ... 直到达到浮点限制。下面是一些说明这个想法的示例 MATLAB 代码 - 我只是把它整理出来，所以我不能保证它没有错误，但它似乎表现出一些基本的预期行为。

    function [result] = Periodic_Gaussian(x, L, sigma)

    gaussian = @(y) 1/(2*pi*sigma)*exp(-y.^2 ./ 4 ./sigma^2);

    x = mod(x, L);

    oldresult = NaN;
    newresult = gaussian(x);
    k = 1;
    while any(newresult ~= oldresult)
        oldresult = newresult;
        newresult = oldresult + gaussian(x-k*L) + gaussian(x+k*L);
        k = k+1;
    end
    result = newresult;

希望对您有所帮助！

编辑：在指数参数的分母中遗漏了 4 倍，如果需要，更新代码以采用 x 向量。

【讨论】：

谢谢，但是这个函数对于任何给定的 L 都返回相同的结果
我不知道你测试了哪些 L 值，但它不应该，除非我在复制它时出错：例如，Periodic_Gaussian(11, 2, 1) = 0.2821, Periodic_Gaussian(11, 8, 1) =0.0171。
我的意思是固定 L 和变化 x 和 sigma 的结果相同，在您的示例中，您保持 x 恒定（= 11）顺便说一句，这就是我得到的 Periodic_Gaussian(11, 8, 1) = 1.9641e-05 和 PeriodicGaussian(11, 2, 1) = 0.1171.... 对于 L = 1，我得到 PeriodicGaussian(11, 1, 5) = PeriodicGaussian(0, 1, 5) = 0.2821
你的第二组输出对我来说很有意义 - 高斯是周期性的。如果 L = 1，那么函数应该返回相同的值 11 或 0，因为 11 是 0 + 周期 L 的整数倍。PeriodicGaussian(x1, L, sigma) == PeriodicGaussian(x2, L, sigma) 对于任意 x2 = x1+k*L，这是周期性的定义.另请注意，当 sigma >> L 时，如您的示例所示，此函数将接近所有 x 的常数值。虽然我确实注意到上面的代码中有一个错误解释了我们得到的不同值，但我现在已经对其进行了编辑以修复。
是的，你是对的..你的代码中有一个小错字，根据定义它是 2 而不是 4。最后，我有一个代码可以在没有下面发布的 while 循环的情况下执行相同的操作。谢谢！

【解决方案2】：

function [result] = PeriodicGaussian(x, L, sigma)                                                                                                                                                                                             



    gaussian = @(y, sigma) 1/(2*pi*sigma)*exp(-y.^2 ./ 2 ./sigma^2);                                                                                                                                                                           


    x0 = mod(x, L)                                                                                                                                                                                                                             
    x1 = mod(x, -1 * L)                                                                                                                                                                                                                        
    result = gaussian(x0, sigma) + gaussian(x1, sigma);                                                                                                                                                                                        
    correctionIdx = (x0 == 0 & x1 == 0);                                                                                                                                                                                                       
    result(correctionIdx) = 0.5 * result(correctionIdx);                                                                                                                                                                                       

end

【讨论】：