數理統計分析的MATLAB實現MATLAB基礎知識統計分析的基本概念、工具及推理基礎統計估計第一章MATLAB基礎知識1.1MATLAB概況1.1.1MATLAB發展歷程

1980年，美國的CleverMoler博士建立了MATLAB(MatrixLaboratory)，即矩陣實驗室。隨后1984年建立了名為MathWorks的軟件開發公司。早期的MATLAB只能做矩陣運算，繪圖也只能用星號描點的形式畫圖。而后經過不斷的改進發展，現在MATLAB已經成為國際最為流行的科學與工程計算軟件之一，它以其模塊化的計算方法、可視化與智能化的人機交互功能、豐富的矩陣運算、圖形繪制和數據處理函數，以及模塊化圖形組態的動態系統仿真工具Simulink，成為控制系統設計和仿真領域最受歡迎的軟件系統。1.1.2MATLAB的語言特點

MATLAB命令和數學中的符號、公式非常接近，可讀性強，容易掌握。MATLAB語言除了具有強大數值計算和圖形功能以外，還提供了應用于許多領域的工具箱。它與其他語言的接口能夠保證其與各種強大的計算機軟件相結合，可擴展性很強。MATLAB目前為止可以在各種類型的計算機上運行，程序也可以直接移植到其他機型上使用。可以說MATLAB是和機器類型及操作系統基本無關的軟件。MATLAB語言具有較高的運算精度，矩陣類運算可以達到10-15數量級的精度，符合一般科學與工程運算的要求。1.2MATLAB的功能

MATLAB產品族被廣泛地應用于信號圖像處理、控制系統設計、通信、系統仿真、虛擬現實等諸多領域。它的一大特性是有眾多面向具體應用的工具箱和仿真塊，包含了完整的函數集用來對信號和圖像處理、控制系統設計、神經網絡等特殊應用進行分析和設計。圖1-1世界正則矩陣地圖圖1-2朝鮮半島的數字高程模型圖1-3經過光照處理的三維圖MATLAB主要功能如下：（1）MATLAB。MATLAB是MathWorks公司所有產品的數值分析和圖形基礎環境，它將二維和三維圖形、MATLAB語言能力集成到一個單一的、易學易用的環境之中。（2）MATLABToolbox。工具箱是一系列專用的MATLAB函數庫，以解決特定領域的問題，它是開放的、可擴展的，也就是說用戶可以查看其中的算法或開發自己的算法。（3）MATLABCompiler。編譯器可以將MATLAB語言編寫的M-文件自動轉換成C或C++文件，支持用戶進行獨立應用開發。（4）Simulink。Simulink是結合了框圖界面和交互仿真能力的非線性動態系統仿真工具，它以MATLAB的核心數學、圖形和語言為基礎。（5）Stateflow。Stateflow與Simulink框圖模型相結合，描述復雜事件驅動系統的邏輯行為，驅動系統在不同的模式之間進行切換。（6）Real-TimeWorkshop。直接從Simulink框圖自動生成C代碼，整個代碼的生成可以根據需要進行完全定制。（7）SimulinkBlockset。專門為特定領域設計的Simulink功能塊的集合，用戶也可以利用已有的塊或自動編寫的C和MATLAB程序建立自己的塊。1.3MATLAB的開發環境1.3.1MATLAB桌面平臺

桌面平臺是各桌面組件的展示平臺，默認設置情況下的桌面平臺包括4個窗口，即命令窗口(CommandWindow)、命令歷史窗口(CommandHistory)、當前目錄窗口(CurrentDirectory)和工作空間窗口(Workspace)。此外，MATLAB還有編譯窗口、圖形窗口和幫助窗口等其他種類的窗口。1.3.2.運行方式

