diff --git a/p13.py b/p13.py
new file mode 100644
index 0000000..0d4b5e4
--- /dev/null
+++ b/p13.py
@@ -0,0 +1,21 @@
+#-*- coding: utf-8 -*-
+import numpy as np
+
+x = np.array([7,38,4,23,18])
+
+def compMean(x):
+	mu = np.sum(np.sum(x))/x.shape[0]
+	return mu
+
+def compVariance(x):
+	x1 = np.square(x-mu)
+	sigma2 = compMean(x1)
+	return sigma2
+
+mu=compMean(x)
+sigma2 = compVariance(x)
+
+if __name__=="__main__":
+	print('The Expectation / Mean:\t %.1f \n'%mu)
+	print('The Variance is:\t %.1f \n'%sigma2)
+	print('The Standard Deviation is: %.1f \n'%np.sqrt(sigma2))
diff --git a/p14.py b/p14.py
new file mode 100644
index 0000000..932952c
--- /dev/null
+++ b/p14.py
@@ -0,0 +1,15 @@
+#-*- coding: utf-8 -*-
+import numpy as np
+from p13 import compMean
+
+a = 2
+b = 4
+x = np.array([7,38,4,23,18])
+y = a*x + b
+
+if __name__ == "__main__":
+	xmu = compMean(x)
+	ymu = compMean(y)
+
+print('a * Xmu + b \t= %.2f \n'%(a * xmu + b))
+print('Ymu \t\t= %.2f \n'%ymu)
diff --git a/p15.py b/p15.py
new file mode 100644
index 0000000..08c9ad9
--- /dev/null
+++ b/p15.py
@@ -0,0 +1,12 @@
+#-*- coding: utf-8 -*-
+import numpy as np
+import p13
+
+x = np.array([7,38,4,23,18])
+x2 = np.square(x)
+x2mu = p13.compMean(x2)
+xmu = p13.compMean(x)
+xvar = p13.compVariance(x)
+
+print('The Variance of X\t= %.2f \n'%xvar)
+print('E[X^2]-E[X]^2 \t\t= %.2f \n'%(x2mu - np.square(xmu)))
diff --git a/p17.py b/p17.py
new file mode 100644
index 0000000..8d7488a
--- /dev/null
+++ b/p17.py
@@ -0,0 +1,26 @@
+#-*- coding: utf-8 -*-
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+
+def gaussian1d(mu, sigma, x):
+	p = 1/(sigma*np.sqrt(2*np.pi))*np.exp((-np.square(x-mu))/(2*(sigma**2)))
+	return p
+	
+x = np.arange(-2,12.01,0.01)
+mu = 5
+sigma = 0.5
+p = gaussian1d(mu, sigma, x)
+if __name__ == "__main__":
+	print('The integration of p with respect to x is %.f\n' %np.trapz(p,x,axis=0))
+	plt.figure('One D Gaussian distribution')
+	plt.plot(x,p,'b-',linewidth=4) 
+	plt.plot([mu,mu],[np.min(p),np.max(p)],'r-',linewidth=4)
+ 
+	plt.title(r'$mu=%s,sigma=%s$'%(mu,sigma))
+	plt.xlabel('x')
+	plt.ylabel('P(x)') 
+	plt.grid('on') 
+	plt.show()
+
diff --git a/p18.py b/p18.py
new file mode 100644
index 0000000..bb3cab4
--- /dev/null
+++ b/p18.py
@@ -0,0 +1,67 @@
+#-*- coding: utf-8 -*-
+
+import matplotlib.pyplot as plt
+import numpy as np
+from mpl_toolkits.mplot3d import Axes3D
+
+def gaussian2d(mu, Sigma, x, y):
+	xLen = len(x)
+	yLen = len(y)
+	if xLen!=yLen:
+		error('The size of x and y should be the same!')
+	if len(mu)!=2:
+		error('The size of mu should be (2 * 1)!')
+	if Sigma.shape[0]!=Sigma.shape[1]:
+		error('The covariance matrix should be symmetric!')
+	if Sigma.shape[0]!=2:
+		error('The size of covariance matrix is invalid!')
+	if np.linalg.det(Sigma) < 2.2204*10**(-16):
+		error('Covariance must be positive semidefinite!')
+	pMat = np.zeros((xLen,yLen))
+	[xMat,yMat] = np.meshgrid(x,y)
+	[xMat,yMat] = [np.transpose(xMat),np.transpose(yMat)]
+	for i in range(xLen):
+		for j in range(yLen):
+			xMinusMu = np.array([[xMat[i][j] - mu[0]],[ yMat[i][j] - mu[1]]])
+			p = np.exp(-0.5*np.dot(np.dot(np.transpose(xMinusMu),np.linalg.inv(Sigma)),xMinusMu))
+			p = p*np.linalg.det(2*np.pi*Sigma)**(-0.5)
+			pMat[i][j] = p
+	return pMat, xMat, yMat
+			
+mu = np.array([6,5])
+sigmax = 2
+sigmay = 1
+rho = -0.2
+sigmaxy = sigmax * sigmay * rho
+Sigma = np.array([[sigmax**2,sigmaxy],[sigmaxy,sigmay**2]])
+x = np.arange(0,12.01,0.01)
+y = x
+[pMat, xMat, yMat] = gaussian2d(mu, Sigma, x, y)
+if __name__ == "__main__":
+	print('The integration of p with respect to x and y is %.4f\n'%np.trapz(np.trapz(pMat,x),y))
+
+	fig = plt.figure('2D Gaussian distribution1')
+	ax = Axes3D(fig)
+	ax.view_init(30,35)
+	ax.plot_surface(xMat, yMat, pMat, rstride=1, cstride=1, cmap='rainbow')
+	ax.set_xlabel("x")
+	ax.set_ylabel("y")
+	ax.set_zlabel("Probability")
+	plt.axis([max(x),min(x),max(y),min(y)])
+
+	fig = plt.figure('2D Gaussian distribution2')
+	ax = Axes3D(fig)
+	ax.view_init(100,-270)
+	ax.plot_surface(xMat, yMat, pMat, rstride=1, cstride=1, cmap='rainbow')
+	plt.xlabel('x')
+	plt.ylabel('y')
+	plt.axis([max(x),min(x),max(y),min(y)])
+
+	plt.figure('The contour of 2D Gaussian distribution')
+	plt.plot (mu[0], mu[1], 'xg', markersize= 6, linewidth= 1)
+	C = plt.contour(xMat, yMat, pMat, 8)
+	plt.clabel(C, inline = True, fontsize = 10)
+	plt.xlabel('x')
+	plt.ylabel('y')
+	plt.axis([min(x),max(x),min(y),max(y)])
+	plt.show()
diff --git a/p19.py b/p19.py
new file mode 100644
index 0000000..94a16fc
--- /dev/null
+++ b/p19.py
@@ -0,0 +1,51 @@
+#-*- coding: utf-8 -*-
+import matplotlib.pyplot as plt
+import numpy as np
+from p17 import gaussian1d
+
+
+def product1d(fx,fy,x):
+	fxLen = len(fx)
+	fyLen = len(fy)
+	xLen = len(x)
+	if fxLen!=fyLen:
+		error('Probability vectors must be of the same length')
+	if xLen!= fxLen:
+		error('X range and probability vectors must be of the same length')
+	fxy = np.multiply(fx, fy)
+	fxy = fxy/np.trapz(fxy,x)
+	return fxy
+
+def convolute1d(fx,fy,x):
+	fxy = np.convolve(fx,fy,mode='same')
+	fxy = fxy / np.trapz(fxy,x)
+	return fxy
+	
+x = np.arange(-10,10.01,0.01)
+if __name__ == "__main__":
+	prior         = gaussian1d(2, 1, x) #a Gaussian Prior of mean 2 and SD 1
+	motionModel   = gaussian1d(1, 2, x) #a Gaussian motion model of mean 1 and SD 2
+	sensorModel   = gaussian1d(5, 1, x) # a Gaussian observation model of mean 5 and SD 1
+	probAfterMove  = convolute1d(motionModel, prior, x) #prediction is convolution
+	probAfterSense = product1d(sensorModel, probAfterMove, x) #update is multiplication
+	
+	plt.figure('Bayes filtering')
+	plt.subplot(211)
+	plt.plot(x, prior, 'b-', linewidth = 2,label="Prior")
+	plt.plot(x, motionModel, 'g-', linewidth = 2,label="Motion model")
+	plt.plot(x, probAfterMove, 'r-', linewidth= 2,label="Predicted prob")
+	plt.legend(loc=1,ncol=1)
+	plt.title('Bayes filtering (Prediction)')
+	plt.ylabel('Probability')
+	plt.axis([-10,10,0,0.5])
+	plt.grid('on')
+	plt.subplot(212)
+	plt.title('Update')
+	plt.plot(x, probAfterMove, 'r-', linewidth= 2,label="Predicted prob")
+	plt.plot(x, sensorModel, 'c-', linewidth= 2,label="Sensor model")
+	plt.plot(x, probAfterSense, 'k-', linewidth=2,label="Posterior")
+	plt.legend(loc=1,ncol=1)
+	plt.ylabel('Probability')
+	plt.axis([-10,10,0,0.5])
+	plt.grid('on') 
+	plt.show()
diff --git a/p20.py b/p20.py
new file mode 100644
index 0000000..fcf83d6
--- /dev/null
+++ b/p20.py
@@ -0,0 +1,50 @@
+#-*- coding: utf-8 -*-
+import matplotlib.pyplot as plt
+import numpy as np
+from p17 import gaussian1d
+import p19
+
+def convolute1dParam(muPrior, sdPrior, muMotion, sdMotion):
+	muAfterMove = muPrior + muMotion;
+	sdAterMove  = np.sqrt(sdPrior**2 + sdMotion**2)
+	return muAfterMove,sdAterMove
+
+def product1dParam(muAfterMove, sdAfterMove, muSensor, sdSensor):
+	muPost = (muAfterMove * sdSensor**2 + muSensor * sdAfterMove**2)/(sdAfterMove**2 + sdSensor**2)
+	sdPost = (sdSensor**2 * sdAfterMove**2)/(sdAfterMove**2 + sdSensor**2)
+	return muPost,sdPost
+
+x = np.arange(-100,100.001,0.001)
+prior         = gaussian1d(2, 1, x)    # a Gaussian Prior of mean 2 and SD 1
+motionModel   = gaussian1d(1, 2, x)    # a Gaussian motion model of mean 1 and SD 2
+sensorModel   = gaussian1d(5, 1, x)    # a Gaussian observation model of mean 5 and SD 1   
+
+probAfterMove  = p19.convolute1d(motionModel, prior, x)      #prediction is convolution
+probAfterSense = p19.product1d(sensorModel, probAfterMove, x) #update is multiplication
+[muAfterMove, sdAfterMove]   = convolute1dParam(2, 1, 1, 2)
+[muAfterSense, sdAfterSense] = product1dParam(muAfterMove, sdAfterMove, 5, 1)
+
+probAfterMoveParam  = gaussian1d(muAfterMove, sdAfterMove, x)
+probAfterSenseParam = gaussian1d(muAfterSense, sdAfterSense, x)
+probAfterSenseParamSeq = p19.product1d(sensorModel, probAfterMoveParam, x)
+
+plt.figure('Bayes filtering')
+plt.subplot(211)
+plt.grid('on')
+plt.plot(x, probAfterMove, 'r-', linewidth= 6,label="Predicted prob")
+plt.plot(x, probAfterMoveParam, 'g:', linewidth=3,label="Predicted prob \n with params")
+plt.legend(loc=1,ncol=1)
+plt.title('Bayes filtering (Prediction)')
+plt.ylabel('Probability')
+plt.axis([-5,10,0,0.2])
+
+plt.subplot(212)
+plt.grid('on')
+plt.title('Update')
+plt.plot(x, probAfterSense, 'k-', linewidth= 6,label="Posterior")
+plt.plot(x, probAfterSenseParam, 'c-', linewidth= 2,label="Posterior with params as input and output")
+plt.plot(x, probAfterSenseParamSeq, 'y-',linewidth=2,label="Posterior with params as input")
+plt.legend(loc=0,ncol=1)
+plt.ylabel('Probability')
+plt.axis([-5,10,0,0.5])
+plt.show()
diff --git a/p21.py b/p21.py
new file mode 100644
index 0000000..92f2e11
--- /dev/null
+++ b/p21.py
@@ -0,0 +1,111 @@
+#-*- coding: utf-8 -*-
+import matplotlib.pyplot as plt
+import numpy as np
+from p18 import gaussian2d
+from mpl_toolkits.mplot3d import Axes3D
+from scipy import signal
+from scipy import misc
+
+def product2d(p1, p2, x, y):
+	p = np.multiply(p1, p2)
+	p = p / np.trapz(np.trapz(p,x),y)
+	return p
+	
+def makeSigma(sigmax, sigmay, rho):
+	sigmaxy = sigmax * sigmay * rho
+	Sigma = np.array([[sigmax**2,sigmaxy],[sigmaxy,sigmay**2]])
+	return Sigma
+	
+def convolute2d(p1, p2, x, y):
+	p= signal.convolve2d(p1,p2,boundary='symm',mode='same')
+	p = p / np.trapz(np.trapz(p,x),y)
+	return p
+
+def maxMat(g):
+	r=g.argmax(axis=0)
+	temp=g[r,range(g.shape[1])]
+	index_Row=np.argsort(temp)
+	c=index_Row[len(index_Row)-1]
+	r=r[c]
+	return r,c
+	
+x = np.arange(-2,10.02,0.02)
+y = x
+mu1 = np.array([5,6])
+mu2 = np.array([3,2])
+mu3 = np.array([6,4])
+#Prior
+[priorMat, xMat, yMat] = gaussian2d(mu1, makeSigma(0.6, 0.9, 0.7), x, y)
+#Motion model
+[motMat,xMat1,yMat1] = gaussian2d(mu2, makeSigma(0.4, 0.7, 0.3), x, y)
+#Observation model
+[obsMat,xMat2,yMat2] = gaussian2d(mu3, makeSigma(0.6, 0.5, 0.5), x, y)
+#Prediction
+predMat = convolute2d(priorMat, motMat, x, y)
+#Update
+postMat = product2d(predMat, obsMat, x, y)
+
+fig = plt.figure('Prediction')
+ax = Axes3D(fig)
+ax.view_init(30,35)
+ax.plot_surface(xMat, yMat, priorMat, rstride=1, cstride=1, cmap='rainbow')
+ax.text(5,6, np.max(np.max(priorMat)),'Prior',family = 'monospace', fontsize = 10)
+ax.plot_surface(xMat, yMat, motMat, rstride=1, cstride=1, cmap='rainbow')
+ax.text(3,2, np.max(np.max(motMat)),'Motion model',family = 'monospace', fontsize = 10)
+ax.plot_surface(xMat, yMat, predMat, rstride=1, cstride=1, cmap='rainbow')
+[r, c] = maxMat(predMat)
+ax.text(xMat[r,c],yMat[r,c], predMat[r, c],'Predicted',family = 'monospace', fontsize = 10)
+ax.set_xlabel("x")
+ax.set_ylabel("y")
+ax.set_zlabel("Probability")
+plt.axis([max(x),min(x),max(y),min(y)])
+
+plt.figure('Prediction contour')
+plt.grid('on')
+C = plt.contour(xMat, yMat, priorMat, 8)
+plt.clabel(C, inline = True, fontsize = 10)
+plt.text(5,6,'Prior')
+C = plt.contour(xMat, yMat, motMat, 8)
+plt.clabel(C, inline = True, fontsize = 10)
+plt.text(3,2,'Motion model')
+C = plt.contour(xMat, yMat,predMat, 8)
+plt.clabel(C, inline = True, fontsize = 10)
+[r, c] = maxMat(predMat)
+plt.text(xMat[r,c],yMat[r,c], 'Predicted')
+plt.xlabel('x')
+plt.ylabel('y')
+plt.axis([min(x),max(x),min(y),max(y)])
+
+fig = plt.figure('Update')
+ax = Axes3D(fig)
+ax.view_init(30,35)
+ax.plot_surface(xMat, yMat,predMat, rstride=1, cstride=1, cmap='rainbow')
+[r, c] = maxMat(predMat)
+ax.text(xMat[r,c],yMat[r,c], predMat[r, c],'Predicted',family = 'monospace', fontsize = 10)
+ax.plot_surface(xMat, yMat,obsMat, rstride=1, cstride=1, cmap='rainbow')
+ax.text(6,4, np.max(np.max(obsMat)),'Observation model',family = 'monospace', fontsize = 10)
+ax.plot_surface(xMat, yMat,postMat, rstride=1, cstride=1, cmap='rainbow')
+[r, c] = maxMat(postMat)
+ax.text(xMat[r,c],yMat[r,c], postMat[r, c],'Updated',family = 'monospace', fontsize = 10)
+ax.set_xlabel("x")
+ax.set_ylabel("y")
+ax.set_zlabel("Probability")
+plt.axis([max(x),min(x),max(y),min(y)])
+
+plt.figure('Update contour')
+plt.grid('on')
+C = plt.contour(xMat, yMat,predMat, 8)
+plt.clabel(C, inline = True, fontsize = 10)
+[r, c] = maxMat(predMat)
+plt.text(xMat[r,c],yMat[r,c], 'Predicted')
+C = plt.contour(xMat, yMat,obsMat, 8)
+plt.clabel(C, inline = True, fontsize = 10)
+ax.text(6,4, np.max(np.max(obsMat)),'Observation model',family = 'monospace', fontsize = 10)
+C = plt.contour(xMat, yMat,postMat, 8)
+plt.clabel(C, inline = True, fontsize = 10)
+[r, c] = maxMat(postMat)
+plt.text(xMat[r,c],yMat[r,c], 'Updated')
+plt.xlabel('x')
+plt.ylabel('y')
+plt.axis([min(x),max(x),min(y),max(y)])
+plt.show()