golang 蒙特卡罗方法应用于近似π的值。

Posted 2021-05-24

# coding:utf-8
from __future__ import division

import random

import numpy as np
import matplotlib.pyplot as plt


def dis(x, y):
    return np.sqrt(x**2 + y**2)


def main():
    M = int(4e4)
    pi_n = np.zeros(M)
    count_inner = 0
    for n in range(M):
        x = random.random()
        y = random.random()
        if dis(x, y) < 1:
            count_inner += 1
        pi_n[n] = count_inner / (n + 1) * 4
    error = pi_n - np.pi
    error = np.abs(error)
    plt.plot(error)
    plt.ylabel('Error')
    plt.xlabel('n')
    plt.show()

if __name__ == '__main__':
    main()

package main

import (
	"fmt"
	"math"
	"math/rand"

	"github.com/gonum/plot"
	"github.com/gonum/plot/plotter"
	"github.com/gonum/plot/plotutil"
	"github.com/gonum/plot/vg"
)

type points [][]float64
type point []float64

func appendIfMissing(slice []float64, i float64) []float64 {
	for _, ele := range slice {
		if ele == i {
			return slice
		}
	}
	return append(slice, i)
}

func unique(array []float64) (rv []float64) {
	for _, v := range array {
		rv = appendIfMissing(rv, v)
	}
	return
}

func randUniform(low, high float64) float64 {
	return rand.Float64()*(high-low) + low
}

func randsign() (sign float64) {
	if rand.Float64() > 0.5 {
		sign = 1
	} else {
		sign = -1
	}
	return
}

func ranpt(low, high point) (pt point) {
	length := len(low)
	pt = make([]float64, length)
	for i := 0; i < length; i++ {
		pt[i] = randUniform(low[i], high[i])
	}
	return
}

func meanAndVariance(array []float64) (mean, variance float64) {
	length := len(array)
	if length == 0 {
		return
	}
	var a2sum, asum float64
	for i := 0; i < length; i++ {
		asum += array[i]
		a2sum += array[i] * array[i]
	}
	mean = asum / float64(length)
	variance = a2sum/float64(length) - mean*mean
	return
}

