Клен: ГСЧ не случайно

Question

Я "находил Пи" методом Монте-Карло, но ответ был неверным. Оригинальный код был:

RandomTools[MersenneTwister]: with(Statistics):

tries := 10000:

s := 0; 
for i to tries do 
    if GenerateFloat()^2+GenerateFloat()^2 < 1 then s := s+1 end if; 
end do: 
evalf(4*s/tries)

Это дает ответ aroud 2.8-2.85

когда я изменяю код на

s := 0; 
x := Array([seq(GenerateFloat(), i = 1 .. tries)]); 
y := Array([seq(GenerateFloat(), i = 1 .. tries)]); 
for i to tries do 
if x[i]^2+y[i]^2 < 1 then s := s+1 end if;
end do:
evalf(4*s/tries)

Тогда ответ правильный. Я понятия не имею, почему я не могу сгенерировать число в цикле "for".

Я обнаружил, что среднее значение одно и то же, но разница отличается. За:

tries := 100000; 
A := Array([seq(GenerateFloat(), i = 1 .. 2*tries)]); 
s1 := Array([seq(A[i]^2+A[tries+i]^2, i = 1 .. tries)]); 
Mean(s1); 
Variance(s1);
s2 := Array([seq(GenerateFloat()^2+GenerateFloat()^2, i = 1 .. tries)]);
Mean(s2); 
Variance(s2);

вывод:

0.6702112097021581
0.17845439723457215
0.664707674135025
0.35463131700965245

Что с этим не так? GenerateFloat() должен быть максимально равномерным.

Answer 1

Автоматическое упрощение превращается,

GenerateFloat()^2+GenerateFloat()^2

в,

2*GenerateFloat()^2

до GenerateFloat() оценивается.

Одно простое изменение, чтобы заставить его работать так, как вы ожидали, будет разделять их. Например,

restart:
with(RandomTools[MersenneTwister]):
tries := 10^4:
s := 0:
for i to tries do
  t1,t2 := GenerateFloat(),GenerateFloat();
  if t1^2+t2^2 < 1 then s := s+1 end if;
end do:
evalf(4*s/tries);

Другой способ - использовать немного другую конструкцию, которая не упрощается автоматически. Учтите, что одинарные правильные кавычки (неравные кавычки) не останавливают автоматическое упрощение (которое, если хотите, является определением термина).

'f()^2 + f()^2';                                                
                                     2
                                2 f()

Но следующее автоматически не упрощается,

a:=1:
'f()^2 + a*f()^2';
                                2        2
                             f()  + a f()

Поэтому еще один простой обходной путь,

restart:
with(RandomTools[MersenneTwister]):
tries := 10^4:
s := 0:
a := 1;
for i to tries do
  if GenerateFloat()^2 + a*GenerateFloat()^2 < 1 then s := s+1 end if;
end do:
evalf(4*s/tries);