Instantly share code, notes, and snippets.

Embed
What would you like to do?
Monte Carlo method applied to approximating the value of π.
# 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"
"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")
}
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))
}
}
# 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()
@linw1995

This comment has been minimized.

Owner

linw1995 commented May 21, 2017

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment