{{ message }}

Instantly share code, notes, and snippets.

# linw1995/Monte Carle Intergation.go

Last active May 27, 2017
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()