MATLAB提供了兩種運行方式，即命令行方式和M文件方式。1.3.3.MATLAB幫助系統1.3.4.工具箱第二章統計分析的基本概念、工具及推理基礎2.1變量與數據的基本概念2.1.1變量及其概率分布1.繪制正態分布的密度函數、分布函數程序如下：clearmu=2.5;sigma=0.6;x=(mu-4*sigma):0.005:(mu+4*sigma);y=normpdf(x,mu,sigma);f=normcdf(x,mu,sigma);plot(x,y,'-g',x,f,':b')legend('pdf','cdf',-1)表2-1線型和顏色控制符線形點標記顏色-實線:虛線-.點劃線--間斷線.點o小圓圈x叉子符+加號*星號s方格d菱形^朝上三角v朝下三角>朝右三角<朝左三角p五角星h六角星y黃m棕色c青色r紅色g綠色b藍色w白色k黑色2.正態分布參數μ和σ對密度曲線的影響程序如下：clearmu1=2.5;mu2=3;sigma1=0.5;sigma2=0.6;x=(mu2-4*sigma2):0.01:(mu2+4*sigma2);y1=normpdf(x,mu1,sigma1);y2=normpdf(x,mu2,sigma1);y3=normpdf(x,mu1,sigma1);y4=normpdf(x,mu1,sigma2);subplot(1,2,1)plot(x,y1,'-g',x,y2,'-b')xlabel('\fontsize{12}mu1<mu2,sigma1=sigma2')legend('mu1','mu2')subplot(1,2,2)plot(x,y3,'-g',x,y4,'-b')xlabel('\fontsize{12}mu1=mu2,sigma1<sigma2')legend('sigma1','sigma2')3.正態分布參數μ和σ對變量X取值規律的約束—3σ原則clear,clfX=linspace(-5,5,100);Y=normpdf(X,0,1);yy=normpdf([-3,-2,-1,0,1,2,3],0,1);plot(X,Y,'k-',[0,0],[0,yy(4)],'c-.')holdonplot([-2,-2],[0,yy(2)],'m:',[2,2],[0,yy(6)],'m:',[-2,-0.5],[yy(6),yy(6)],'m:',[0.5,2],[yy(6),yy(6)],'m:')plot([-1,-1],[0,yy(3)],'g:',[1,1],[0,yy(5)],'g:',[-1,-0.5],[yy(5),yy(5)],'g:',[0.5,1],[yy(5),yy(5)],'g:')plot([-3,-3],[0,yy(1)],'b:',[3,3],[0,yy(7)],'b:',[-3,-0.5],[yy(7),yy(7)],'b:',[0.5,3],[yy(7),yy(7)],'b:')holdofftext(-0.5,yy(6)+0.005,'\fontsize{14}95.44%')text(-0.5,yy(5)+0.005,'\fontsize{14}68.26%')text(-0.5,yy(7)+0.005,'\fontsize{14}99.74%')text(-3.2,-0.03,'\fontsize{10}μ-3σ')text(-2.2,-0.03,'\fontsize{10}μ-2σ')text(-1.2,-0.03,'\fontsize{10}μ-σ')text(-0.05,-0.03,'\fontsize{10}μ')text(0.8,-0.03,'\fontsize{10}μ+σ')text(1.8,-0.03,'\fontsize{10}μ+2σ')text(2.8,-0.03,'\fontsize{10}μ+3σ')

由圖可以看出盡管正態變量的取值范圍是(-∞,+∞)，但它的值落在(μ-3σ,μ+3σ)內幾乎是肯定的事。2.2統計分析的基本工具2.2.1數據集中度的度量1.樣本均值函數名稱：mean調用格式：m=mean(X)2.樣本中值函數名稱：median調用格式：m=median(X)3.幾何平均函數名稱：geomean調用格式：m=geomean(X)4.調和平均函數名稱：harmmean調用格式：m=harmmean(X)2.2.2數據差異性的度量1.樣本方差函數名稱：var調用格式：y=var(X)2.樣本標準差函數名稱：std調用格式：y=std(X)3.樣本極差函數名稱：range調用格式：y=range(X)2.2.3數據分布特征的度量1.樣本的經驗分位數函數名稱：prctile調用格式：y=prctile(X,p)2.樣本峰度函數名稱：kurtosis調用格式：y=kurtosis(X)3.樣本偏度函數名稱：skewness調用格式：y=skewness(X)2.2.4兩組數據線性相依程度的度量1.樣本協方差函數名稱：cov調用格式：y=cov(X)2.樣本相關系數函數名稱：corrcoef調用格式：y=corrcoef(X)2.3.1三大分布與分位數1.t分布和標準正態分布的比較（n=2,4,45）clear,clfX=linspace(-4,4,100);Y0=normpdf(X,0,1);Y1=tpdf(X,45);Y2=tpdf(X,4);Y3=tpdf(X,2);YY0=normpdf(0,0,1);plot(X,Y0,'.b',X,Y1,'-c',X,Y2,'-m',X,Y3,'-k',[0,0,],[0,YY0],':r')title('\fontsize{18}\fontname{華文新魏}不同自由度的t分布概率密度曲線')legend('N(0,1)','df:n=45','df:n=4','df:n=2')2.3統計分析的理論基礎2.α分位數的定義、求上側、下側、雙側α分位數介紹程序中幾個關鍵的MATLAB函數：（1）計算標準正態分布的0.05分位點。如上側α/2分位點的計算指令是:xalpha=norminv(0.975,0,1)下面幾個數值是標準正態分布的0.05分位點：下側分位點：-1.6449；上側分位點：1.6449；雙側分位點：-1.9600，1.9600。（2）生成樣本數據。如生成300個標準正態分布隨機數的計算指令是:data=normrnd(0,1,300,1)（3）繪制工序能力圖（繪制由分位點控制的密度曲線下的面積圖，用陰影表示，并計算樣本數據落入控制區域的概率，顯示在標題位置上）。計算指令是:capaplot(data,[xalpha,inf])程序如下：%alpha分位點示意圖（alpha=0.05）clearclfdata=normrnd(0,1,300,1);xalpha1=norminv(0.05,0,1);xalpha2=norminv(0.95,0,1);xalpha3=norminv(0.025,0,1);xalpha4=norminv(0.975,0,1);subplot(3,1,1)capaplot(data,[-inf,xalpha1]);axis([-3,3,0,0.45])subplot(3,1,2)capaplot(data,[xalpha2,inf]);axis([-3,3,0,0.45])subplot(3,1,3)capaplot(data,[-inf,xalpha3]);axis([-3,3,0,0.45])holdoncapaplot(data,[xalpha4,inf]);axis([-3,3,0,0.45])holdoff第三章統計估計3.1變量分布形態的估計3.1.1 頻率分布表

