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.

Show comment
Hide comment
Owner

linw1995 commented May 21, 2017

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