func miser(f func(vals ...float64) float64, regn [2]point, N int, dith float64) (pts points, average, variance float64, rvN int) {
	const (
		MNBS = 60
		MNPT = 15
		PFAC = 0.1
		BIG  = math.MaxFloat64
		TINY = 1 / math.MaxFloat64
	)
	low := regn[0]
	high := regn[1]
	length := len(low)
	Npre := int(math.Max(float64(MNPT), float64(N)*PFAC))
	V := high[0] - low[0]
	for i := 1; i < length; i++ {
		V *= high[i] - low[i]
	}
	if N < MNBS {
		pts = make(points, N)
		fvals := make([]float64, N)
		for i := 0; i < N; i++ {
			pt := ranpt(low, high)
			pts[i] = pt
			fvals[i] = f(pt...)
		}
		average, variance = meanAndVariance(fvals)
		variance = math.Max(TINY, variance)
		rvN = N
		return pts, average, variance, rvN
	}
	//else
	mid := make(point, length)
	for i := 0; i < length; i++ {
		s := randsign() * dith
		mid[i] = (0.5+s)*low[i] + (0.5-s)*high[i]
	}
	fvalsL := make([][]float64, length)
	fvalsR := make([][]float64, length)
	for i := 0; i < Npre; i++ {
		pt := ranpt(low, high)
		fval := f(pt...)
		for j := 0; j < length; j++ {
			if pt[j] <= mid[j] {
				fvalsL[j] = append(fvalsL[j], fval)
			} else {
				fvalsR[j] = append(fvalsR[j], fval)
			}
		}
	}
	var (
		sigmaL, sigmaR, sigmaLB, sigmaRB float64
		NL, NR                           int
	)
	sumB := BIG
	iB := -1
	sSums := make([]float64, length)
	variancesL := make([]float64, length)
	variancesR := make([]float64, length)
	for i := 0; i < length; i++ {
		_, variancesL[i] = meanAndVariance(fvalsL[i])
		_, variancesR[i] = meanAndVariance(fvalsR[i])
		sigmaL = math.Max(TINY, math.Pow(variancesL[i], 2.0/3))
		sigmaR = math.Max(TINY, math.Pow(variancesR[i], 2.0/3))
		sSums[i] = sigmaL + sigmaR
		if sSums[i] <= sumB {
			sumB = sSums[i]
			iB = i
			sigmaLB = sigmaL
			sigmaRB = sigmaR
		}
	}
	if iB == -1 || len(unique(sSums)) == 1 {
		iB = int(randUniform(0, float64(length)))
	}
	regnL := low[iB]
	regnMid := mid[iB]
	regnR := high[iB]
	fracL := math.Abs((regnMid - regnL) / (regnR - regnL))
	fracR := 1.0 - fracL
	if sigmaLB > 0 || sigmaRB > 0 {
		L := fracL * sigmaLB
		R := fracR * sigmaRB
		NL = MNPT + int(float64(N-Npre-2*MNPT)*L/(L+R))
	} else {
		NL = MNPT + int((N-Npre-2*MNPT)/2)
	}
	NR = N - Npre - NL
	lowTmp := make(point, length)
	highTmp := make(point, length)
	copy(lowTmp, low)
	copy(highTmp, high)
	highTmp[iB] = mid[iB]
	ptsL, aveL, varL, NL := miser(f, [2]point{lowTmp, highTmp}, NL, dith)
	lowTmp[iB] = mid[iB]
	highTmp[iB] = high[iB]
	ptsR, aveR, varR, NR := miser(f, [2]point{lowTmp, highTmp}, NR, dith)
	pts = append(ptsL, ptsR...)
	average = fracL*aveL + fracR*aveR
	variance = fracL*fracL*varL + fracR*fracR*varR
	rvN = NL + NR
	return pts, average, variance, rvN
}

func f(vals ...float64) float64 {
	x, y := vals[0], vals[1]
	if math.Sqrt(x*x+y*y) < 1 {
		return 1
	}
	return 0
}

func main() {
// 	for i := 3; i < 10; i++ {
// 		N := int(math.Pow10(i))
// 		_, average, variance, rvN := miser(f, [2]point{point{0, 0}, point{1, 1}}, N, 0.1)
// 		fmt.Printf("%d %.6f %.4e %.4e %d %d\n",i, average*4, 4*math.Sqrt(variance)/math.Sqrt(float64(N)), math.Sqrt(variance), N, rvN)
// 	}
	N := int(1e7)
	pts, average, variance, rvN := miser(f, [2]point{point{0, 0}, point{1, 1}}, N, 0.1)
	fmt.Printf("%.6f %.4e %.4e %d %d\n", average*4, 4*math.Sqrt(variance)/math.Sqrt(float64(N)), math.Sqrt(variance), N, rvN)
	plt, _ := plot.New()
	inner := make(plotter.XYs, len(pts))
	outter := make(plotter.XYs, len(pts))
	var innerN, outterN int
	for _, v := range pts {
		if f(v...) > 0 {
			inner[innerN].X = v[0]
			inner[innerN].Y = v[1]
			innerN++
		} else {
			outter[outterN].X = v[0]
			outter[outterN].Y = v[1]
			outterN++
		}
	}
	s1, _ := plotter.NewScatter(inner[:innerN])
	s2, _ := plotter.NewScatter(outter[:outterN])
	s1.Shape = plotutil.Shape(4)
	s2.Shape = plotutil.Shape(4)
	s1.Color = plotutil.Color(0)
	s2.Color = plotutil.Color(1)
	plt.Title.Text = "MISER Monte Carlo"
	plt.X.Label.Text = "X"
	plt.Y.Label.Text = "Y"
	plt.Add(s1, s2)
	plt.Save(4*vg.Inch, 4*vg.Inch, "distributions.png")
}

