We have 3 types of Distribution which we have talked about:-
- Uniform Distribution
- Normal Distribution
- Exponential Distribution
A Uniform Distribution will look like this:-
The probability will be defined as
We can plot two type of graph for uniform distribution, and there respective plots
- Continuous Graph
- Discrete Graph.
Both the graph looks a quite alike, the only difference being the value of
X axis and
A Normal Distribution will look like this, we some empirical rules:-
A Normal Distribution have these property:-
- It always peeks at the mean.
- Falls of symmetrically.
- Normal Distribution can be nicely characterized by 2 parameters, Mean and Standard Deviation.
The function for normal distribution is denoted by:-
The parameter in this definition is the mean or expectation of the distribution (and also its median and mode). The parameter is its standard deviation; its variance is therefore . A random variable with a Gaussian distribution is said to be normally distributed and is called a normal deviate.
In probability theory and statistics, the exponential distribution (a.k.a. negative exponential distribution) is the probability distribution that describes the time between events in a Poisson process, i.e. a process in which events occur continuously and independently at a constant average rate.
The function for Exponential distribution will be:-
Monte Carlo methods (or Monte Carlo experiments) are a broad class of computational algorithms that rely on repeated random sampling to obtain numerical results.
Monty Hall Problem
One of the good application of Monte Carlo simulation is to solve the Monty Hall Problem, as the code shows below:-
import random def chooseDoor(): return random.choice([1,2,3]) def playMontyHall(numTrails = 1000): stayWins = 0 switchWins = 0 for trails in range(numTrails): prizeDoor = chooseDoor() playerDoor = chooseDoor() if prizeDoor == playerDoor: stayWins += 1 elif prizeDoor != playerDoor: switchWins += 1 print "Stay Wins: ", stayWins/float(numTrails) print "Switch Wins: ", switchWins/float(numTrails) playMontyHall()
Estimation of PI.
We can estimate PI using the below code:-
import math import random import pylab ##### # # Computing Pi # # Square with 2r length Sides # Inscribe circle with radius r # Area of square = (2r)*2 = 4r^2 # Area of circle = pi*r^2 # Ratio of circle area to sqare area is (pi*r^2)/(4r^2) = pi/4 # # Implication: Of N random points picked from inside square, N*pi/4 # will be inside circle # # So if M = number of points inside circle # M = N*pi/4 # pi = M/N * 4 def randomPoints(r): x = random.uniform(-r,r) y = random.uniform(-r,r) return (x,y) def makePoints(r,n): """ Make n random points inside square with 2*r side. """ points =  for i in range(n): points.append(randomPoints(r)) return points def inCircle(r,points): x = points y = points return x ** 2 + y ** 2 <= r ** 2 def numInCircle(r,points): """ Figure out no of points inside circle of radius r """ count = 0 for point in points: if inCircle(r,point): count += 1 return count def computePi(numPoints,points = None): """ Computes Pi using Monte Carlo simulation of n points """ if points is None: points = makePoints(1.0,numPoints) inCircle = numInCircle(1.0, points) return float(inCircle)/float(numPoints) * 4.0 def runTrails(numTrailsPerPoints,numPointsList): results =  for numPoints in numPointsList: print numPoints for trails in range(numTrailsPerPoints): results.append((numPoints,computePi(numPoints))) return results def plotPi(trails,trailsResults): numPoints =  results =  for result in trailsResults: numPoints.append(result) results.append(result) pylab.figure() pylab.clf() pylab.scatter(numPoints, results, c="r") pylab.plot(trails,[math.pi for trails in trails], c="b") pylab.xlabel("Number of Points") pylab.ylabel("Pi") pylab.title("Pi Vs Number of Points") pylab.show() numTrailsPerPoints = 50 numPointsList = range(10,10000,1000) trailsResults = runTrails(numTrailsPerPoints, numPointsList) plotPi(numPointsList, trailsResults)
The plot will look like this:-
- The blue line on the plot is the exact mathematical calculation of Pi which a computer can do.
- The Red Dots are the estimation of PI, so if we increase the no of trails we see the deviation is much smaller.
We can also check where the needles are dropped.
Here is the code:-
def PlotPiScatter(r,n): points = makePoints(1.0, n) piEstimate = computePi(n,points) squarePoints_X =  squarePoints_Y =  circlePoints_X =  circlePoints_Y =  for point in points: if inCircle(r,point): circlePoints_X.append(point) circlePoints_Y.append(point) else: squarePoints_X.append(point) squarePoints_Y.append(point) pylab.figure() pylab.clf() pylab.scatter(squarePoints_X, squarePoints_Y,c="r") pylab.scatter(circlePoints_X, circlePoints_Y,c="b") pylab.title("With " +str(n) +" points ") pylab.axis([-1.5,1.5,-1.5,1.5]) pylab.text(-1.4, -1.4, "Pi is estimated to be " +str(piEstimate)) pylab.show()
- No of needles with just 10 points
- No of needles with just 100 points
- No of needles with just 1000 points
- No of needles with just 10000 points
- No of needles with just 1000000 points
We can use the same method of needle drop to find the area of a function using integration.
Kindly understand this code:-
import math import random import numpy import pylab # float range def frange(start,stop,step): l =  for i in range(int((stop - start)/step)): l.append(start + step*i) return l # # rSquared # def mse(measured,predicted): """Compute Sum of Residual Square""" sum_sq = 0 for i in xrange(len(measured)): sum_sq = (measured[i] - predicted[i]) ** 2 return float(sum_sq) def sstot(measured): """Compute total of sum square""" nMean = sum(measured) /float(len(measured)) ntot = 0 for m in measured: ntot += (m - nMean) ** 2 return ntot def rSquared(measured,predicted): """ measured: list of measured values predicted: list of predicted values """ SSerr = mse(measured, predicted) SStot = sstot(measured) return 1.0 - SSerr/SStot def makeCurve(f,xs): ys =  for x in xs: ys.append(f(x)) return ys def addNoise(ys, stddev = 1): ysp =  for y in ys: ysp.append(random.gauss(y,stddev)) return ysp def makeObservations(xs,f,noise): # Make a theoritical curve ys = makeCurve(f,xs) # Add some noise to the simulate environment nys = addNoise(ys) return nys def showExamplePolyFit(xs,ys,fitDegree1 = 1,fitDegree2 = 2): pylab.figure() pylab.plot(xs,ys,'r.',ms=2.0,label = "measured") # poly fit to noise coeeff = numpy.polyfit(xs, ys, fitDegree1) # Predict the curve pys = numpy.polyval(numpy.poly1d(coeeff), xs) se = mse(ys, pys) r2 = rSquared(ys, pys) pylab.plot(xs,pys, 'g--', lw=5,label="%d degree fit, SE = %0.10f, R2 = %0.10f" %(fitDegree1,se,r2)) # Poly fit to noise coeeffs = numpy.polyfit(xs, ys, fitDegree2) # Predict the curve pys = numpy.polyval(numpy.poly1d(coeeffs), xs) se = mse(ys, pys) r2 = rSquared(ys, pys) pylab.plot(xs,pys, 'b--', lw=5,label="%d degree fit, SE = %0.10f, R2 = %0.10f" %(fitDegree2,se,r2)) pylab.legend() def f(x): return x**3 + 5*x xs = frange(-2,2,0.001) ys = makeObservations(xs, f, 3) showExamplePolyFit(xs, ys, 1,2) showExamplePolyFit(xs, ys, 2,3) showExamplePolyFit(xs, ys, 3,4) showExamplePolyFit(xs, ys, 4,5) pylab.show()
These are graph which we will get.
- Regression w.r.t 1-degree and 2-Degree
- Regression w.r.t 2-degree and 3-Degree
- Regression w.r.t 3-degree and 4-Degree
- Regression w.r.t 4-degree and 5-Degree
The Standard error are not very good for the range
2, change it to
5 and see the difference.