例3-1鋼材中的含硅量X是影響材料性能的一項重要因素。在煉鋼生產過程中，由于各種隨機因素的影響，各爐鋼的含硅量X是有差異的。對含硅量X概率分布的了解是有關鋼材料性能分析的重要依據。某煉鋼廠120爐正常生產的25MnSi鋼的含硅量（單位：%）如下：0.860.830.770.810.810.800.790.820.820.810.820.780.800.810.870.810.770.780.770.780.770.710.950.780.810.790.800.770.760.820.840.790.900.820.790.820.790.860.810.780.820.780.730.840.810.810.830.890.780.860.780.840.840.750.810.810.740.780.760.800.750.790.850.780.740.710.880.820.760.850.810.790.770.810.810.870.830.650.640.780.800.800.770.840.750.830.900.800.850.810.820.840.850.840.820.850.840.820.850.840.810.770.820.830.820.740.730.750.770.780.870.770.800.750.820.780.780.820.780.78編制25MnSi鋼含硅量數據的頻率分布表。頻率分布表的編制步驟：1.數據分組。（1）確定數據組個數。根據樣本容量n確定分組數k，推薦公式為k=1.87(n-1)2/5

。（2）計算極差。（3）確定組距。（4）確定各組端點。2.統計各組頻數。各組頻數就是數據落入各個小組中的個數，記為ni。3.計算頻率。（1）計算各組頻率。計算公式為fi=ni/n。（2）計算各組累積頻率。4.編制頻率分布表。程序如下：clearloadE:\matlab7\mnsi.txtk=ceil(1.87*(length(mnsi)-1)^0.4);[ni,ak]=hist(mnsi,k);fi=ni/length(mnsi);mfi=cumsum(fi);stats=[[1:k]',ak',ni',fi',mfi']表3-1120爐25MnSi鋼的含硅量數據頻率分布表3.1.2頻率直方圖繪制例3-1的頻率直方圖clearloadE:\matlab7\mnsi.txthist(mnsi)h=findobj(gca,'Type','patch');set(h,'FaceColor','y','EdgeColor','b')h=histfit(mnsi,13);%畫附正態參考曲線的直方圖，并提取圖形句柄hset(h(1),'FaceColor','c','EdgeColor','w')set(h(2),'Color','r')3.1.3經驗分布函數求出例3-1中25MnSi鋼含硅量的經驗分布函數clearloadE:\matlab7\mnsi.txt[h,stats]=cdfplot(mnsi)

由圖可以看出，樣本經驗分布函數圖像上升速度較快，均值與中值接近，圖像的S形狀均衡對稱，均值處函數值為0.5.這些特征表明，25MnSi鋼的含硅量可能服從均值為0.8026、標準差為0.0450的正態分布。接下來可以通過正態分布擬合檢驗進一步證實這種推測。3.2 變量分布參數的估計3.2.1 極大似然估計表3-2重復觀測樣本數據20124212021511311155例3-2對于表3-2中的數據，請利用MATLAB函數解答下列問題。（1）假設表3-2中的數據為X~P(λ)的重復觀測樣本數據，計算X方差的極大似然估計值；（2）假設表3-2中的數據為X~N(μ,σ2)的重復觀測樣本數據，計算μ的極大似然估計值；（3）假設表3-2中的數據為X~N(μ,σ2)的重復觀測樣本數據，計算μ和σ的極大似然估計值；（4）假設表3-2中的數據為X~U(a,b)的重復觀測樣本數據，計算a的極大似然估計值；（5）假設表3-2中的數據為X~U(a,b)的重復觀測樣本數據，計算a和b的極大似然估計值；解：先