# coding:utf-8
from __future__ import division

import random

import matplotlib.animation as animation
import matplotlib.pyplot as plt
import numpy as np


def dis(x, y):
    return np.sqrt(x**2 + y**2)


def main():
    M = int(4e4)
    inner_x = np.zeros(M)
    inner_y = np.zeros(M)
    outer_x = np.zeros(M)
    outer_y = np.zeros(M)
    count_inner = count_outer = 0
    fig = plt.figure()
    ims =[]
    count = 0
    for n in range(M):
        x = random.random()
        y = random.random()
        if dis(x, y) < 1:
            inner_x[count_inner] = x
            inner_y[count_inner] = y
            count_inner += 1
        else:
            outer_x[count_outer] = x
            outer_y[count_outer] = y
            count_outer += 1
        if count >= M / 10 ** 2:
            im = (
                plt.scatter(
                  inner_x[:count_inner], inner_y[:count_inner],
                  c='r', marker='.', s=1),
                plt.scatter(
                  outer_x[:count_outer], outer_y[:count_outer],
                  c='b', marker='.', s=1),
                plt.text(-0.2, 0, '$n=$%5d' % n),
                plt.text(
                  -0.2, -0.035, '$\pi \\approx$%.7f' % (count_inner / n * 4)))
            ims.append(im)
            count = 0
        count += 1
    im_ani = animation.ArtistAnimation(fig, ims, interval=50, blit=True)
    plt.axis('equal')
    plt.show()


if __name__ == '__main__':
    main()

package main

import (
	"fmt"
	"math"
	"math/rand"
)

func f(x float64) float64 {
	return math.Sqrt(1 - x*x)
}
func main() {
	var pn, x, fn, fn2sum, fnsum, fsum, fnave, fn2ave float64
	for j := 1; j < 10; j++ {
		N := math.Pow10(j)
		for i := pn; i < N; i++ {
			x = rand.Float64()
			fsum += f(x)
			fn = fsum / (i + 1)
			fnsum += fn
			fn2sum += fn * fn
		}
		fnave = fnsum / N
		fn2ave = fn2sum / N
		sigmaN := math.Sqrt(fn2ave - fnave*fnave)
		sigma := 4 * sigmaN / math.Sqrt(N)
		pn = N
		fmt.Printf("%d,%.4e,%.4e,%.6f,%.6f\n", j, sigma, sigmaN, 4*fn, math.Abs(4*fn-math.Pi))
	}
}

5.15求π近似值

”正多边形逼近“法求π：核心思想是极限的思想。假设一个直径d为1的圆，只要求出该圆的周长C，就可以通过π=C/d方法求出π的值。所以关键是求出该圆的周长C。

 

”正多边形逼近“也叫做”割圆术“，当一个圆的内接正多边形边数越多时，其边长就越接近外接的圆周长。

设一个直径为1的圆的内接多边形边长为b，边数为i，周长为C=bi；当多边形的边数加倍后，新多边形边长为x=sqrt(2-2sqrt(1-b*b))/2,新多边形周长为：C=2ix。

#include <iostream>
#include<cmath>
using namespace std;

double getPI(int n) {//n指的是要迭代的次数（迭代次数越多越精确）
	int div,i=4;//div记录迭代次数
	double b=sqrt(2)/2.0;　　//直径为1的圆，内接正四边形的初始化，其边长为sqrt(2)/2，边数为4，即i
	double c=0.0;
	for(div=0;div<n;div++) {
		b=sqrt(2.0-2.0*sqrt(1.0-b*b))*0.5;//边数加倍后的边长
		i=i*2;　　　　　　　　　　　　　　　　//边数加倍
	}
	c=b*i;　　　　//计算多边形周长
	return c;
}

int main() {
	int n;
	double PI;
	cin>>n;
	PI=getPI(n);
	cout<<PI;
	return 0;
}