From fc59aecb1473e719cd673f53901d212545e305b4 Mon Sep 17 00:00:00 2001 From: lifeyst <33717659+lifeyst@users.noreply.github.com> Date: Wed, 6 Dec 2017 22:07:00 +0800 Subject: [PATCH] python(13~15,17~21) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 于晟焘 学号：201730310131 邮箱：lifeyst0824@163.com --- p13.py | 21 +++++++++++ p14.py | 15 ++++++++ p15.py | 12 +++++++ p17.py | 26 ++++++++++++++ p18.py | 67 ++++++++++++++++++++++++++++++++++ p19.py | 51 ++++++++++++++++++++++++++ p20.py | 50 ++++++++++++++++++++++++++ p21.py | 111 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 8 files changed, 353 insertions(+) create mode 100644 p13.py create mode 100644 p14.py create mode 100644 p15.py create mode 100644 p17.py create mode 100644 p18.py create mode 100644 p19.py create mode 100644 p20.py create mode 100644 p21.py 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()