我尝试将代码直接从 java 源代码移植到 pascal,但是它引发了运行时错误。
如何获得正确的高斯曲线? 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;
}
}
第一次尝试 pascal 端口(引发运行时错误):
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;
最佳答案
您从 Java 到 Pascal 的移植有一个重要部分出现错误
do {...} while (s >= 1 || s == 0);
应翻译为
repeat {...} until ((s<1) and (s<>0));
所以你的终止条件是错误的。 Java 终止循环 if 0 < s < 1
,但是如果 w >= 1
则循环结束.
如果w > 1
你有-2*ln(w) < 0
浮点异常来自于取负数的平方根!
对于大多数 Pascal 版本,标准函数的命名是不寻常的, IMO 它应该阅读
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 - 以帕斯卡为单位的高斯,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/38364880/