使用 Monte Carlo 模拟限制计算 π

Question

我问过一个与此非常相似的问题所以我会在最后提到以前的解决方案，我有一个 website 可以用客户端的 CPU 计算 π 并将其存储在服务器上，到目前为止我有：

'701.766.448.388'点在圆圈内，总共'893.547.800.000'，这些数字是使用此代码计算的。（工作示例位于：https://jsfiddle.net/d47zwvh5/2/）

let inside = 0;
let size = 500;

for (let i = 0; i < iterations; i++) {
  var Xpos = Math.random() * size;
  var Ypos = Math.random() * size;

  var dist = Math.hypot(Xpos - size / 2, Ypos - size / 2);

  if (dist < size / 2) {
    inside++;
  }
}

问题

(4 * 701.766.448.388) / 893.547.800.000 = 3,141483638

这是我们得到的结果，直到第四位都是正确的，4应该是5。

以前的问题：

我搞砸了距离计算。
我从 0...499 开始放置圆，应该是 0...500
我没有使用浮动，这减少了 'resolution'

免责声明

可能只是我达到了极限，但是这个演示使用了100万个点，得到了3.16。考虑到我有大约 9000 亿，我认为它可能更准确。

我明白如果我想计算 π 这不是正确的方法，但我只是想确保一切正确所以我希望任何人都能发现错误或者我只需要更多 'dots'.

编辑： 有很多人提到这些数字是多么不切实际，但这些都是正确的，我现在已经将它们更新为正确的。

Answer 1

任何 FPU 操作都会降低您的准确性。为什么不做这样的事情：

let inside = 0;
for (let i = 0; i < iterations; i++)
  {
  var X = Math.random();
  var Y = Math.random();
  if ( X*X + Y*Y <= 1.0 ) inside+=4;
  }

如果我们探测单位圆的第一象限，我们不需要通过 size 改变动态范围，我们也可以测试 2 形式的距离，它摆脱了 sqrt .这些更改应该会提高精度和速度。

不是 JAVASCRIPT 编码员，所以我不知道您使用的是什么数据类型，但您需要确保您没有超过它的精度。在这种情况下，您需要添加更多计数器变量以减轻其负载。有关详细信息，请参阅：.

因为你的数字相当大，我敢打赌你已经越界了（应该没有小数部分，尾随零也很可疑）例如 32 位 float 只能存储整数

2^23 = 8388608

而你的 698,565,481,000,000 远高于此，因此即使对此类变量进行 ++ 操作也会导致精度损失，当指数太大时，它甚至会停止添加...

在整数上这不是问题，但是一旦你根据内部格式越过边界，值就会绕零或取反......但我怀疑是这种情况，因为结果将与 PI 相去甚远。

Answer 2

你可以很容易地估计出你应该得到什么样的错误（错误条），这就是 Monte Carlo 的美妙之处。为此，您必须计算第二动量并估计方差和 std.deviation。好的是，收集的值将与您收集的平均值相同，因为您只是在 1 后 1 后加 1。

然后您可以获得模拟西格玛的估计值，以及所需值的误差线。抱歉，我了解的不够Javascript，所以这里的代码是 C#:

using System;

namespace Pi
{
    class Program
    {
        static void Main(string[] args)
        {
            ulong N = 1_000_000_000UL; // number of samples
            var rng = new Random(312345); // RNG

            ulong v  = 0UL; // collecting mean values here
            ulong v2 = 0UL; // collecting squares, should be the same as mean
            for (ulong k = 0; k != N; ++k) {
                double x = rng.NextDouble();
                double y = rng.NextDouble();

                var r = (x * x + y * y < 1.0) ? 1UL : 0UL;

                v  += r;
                v2 += r * r;
            }

            var mean = (double)v / (double)N;
            var varc = ((double)v2 / (double)N - mean * mean ) * ((double)N/(N-1UL)); // variance
            var stdd = Math.Sqrt(varc); // std.dev, should be sqrt(Pi/4 (1-Pi/4))
            var errr = stdd / Math.Sqrt(N);

            Console.WriteLine($"Mean = {mean}, StdDev = {stdd}, Err = {errr}");

            mean *= 4.0;
            errr *= 4.0;

            Console.WriteLine($"PI (1 sigma) = {mean - 1.0 * errr}...{mean + 1.0 * errr}");
            Console.WriteLine($"PI (2 sigma) = {mean - 2.0 * errr}...{mean + 2.0 * errr}");
            Console.WriteLine($"PI (3 sigma) = {mean - 3.0 * errr}...{mean + 3.0 * errr}");
        }
    }
}

在 10⁹ 个样本之后我得到了

Mean = 0.785405665, StdDev = 0.410540627166729, Err = 1.29824345388086E-05
PI (1 sigma) = 3.14157073026184...3.14167458973816
PI (2 sigma) = 3.14151880052369...3.14172651947631
PI (3 sigma) = 3.14146687078553...3.14177844921447

这看起来是正确的。很容易看出，在理想情况下，方差等于 (Pi/4)*(1-Pi/4)。 v2确实没必要计算，模拟后设置为v即可。

坦率地说，我不知道为什么你得到的不是预期的。求和中的精度损失可能是答案，或者我怀疑，由于播种和重叠序列，您的模拟没有产生独立样本（因此实际 N 远低于 900 万亿）。

但是使用这种方法可以控制错误并检查计算的进展情况。

更新

我输入了您的数字，表明您明显低估了价值。代码

    N  = 893_547_800_000UL;
    v  = 701_766_448_388UL;
    v2 = v;

    var mean = (double)v / (double)N;
    var varc = ((double)v2 / (double)N - mean * mean ) * ((double)N/(N-1UL)); 
    var stdd = Math.Sqrt(varc); // should be sqrt(Pi/4 (1-Pi/4))
    var errr = stdd / Math.Sqrt(N);

    Console.WriteLine($"Mean = {mean}, StdDev = {stdd}, Err = {errr}");

    mean *= 4.0;
    errr *= 4.0;

    Console.WriteLine($"PI (1 sigma) = {mean - 1.0 * errr}...{mean + 1.0 * errr}");
    Console.WriteLine($"PI (2 sigma) = {mean - 2.0 * errr}...{mean + 2.0 * errr}");
    Console.WriteLine($"PI (3 sigma) = {mean - 3.0 * errr}...{mean + 3.0 * errr}");

并输出

Mean = 0.785370909522692, StdDev = 0.410564786603016, Err = 4.34332975349809E-07
PI (1 sigma) = 3.14148190075886...3.14148537542267
PI (2 sigma) = 3.14148016342696...3.14148711275457
PI (3 sigma) = 3.14147842609506...3.14148885008647

所以，很明显你在某处有问题（代码？表示中的精度丢失？求和中的精度丢失？repeated/non-independent 采样？）

使用 Monte Carlo 模拟限制计算 π

Calculating π using a Monte Carlo Simulation limitations

javascript

math

montecarlo

pi