78 lines
2.3 KiB
Go
78 lines
2.3 KiB
Go
// Copyright (c) 2026 Petr Balvín <opensource@petrbalvin.org> (https://petrbalvin.org)
|
|
// SPDX-License-Identifier: MIT
|
|
|
|
// Command qmc compares quasi-random integration against plain Monte
|
|
// Carlo on the same two-dimensional integral. Sobol points are a
|
|
// digital lattice: every block of 2^m points stratifies each
|
|
// coordinate exactly, so the error decays far faster than the
|
|
// 1/sqrt(n) of random sampling, and Halton sits in between.
|
|
//
|
|
// Usage: go run ./examples/qmc
|
|
package main
|
|
|
|
import (
|
|
"fmt"
|
|
"log"
|
|
"math"
|
|
|
|
"sourcedock.dev/petrbalvin/tensor"
|
|
)
|
|
|
|
// f is the integrand: smooth, with its curvature spread over the unit
|
|
// square. The exact value is (1-e^-1)*sqrt(pi)/2*erf(1), the product
|
|
// of the x integral and the error function integral over y.
|
|
func f(x, y float64) float64 { return math.Exp(-x - y*y) }
|
|
|
|
const exact = 0.4720828881800443 // (1-e^-1)*sqrt(pi)/2*erf(1)
|
|
|
|
// estimate integrates f over [0,1]^2 from an (n,2) point set.
|
|
func estimate(pts *tensor.Array, n int) float64 {
|
|
s := 0.0
|
|
for i := range n {
|
|
x, err := tensor.FloatAt(pts, i, 0)
|
|
if err != nil {
|
|
log.Fatal(err)
|
|
}
|
|
y, err := tensor.FloatAt(pts, i, 1)
|
|
if err != nil {
|
|
log.Fatal(err)
|
|
}
|
|
s += f(x, y)
|
|
}
|
|
return s / float64(n)
|
|
}
|
|
|
|
func main() {
|
|
fmt.Printf("integral of exp(-x - y^2) over the unit square, exact %.10f\n\n", exact)
|
|
fmt.Println(" points Monte Carlo Halton Sobol")
|
|
for _, n := range []int{64, 256, 1024, 4096, 16384} {
|
|
// Monte Carlo: uniform draws from the seeded generator.
|
|
g := tensor.NewGenerator(int64(n))
|
|
mc, err := tensor.Floats(g, 2*n)
|
|
if err != nil {
|
|
log.Fatal(err)
|
|
}
|
|
mcPts, err := tensor.Reshape(mc, n, 2)
|
|
if err != nil {
|
|
log.Fatal(err)
|
|
}
|
|
// Halton and Sobol from the first point on; Sobol skips its
|
|
// origin point exactly as Halton does.
|
|
hal, err := tensor.HaltonPoints(n, 2, 0)
|
|
if err != nil {
|
|
log.Fatal(err)
|
|
}
|
|
sob, err := tensor.SobolPoints(n, 2, 0)
|
|
if err != nil {
|
|
log.Fatal(err)
|
|
}
|
|
eMC := math.Abs(estimate(mcPts, n) - exact)
|
|
eHal := math.Abs(estimate(hal, n) - exact)
|
|
eSob := math.Abs(estimate(sob, n) - exact)
|
|
fmt.Printf(" %6d %.3e %.3e %.3e\n", n, eMC, eHal, eSob)
|
|
}
|
|
fmt.Println()
|
|
fmt.Println("the quasi-random errors collapse with n; the Monte Carlo")
|
|
fmt.Println("error only shrinks as 1/sqrt(n) and stays noisy on top")
|
|
}
|