帕斯卡高斯答案

【问题标题】：Gaussian in pascal帕斯卡高斯
【发布时间】：2016-07-14 22:48:57
【问题描述】：

我尝试将代码直接从 java 源代码移植到 pascal，但是它会引发运行时错误。

如何获得合适的高斯曲线？帕斯卡内置函数呢？

原始源代码：

    synchronized public double nextGaussian() {
    // See Knuth, ACP, Section 3.4.1 Algorithm C.
    if (haveNextNextGaussian) {
        haveNextNextGaussian = false;
        return nextNextGaussian;
    } else {
        double v1, v2, s;
        do {
            v1 = 2 * nextDouble() - 1; // between -1 and 1
            v2 = 2 * nextDouble() - 1; // between -1 and 1
            s = v1 * v1 + v2 * v2;
        } while (s >= 1 || s == 0);
        double multiplier = StrictMath.sqrt(-2 * StrictMath.log(s)/s);
        nextNextGaussian = v2 * multiplier;
        haveNextNextGaussian = true;
        return v1 * multiplier;
    }
}

第一次尝试帕斯卡端口（抛出运行时错误）：

  function log (n : double) : double; 
  begin 
    result := ln(n) / ln(10); 
  end; 

  var hgauss : boolean;
  var ngauss : double;

  function gauss() : double;
  var x1, x2, w : double;
  begin
    if hgauss then
    begin
      result := ngauss;
      hgauss := false;
    end else
    begin
      repeat
        x1 := 2.0 * rand() - 1.0;
        x2 := 2.0 * rand() - 1.0;
        w := x1 * x1 + x2 * x2;
      until w >= 1.0;         

      w := sqrt( (-2.0 * log( w ) ) / w );
      result := x1 * w;
      ngauss := x2 * w;
      hgauss := true;
    end;
  end;

这里的浮点运算无效：

w := sqrt((-2.0 * log( w ) ) / w);

第二次尝试转换（运行但我不确定数学是否正确）：

  function log (n : double) : double; 
  begin 
    result := ln(n) / ln(10); 
  end; 

  var hgauss : boolean;
  var ngauss : double;

  function gauss() : double;
  var x1, x2, w, num : double;
  begin
    if hgauss then
    begin
      result := ngauss;
      hgauss := false;
    end else
    begin
      repeat
        x1 := 2.0 * rand() - 1.0;
        x2 := 2.0 * rand() - 1.0;
        w := x1 * x1 + x2 * x2;
      until w >= 1.0;         

      num := -2.0 * log( w )  / w;
      w := sqrt(abs(num));
      if num < 0 then w := -w;
      result := x1 * w;
      ngauss := x2 * w;
      hgauss := true;
    end;
  end;

【问题讨论】：

您使用的是哪个 pascal 编译器？ java句柄如何除以0？ .你需要在你的帕斯卡中处理它。您也忽略了重复循环中 w = 0 的情况，（Java 中的 s == 0）
我正在使用 SCAR DIVI 3.41。据我所知，Java 不允许除以零，所以我认为编译器应该以同样的方式处理它。还有一个名为 RandG() 的函数，据说它返回一个“高斯”，但它似乎不在正态分布中。
如果你用 Free Pascal 编码，你也可以使用randg 函数来获取高斯随机数。有关如何从大量不同分布中生成随机数的更多信息，请参阅wiki.freepascal.org/Generating_Random_Numbers。

标签： java logging random pascal gaussian

【解决方案1】：

你从 Java 到 Pascal 的移植有一个重要部分有问题

do {...} while (s >= 1 || s == 0);

应该翻译成

repeat {...} until ((s<1) and (s<>0));

所以你有错误的终止条件。如果0 < s < 1，Java 将终止循环，但如果w >= 1，则您的循环结束。

如果w > 1 你有-2*ln(w) < 0 并且浮点异常来自取负数的平方根！

对于大多数 Pascal 版本，您对标准函数的命名是不寻常的，国际海事组织它应该读

repeat
  x1 := 2.0 * random - 1.0;
  x2 := 2.0 * random - 1.0;
  w := x1 * x1 + x2 * x2;
until (w<1.0) and (w>0.0);         

w := sqrt(-2.0*ln(w)/w);
result := x1 * w;
ngauss := x2 * w;

请注意，您确实必须使用ln 而不是您自制的以10 为底的对数log。使用的方法是Marsaglia's polar method.

【讨论】：

确实：java 的 Math.log 有自然底 e，所以 ln 已经是。

【解决方案2】：

rand() 在[0,1) 范围内（0 <= rand() < 1）
所以2.0 * rand() - 1.0 在[-1,1) 范围内
所以x1 和x2 在[-1,1) 范围内
所以w := x1 * x1 + x2 * x2 在[0,2) 范围内

在sqrt( -2.0 * ln( w ) / w ) 中，w 是肯定的
所以自然对数：ln(w) 应该是负数
所以 w 应该在(0,1) 范围内
所以循环不应该退出until (w > 0.0)and (w < 1.0);

工作示例代码（使用 SCAR Divi 3.41.00）：

program New; 

 var hgauss : boolean;
 var ngauss : double;

  function gauss() : double;
  var x1, x2, w : double;
  begin
    if hgauss then
    begin
      result := ngauss;
      hgauss := false;
    end else
    begin
      repeat
        x1 := 2.0 * rand() - 1.0;
        x2 := 2.0 * rand() - 1.0;
        w := x1 * x1 + x2 * x2;
      until (w > 0.0)and (w < 1.0);   
      w := sqrt( -2.0 * ln( w ) / w );
      result := x1 * w;
      ngauss := x2 * w;
      hgauss := true;
    end;
  end;

begin
  writeln( gauss() );
  writeln( gauss() );
end.

【讨论】：