Mathematica科學(xué)計算與程序設(shè)計 課件 第6、7章 Mathematica概率計算、Mathematica程序設(shè)計

第六章Mathematica概率計算6.1概率分布6.2隨機(jī)變量數(shù)字特征6.3事件概率本章介紹Mathematica概率論方面的應(yīng)用，主要包括概率密度函數(shù)、分布函數(shù)、隨機(jī)變量數(shù)字特征、概率計算、參數(shù)估計、分布檢驗(yàn)和回歸分析等。概率論是重要的基礎(chǔ)數(shù)學(xué)理論，廣泛應(yīng)用于工程科學(xué)和經(jīng)濟(jì)學(xué)中，也是量子計算的基礎(chǔ)數(shù)學(xué)理論之一。Mathematica集成了現(xiàn)有概率論與數(shù)理統(tǒng)計方面的全部算法實(shí)現(xiàn)函數(shù)，可以幫助學(xué)生快速掌握概率論的基礎(chǔ)理論，并且為學(xué)生在實(shí)際工程中應(yīng)用概率論奠定了算法基礎(chǔ)。第六章Mathematica概率計算6.4分布參數(shù)估計6.5線性回歸分析6.6蒙特卡羅實(shí)驗(yàn)本章介紹Mathematica概率論方面的應(yīng)用，主要包括概率密度函數(shù)、分布函數(shù)、隨機(jī)變量數(shù)字特征、概率計算、參數(shù)估計、分布檢驗(yàn)和回歸分析等。概率論是重要的基礎(chǔ)數(shù)學(xué)理論，廣泛應(yīng)用于工程科學(xué)和經(jīng)濟(jì)學(xué)中，也是量子計算的基礎(chǔ)數(shù)學(xué)理論之一。Mathematica集成了現(xiàn)有概率論與數(shù)理統(tǒng)計方面的全部算法實(shí)現(xiàn)函數(shù)，可以幫助學(xué)生快速掌握概率論的基礎(chǔ)理論，并且為學(xué)生在實(shí)際工程中應(yīng)用概率論奠定了算法基礎(chǔ)。6.1概率分布FirstImpressionMaurisquamdolor,cursusatportaet,luctusegetpurus.Nunctemporluctusinterdum.概率是對隨機(jī)事件發(fā)生可能性的量度。對于連續(xù)型隨機(jī)變量，概率密度函數(shù)為描述隨機(jī)變量在某個點(diǎn)附近取值的可能性的函數(shù)，隨機(jī)變量的值落在某個區(qū)間內(nèi)的概率為其概率密度函數(shù)在該區(qū)間上的積分。對于離散型隨機(jī)變量，直接考察隨機(jī)變量在離散點(diǎn)上的概率。在Mathematica中，使用PDF函數(shù)（概率密度函數(shù)）表征隨機(jī)變量的概率分布情況，PDF函數(shù)的基本語法為：FirstImpressionMaurisquamdolor,cursusatportaet,luctusegetpurus.Nunctemporluctusinterdum.（1）PDF[概率分布,隨機(jī)變量x]得到關(guān)于x的“概率分布”密度函數(shù)（2）PDF[概率分布,{x1,x2,…}]得到聯(lián)合“概率分布”密度函數(shù)6.1概率分布FirstImpressionMaurisquamdolor,cursusatportaet,luctusegetpurus.Nunctemporluctusinterdum.這里以正態(tài)分布和泊松分布為例，介紹概率密度函數(shù)的形態(tài)。其中，正態(tài)分布的PDF函數(shù)如圖6-1所示，泊松分布的PDF函數(shù)如圖6-2所示。在圖6-1中，“In[1]”得到標(biāo)準(zhǔn)正態(tài)分布的概率密度函數(shù)，如“Out[1]”所示；然后，“In[2]”繪制了該函數(shù)的圖形，如“Out[2]”所示。由圖6-1可知，正態(tài)分布的概率密度函數(shù)是典型的鐘形曲線，曲線關(guān)于y軸對稱，向x軸的兩側(cè)快速衰減。正態(tài)分布又稱高斯分布，是最重要的概率密度分布類型。FirstImpressionMaurisquamdolor,cursusatportaet,luctusegetpurus.Nunctemporluctusinterdum.6.1概率分布FirstImpressionMaurisquamdolor,cursusatportaet,luctusegetpurus.Nunctemporluctusinterdum.在圖6-2中，“In[3]”得到了均值為8的泊松分布的概率密度函數(shù)，如“Out[3]”所示；然后，“In[4]”繪制了其函數(shù)圖形，如“Out[4]”所示。從泊松分布的概率密度函數(shù)中可顯式獲得其隨機(jī)變量的均值，泊松分布在排隊(duì)論中有廣泛的應(yīng)用。FirstImpressionMaurisquamdolor,cursusatportaet,luctusegetpurus.Nunctemporluctusinterdum.在圖6-3中，“In[5]”計算了參數(shù)為5的指數(shù)分布的概率密度函數(shù)，如“Out[5]”所示；“In[8]”繪制了該指數(shù)分布的概率密度函數(shù)的圖形，如“Out[8]”所示。“In[9]”計算了參數(shù)為{0,4}的均勻分布的概率密度函數(shù)，如“Out[9]”所示；“In[10]”繪制了該均勻分布的概率密度函數(shù)的圖形，如“Out[10]”所示。在Mathematica中，計算隨機(jī)變量的概率分布函數(shù)使用CDF函數(shù)，其語法為：CDF[概率分布,隨機(jī)變量]下面以指數(shù)分布和均勻分布為例，介紹CDF函數(shù)的典型用法。指數(shù)分布（參數(shù)為5）的概率分布函數(shù)如圖6-4所示。在圖6-4中，繪制了指數(shù)分布（參數(shù)為5）的概率密度函數(shù)和概率分布函數(shù)?！癐n[11]”得到了參數(shù)為5的指數(shù)分布的概率密度函數(shù)pdf1，“In[12]”得到了參數(shù)為5的指數(shù)分布的累積分布函數(shù)cdf1；然后，“In[14]”中繪制了pdf1和cdf1的圖形，如“Out[14]”所示。對于連續(xù)型隨機(jī)變量，累積分布函數(shù)上的每個點(diǎn)x的值表示概率，這里X表示隨機(jī)變量。均勻分布（參數(shù)為{0,4}）的概率分布函數(shù)如圖6-5所示在圖6-5中，繪制了均勻分布（參數(shù)設(shè)為{0,4}）的概率密度函數(shù)和概率分布函數(shù)。“In[15]”得到了均勻分布（參數(shù)為{0,4}）的概率密度函數(shù)pdf2，“In[16]”得到了均勻分布（參數(shù)為{0,4}）的累積分布函數(shù)cdf2；然后，“In[18]”中繪制了pdf2和cdf2的圖形，如“Out[18]”所示。通過查看概率密度函數(shù)和概率分布函數(shù)的圖形，可以直觀地理解隨機(jī)變量的概率取值規(guī)律。6.2隨機(jī)變量數(shù)字特征FirstImpressionMaurisquamdolor,cursusatportaet,luctusegetpurus.Nunctemporluctusinterdum.給定隨機(jī)變量的概率分布，可計算服從該分布的隨機(jī)變量數(shù)字特征，一些常用的隨機(jī)變量數(shù)字特征包括：FirstImpressionMaurisquamdolor,cursusatportaet,luctusegetpurus.Nunctemporluctusinterdum.（1）均值Mean[概率分布]返回服從“概率分布”的隨機(jī)變量的平均值（2）方差Variance[概率分布]返回服從“概率分布”的隨機(jī)變量的方差（3）標(biāo)準(zhǔn)差StandardDeviation[概率分布]返回服從“概率分布”的隨機(jī)變量的標(biāo)準(zhǔn)差上述計算隨機(jī)變量的均值、方差和標(biāo)準(zhǔn)差的典型用法實(shí)例如圖6-6所示在圖6-6中，“In[22]”為均值為5的泊松分布；“In[23]”為參數(shù)為5的指數(shù)分布；“In[24]”和“In[25]”分別計算這兩個分布的均值，其結(jié)果如“Out[24]”和“Out[25]”所示；“In[26]”和“In[27]”分別計算這兩個分布的方差，其結(jié)果如“Out[26]”和“Out[27]”所示；“In[28]”和“In[29]”分別計算這兩個分布的標(biāo)準(zhǔn)差，其結(jié)果如“Out[28]”和“Out[29]”所示。6.2隨機(jī)變量數(shù)字特征FirstImpressionMaurisquamdolor,cursusatportaet,luctusegetpurus.Nunctemporluctusinterdum.除了均值、方差和標(biāo)準(zhǔn)差外，常用的隨機(jī)變量特征量還有：FirstImpressionMaurisquamdolor,cursusatportaet,luctusegetpurus.Nunctemporluctusinterdum.（1）中位數(shù)Median[概率分布]返回服從“概率分布”的隨機(jī)變量的中位數(shù)m，滿足P(X≤m)=1/2。（2）期望值Expectation[含隨機(jī)變量x的表達(dá)式，x服從的概率分布]計算當(dāng)x服從某一“概率分布”時，“含隨機(jī)變量x的表達(dá)式”的期望值。如果“表達(dá)式”中含有多個隨機(jī)變量，要求這些隨機(jī)變量是相互獨(dú)立的。計算中位數(shù)與期望值的典型實(shí)例如圖6-7所示在圖6-7中，“po”和“ex”分別為圖6-6中均值為5的泊松分布和參數(shù)為5的指數(shù)分布?！癐n[30]”和“In[31]”分別計算po和ex的中位數(shù)，其結(jié)果如“Out[30]”和“Out[31]”所示；“In[38]”計算x服從均值為5的泊松分布時x2+1的期望值，如“Out[38]”所示；“In[39]”計算x服從參數(shù)為5的指數(shù)分布時x2+1的期望值，如“Out[39]”所示。在“In[38]”中，“”表示變量服從某個概率分布的符號，其輸入方式為“Esc鍵+dist+Esc鍵”?！癐n[40]”計算在0至10上均勻分布的隨機(jī)變量的中位數(shù)，其結(jié)果如“Out[40]”所示；“In[41]”計算當(dāng)x滿足0至10上的均勻分布時x2+1的期望值，其結(jié)果如“Out[41]”所示。6.3事件概率FirstImpressionMaurisquamdolor,cursusatportaet,luctusegetpurus.Nunctemporluctusinterdum.除了均值、方差和標(biāo)準(zhǔn)差外，常用的隨機(jī)變量特征量還有：FirstImpressionMaurisquamdolor,cursusatportaet,luctusegetpurus.Nunctemporluctusinterdum.計算事件概率的典型實(shí)例如圖6-8所示6.3事件概率FirstImpressionMaurisquamdolor,cursusatportaet,luctusegetpurus.Nunctemporluctusinterdum.FirstImpressionMaurisquamdolor,cursusatportaet,luctusegetpurus.Nunctemporluctusinterdum.計算事件概率的典型實(shí)例如圖6-8所示在圖6-8中，“In[42]”和“In[44]”分別定義了標(biāo)準(zhǔn)正態(tài)分布和0至10上的均勻分布?！癐n[45]”計算當(dāng)x服從標(biāo)準(zhǔn)正態(tài)分布時的概率P(x≥0)，其結(jié)果如“Out[45]”所示；“In[46]”計算P(-2≤x≤2)，其結(jié)果如“Out[46]”所示，其中Erf為誤差函數(shù)；“In[47]”計算“Out[46]”的近似值，如“Out[47]”所示?！癐n[48]”計算當(dāng)隨機(jī)變量x服從0至10上的均勻分布時的概率P(x≤3)，如“Out[48]”所示?！癐n[54]”計算當(dāng)隨機(jī)變量x服從標(biāo)準(zhǔn)正態(tài)分布且y服從0至10上的均勻分布時的概率P(x+y2≤3)，其結(jié)果如“Out[54]”所示；“In[55]”計算得到“Out[54]”的近似值，如“Out[55]”所示。6.4分布參數(shù)估計FirstImpressionMaurisquamdolor,cursusatportaet,luctusegetpurus.Nunctemporluctusinterdum.FirstImpressionMaurisquamdolor,cursusatportaet,luctusegetpurus.Nunctemporluctusinterdum.在Mathematica中，給定一組數(shù)據(jù)，借助于DistributionFitTest函數(shù)檢驗(yàn)這組數(shù)據(jù)是否服從某一概率分布，DistributionFitTest函數(shù)的基本語法為“DistributionFitTest[數(shù)據(jù),概率分布]”，測試數(shù)據(jù)是否服從指定的“概率分布”，返回一個p值，當(dāng)p值大于0.05時，認(rèn)為測試數(shù)據(jù)集服從指定的“概率分布”。DistributionFitTest函數(shù)的典型用法實(shí)例如圖6-9所示。6.4分布參數(shù)估計FirstImpressionMaurisquamdolor,cursusatportaet,luctusegetpurus.Nunctemporluctusinterdum.FirstImpressionMaurisquamdolor,cursusatportaet,luctusegetpurus.Nunctemporluctusinterdum.大于0.05，故可以認(rèn)為dat2服從-2至2的均勻分布。工程中獲得的數(shù)據(jù)往往由經(jīng)驗(yàn)可推知其概率分布類型，這時可以使用函數(shù)EstimatedDistribution估計分布的參數(shù)，該函數(shù)的典型語法為“EstimatedDistribution[測試數(shù)據(jù),帶未知參數(shù)的分布]”，其實(shí)例如圖6-10所示。6.5線性回歸分析“Out[104]”所示；“In[105]”將回歸分析結(jié)果轉(zhuǎn)化為多項(xiàng)式的形式，得到y(tǒng)=2.84624+0.93488x，如“Out[105]”所示；“In[106]”繪制dat的散點(diǎn)圖g1，“In[107]”繪制回歸線g2，“In[108]”將g1和g2放在一起顯示，如“Out[108]”所示。由“Out[108]”所示的圖形可知，回歸線的擬合程度較好。然后，“In[109]”做y=a0+a1x2回歸，得到系數(shù)a0=3.69005，a1=0.163184，如“Out[109]”所示；“In[110]”將回歸函數(shù)轉(zhuǎn)化為多項(xiàng)式的形式，得到y(tǒng)=3.69005+0.163184x2，如“Out[110]”已知一組數(shù)據(jù)（二維嵌套列表的形式），尋找這組數(shù)據(jù)中存在的線性關(guān)系，可以借助于線性回歸分析方法實(shí)現(xiàn)。在Mathematica中，線性回歸分析函數(shù)為LinearModelFit，其基本語法為“LinearModelFit[{{x11,x12,…,x1n,y1},{x21,x22,…,x2n,y2},…},{f1,f2,…,fn},{x1,x2,…,xn}]”，表示建立回歸模型y=a0+a1f1(x1)+a2f2(x2)+….+anfn(xn)。對于一元線性回歸而言，其典型實(shí)例如圖6-11所示。在圖6-11中，“In[103]”引入一組數(shù)據(jù)dat，“In[104]”做y=a0+a1x線性回歸分析，得到系數(shù)a0=2.84624，a1=0.93488，如6.5線性回歸分析所示；“In[114]”繪制回歸線g3，“In[115]”將回歸線g3與圖形g1繪制在一起顯示，如“Out[115]”所示。由“Out[115]”可直觀地看出，回歸線的擬合程度較好。對于多元線性回歸而言，其典型實(shí)例如圖6-12所示。在圖6-12中，“In[121]”輸入一組數(shù)據(jù)dat，如“Out[121]”所示；“In[122]”做y=a0+a1x+a2y線性回歸分析，得到回歸系數(shù)a0=3.10189，a1=1.7465，a1=5.81191，如“Out[122]”所示；“In[123]”計算在點(diǎn)(0.1,0.02)處的回歸線的值，其結(jié)果如“Out[123]”所示。6.5線性回歸分析所示；“In[114]”繪制回歸線g3，“In[115]”將回歸線g3與圖形g1繪制在一起顯示，如“Out[115]”所示。由“Out[115]”可直觀地看出，回歸線的擬合程度較好。對于多元線性回歸而言，其典型實(shí)例如圖6-12所示。在圖6-12中，“In[121]”輸入一組數(shù)據(jù)dat，如“Out[121]”所示；“In[122]”做y=a0+a1x+a2y線性回歸分析，得到回歸系數(shù)a0=3.10189，a1=1.7465，a1=5.81191，如“Out[122]”所示；“In[123]”計算在點(diǎn)(0.1,0.02)處的回歸線的值，其結(jié)果如“Out[123]”所示?！癐n[124]”獲得回歸殘差值，即真實(shí)數(shù)據(jù)與回歸線相應(yīng)點(diǎn)處的值的差值，其結(jié)果如“Out[124]”所示?！癐n[125]”做y=a0+a1x2+a2y線性回歸分析，得到回歸系數(shù)為a0=3.3774，a1=7.06987，a1=1.15277，如“Out[125]”所示；“In[126]”獲得回歸殘差值，其結(jié)果如“Out[216]”所示；然后，“In[129]”繪制了兩次擬合的殘差對比圖，如“Out[129]”所示。6.6蒙特卡羅實(shí)驗(yàn)FirstImpressionMaurisquamdolor,cursusatportaet,luctusegetpurus.Nunctemporluctusinterdum.FirstImpressionMaurisquamdolor,cursusatportaet,luctusegetpurus.Nunctemporluctusinterdum.蒙特卡羅實(shí)驗(yàn)基于概率統(tǒng)計理論，使用計算機(jī)模擬生成各種概率分布的偽隨機(jī)數(shù)作為數(shù)學(xué)模型的輸入，從而求得模型的近似解或統(tǒng)計量。隨著計算機(jī)技術(shù)的發(fā)展，蒙特卡羅實(shí)驗(yàn)已成為現(xiàn)代測試技術(shù)中重要的實(shí)驗(yàn)手段。最早的蒙特卡羅實(shí)驗(yàn)是計算圓周率的概率論方法，設(shè)一個邊長為2的正方形和它的內(nèi)切單位圓，其面積之比為π/4，若隨機(jī)地向該正方形內(nèi)投入大量的“針”（設(shè)共n個），則位于單位圓內(nèi)部的針的數(shù)量n與總的針數(shù)量n之比m/n應(yīng)為π/4。6.6蒙特卡羅實(shí)驗(yàn)FirstImpressionMaurisquamdolor,cursusatportaet,luctusegetpurus.Nunctemporluctusinterdum.FirstImpressionMaurisquamdolor,cursusatportaet,luctusegetpurus.Nunctemporluctusinterdum.用Mathematica實(shí)現(xiàn)計算圓周率的方法如圖6-13所示。其中，函數(shù)RandomVariate的語法為“RandomVariate[概率分布,n]”，用于產(chǎn)生n個滿足指定“概率分布”的隨機(jī)變量；函數(shù)Norm的語法為“Norm[向量]”，給出指定“向量”的模；If函數(shù)的用法請參考第7.2.1節(jié)。在圖6-13中，輸入n的值為100000；然后，調(diào)用RandomVariate函數(shù)生成滿足二元均勻分布的n個點(diǎn)，這里的UniformDistribution[{{-1,1},{-1,1}}]為區(qū)域{{-1,1},{-1,1}}內(nèi)二元均勻分布。接著，令m為0；借助于6.6蒙特卡羅實(shí)驗(yàn)FirstImpressionMaurisquamdolor,cursusatportaet,luctusegetpurus.Nunctemporluctusinterdum.FirstImpressionMaurisquamdolor,cursusatportaet,luctusegetpurus.Nunctemporluctusinterdum.Table函數(shù)，局部變量i從1至n，當(dāng)點(diǎn)needle[[i]]的長度（模）小于1時，認(rèn)為該點(diǎn)落在單位圓內(nèi)部，則m自增1，從而m為落入單位圓內(nèi)部的點(diǎn)的數(shù)量。最后，圓周率的值pi=4m/n，計算結(jié)果如“Out[43]”所示，為3.13548。通過上述實(shí)例可知，蒙特卡羅實(shí)驗(yàn)的特點(diǎn)在于需要大量的偽隨機(jī)數(shù)。上述計算圓周率的實(shí)例中，生成了100000個點(diǎn)，即100000對偽隨機(jī)數(shù)，但是計算結(jié)果與理論值（3.14159）間的誤差仍然比較大。下面考慮一種“公平”的賭輸贏問題設(shè)莊家有幣a=100000，玩家有幣b=1000，通過擲骰子進(jìn)行賭輸贏。骰子有6個面，依次鑲嵌著1至6個點(diǎn)。每次擲骰子將得到一個1至6間的一個數(shù)字點(diǎn)數(shù)。設(shè)玩家總是賭大，即點(diǎn)數(shù)為4、5或6時，玩家贏；莊家總是賭小，即點(diǎn)數(shù)為1、2或3時，莊家贏。假設(shè)比賽是公平的，即每次擲骰子得到的點(diǎn)數(shù)不受人為控制。設(shè)每局賭注均為m=10個幣，即贏了得10個幣，輸了失10個幣。設(shè)一共賭n=500局。模擬上述“公平”的賭輸贏問題，如圖6-14所示，dealer表示莊家，player表示玩家，可見玩家是必輸?shù)?。下面考慮一種“公平”的賭輸贏問題在圖6-14中，a=1000表示玩家的持幣數(shù)量，b=100000表示莊家的持幣數(shù)量，n=500表示賭局的次數(shù)；m=10表示每局輸贏的幣數(shù)；dice保存隨機(jī)生成的n個1至6間的隨機(jī)整數(shù)。這里的DiscreteUniformDistribution[{min,max}]表示生成整數(shù)min至整數(shù)max間的隨機(jī)整數(shù)。在Table函數(shù)中，統(tǒng)計dice列表中每個元素是否大于3，對于大于3的情況，玩家a得10個幣，莊家少10個幣；否則，玩家少10個幣，莊家得10個幣。最后，打印莊家dealer的幣數(shù)和玩家player的幣數(shù)。重復(fù)6-14所示的實(shí)驗(yàn)，發(fā)現(xiàn)即使公平情況下的賭輸贏，玩家也是必輸?shù)?。在公鑰密碼RSA中，大素數(shù)是其安全的首要保證。實(shí)際上，尋找大的素數(shù)是非常困難的事情，常用Miller-Rabin概率方法測試一個數(shù)是否可能為素數(shù)，該方法基于如下的定理（摘自C.Paar和J.Pelzl的《深入淺出密碼學(xué)》）：由圖6-15可知，734293417為一個合數(shù)（如“In[158]”和“Out[158]”所示），而573480473為一個素數(shù)（如“In[157]”和“Out[157]”所示）。下面使用Miller-Rabin方法測試這兩個數(shù)的素性，如圖6-16和圖6-17所示。圖6-16中的代碼如下所示：1p=573480473;2u=0;3If[Mod[p-1,2]==1,Print[p,"isnotaprime."],4r=p-1;5While[Mod[r,2]==0,u++;r=r/2];67len=5;8a=RandomVariate[DiscreteUniformDistribution[{2,p-2}],len];910For[i=1,i<=len,i++,11pq=0;12t=a[[i]]^r;If[Mod[t,p]!=1,pq=pq+1];13For[j=0,j<u,j++,14If[Mod[Power[t,2^j],p]!=p-圖6-16中的代碼如下所示：1,pq++]];15If[pq==u+1,Print[p,"isnotaprime."];Break[]]16];17If[pq<=u,Print[p,"isaprime."]]18]上述代碼中，第1行將p設(shè)為大數(shù)734293417；第2行變量u賦為0。第3行至第18行為一條If語句。第3行判斷p-1除以2的余數(shù)是否為1，如果為1，表明p-1為奇數(shù)，即p為偶數(shù)，所以，輸出p不是一個素數(shù)；否則，執(zhí)行第4行至第17行。第4行中令r等于p-1；第5行循環(huán)執(zhí)行r=r/2，直接r為奇數(shù)。這里計算得到r和u即為式（6.1）中的r和u。第7行令len為5，len表示隨機(jī)從{2,3,…,p-2}中取到的a的個數(shù)，len稱為安全參數(shù)，len的值越大，計算量越大，判定結(jié)果越接近真實(shí)情況。第8行調(diào)用RandomVariate函數(shù)生成2至p-2間均勻分布的整數(shù)序列，序列長度為len，即序列包含len個整數(shù)。圖6-16中的代碼如下所示：1,pq++]];15If[pq==u+1,Print[p,"isnotaprime."];Break[]]16];17If[pq<=u,Print[p,"isaprime."]]18]第10行至第16行為For循環(huán)語句。變量i從1循環(huán)遞增至len，每次循環(huán)中針對一個a[[i]]的值，執(zhí)行如下的操作：（1）第11行pq賦為0；（2）第12行計算t等于a[[i]]的r次冪，并判斷t除以p的余數(shù)是否為1，如果不為1，則pq自增1；（3）第13行至第14行為一個For循環(huán)語句，執(zhí)行式（6.2）的算法。第13行循環(huán)變量j從0遞增到u-1；針對每一個j，第14行判斷t的2j冪除以p的余數(shù)是否為p-1，如果不為p-1，則pq自增1。（4）第15行，判斷pq的值是否為u+1，如果pq等于u+1，則打印p不是素數(shù)，并調(diào)用Break函數(shù)終止程序。第17行，判斷pq的值是否小于等于u，如果pq小于等于u，則打印p是一個素數(shù)。在圖6-17中，p輸入大數(shù)573480473，執(zhí)行程序可得到結(jié)果“573480473isaprime.”，表示573480473為一個素數(shù)。本章討論了常用的概率密度函數(shù)及其分布函數(shù)；然后，基于概率密度函數(shù)闡述了隨機(jī)變量常用數(shù)字特征的計算方法，主要包括均值、方差、標(biāo)準(zhǔn)差、期望和中位數(shù)等；接著，分析了隨機(jī)事件的概率計算方法；之后，討論了測試數(shù)據(jù)的概率分布和分布參數(shù)估計方法；最后，討論了線性回歸方法及其典型實(shí)例。概率論已經(jīng)成為現(xiàn)代科學(xué)技術(shù)的重要數(shù)學(xué)工具，在本章的學(xué)習(xí)基礎(chǔ)上，建議借助于Mathematica參考資料中心，進(jìn)一步深入學(xué)習(xí)概率論相關(guān)的函數(shù)工作原理與應(yīng)用方法。06本章小結(jié)第七章Mathematica程序設(shè)計7.1函數(shù)7.2控制語句7.3純函數(shù)Mathematica軟件平臺上用于科學(xué)計算的指令和語句統(tǒng)稱為Wolfram語言，Wolfram語言屬于高級計算機(jī)語言，也是全球頂級的科學(xué)計算語言，是物理學(xué)等自然科學(xué)研究和數(shù)學(xué)與密碼學(xué)等信息科學(xué)研究的利器。本章將詳細(xì)介紹使用Wolfram語言進(jìn)行程序設(shè)計的方法，主要內(nèi)容包括函數(shù)定義與應(yīng)用方法、常用程序設(shè)計控制語句、模塊定義與調(diào)用方法和一些程序設(shè)計實(shí)例。第七章Mathematica程序設(shè)計7.4模塊7.5程序設(shè)計實(shí)例7.6本章小結(jié)Mathematica軟件平臺上用于科學(xué)計算的指令和語句統(tǒng)稱為Wolfram語言，Wolfram語言屬于高級計算機(jī)語言，也是全球頂級的科學(xué)計算語言，是物理學(xué)等自然科學(xué)研究和數(shù)學(xué)與密碼學(xué)等信息科學(xué)研究的利器。本章將詳細(xì)介紹使用Wolfram語言進(jìn)行程序設(shè)計的方法，主要內(nèi)容包括函數(shù)定義與應(yīng)用方法、常用程序設(shè)計控制語句、模塊定義與調(diào)用方法和一些程序設(shè)計實(shí)例。7.1函數(shù)Wolfram語言中的函數(shù)包括內(nèi)置函數(shù)、包函數(shù)和自定義函數(shù)三種，其中，內(nèi)置函數(shù)名以大寫字母開頭，可以直接調(diào)用；包函數(shù)是指存儲在Wolfram資源庫中的線上函數(shù)，使用時需動態(tài)地裝入函數(shù)所在的包文件；自定義函數(shù)為用戶編寫的函數(shù)，可用于擴(kuò)展內(nèi)置函數(shù)的功能。Wolfram語言中，函數(shù)的調(diào)用方法非常靈活，設(shè)函數(shù)名為f，其參數(shù)為x，則基本的調(diào)用方法為f[x]，即使用中括號將參數(shù)包括其中；另一種更簡潔的方法為f@x或x//f，均等價于f[x]。如果多個函數(shù)f、g和h嵌套調(diào)用，基本用法為f[g[h[x]]]，也可使用f@g@h@x或x//h//g//f。函數(shù)調(diào)用的另一種常見表示為“/@”，例如，f/@x，表示函數(shù)f作用于列表x的每一個元素，因此，f/@{1,2,3}等價于{f[1],f[2],f[3]}。7.1函數(shù)Wolfram語言內(nèi)置了常用的科學(xué)計算相關(guān)的函數(shù)。這里以偽隨機(jī)數(shù)相關(guān)的函數(shù)為例，介紹內(nèi)置函數(shù)的用法。偽隨機(jī)數(shù)相關(guān)的常用函數(shù)及其基本用法如表7.1所示。7.1.1內(nèi)置函數(shù)7.1.1內(nèi)置函數(shù)7.1.1內(nèi)置函數(shù)7.1.1內(nèi)置函數(shù)例7.1偽隨機(jī)數(shù)函數(shù)應(yīng)用舉例（1）生成長度為10的0或1偽隨機(jī)序列，并統(tǒng)計其中1的個數(shù)，見表7.2序號1。（2）生成長度為12的10至99間的偽隨機(jī)整數(shù)序列，見表7.2序號2。（3）生成3×4的偽隨機(jī)實(shí)數(shù)矩陣，每個元素的取值在1和2之間，見表7.2序號3。（4）生成位于直角坐標(biāo)系中頂點(diǎn)為(0,0)、(1,0)、(0,3)和(1,3)的矩形內(nèi)的偽隨機(jī)復(fù)數(shù)序列，序列長度為20，見表7.2序號4。（5）隨機(jī)生成一個小于100的素數(shù)，見表7.2序號5。（6）從1至100中隨機(jī)選擇10個數(shù)（有重復(fù)采樣），見表7.2序號6。（7）從1至100中隨機(jī)選擇10個數(shù)（不重復(fù)采樣），見表7.2序號7。（8）生成服從正態(tài)分布的10個偽隨機(jī)數(shù)序列，見表7.2序號8。（9）給定偽隨機(jī)數(shù)種子為20200705，生成長度為10的偽隨機(jī)整數(shù)序列，每個元素在10至99間取值，然后，從小至大進(jìn)行排序。見表7.2序號9。表7.2偽隨機(jī)數(shù)函數(shù)典型應(yīng)用實(shí)例表7.2中序號1至8的執(zhí)行結(jié)果與時間有關(guān)，讀者每次執(zhí)行時得到的結(jié)果是不同的，所以，這里不給出這些語句的執(zhí)行結(jié)果。但是表7.2中序號9的執(zhí)行結(jié)果是不隨時間變化的，執(zhí)行后，變量a的結(jié)果一定為列表{81,67,94,51,22,95,65,17,92,41}，而排序后的結(jié)果為{17,22,41,51,65,67,81,92,94,95}。7.1函數(shù)7.1.1內(nèi)置函數(shù)7.1函數(shù)Mathematica不但是全球最先進(jìn)的科學(xué)計算軟件，而且也是全球最先進(jìn)的科學(xué)數(shù)據(jù)庫，實(shí)時地記錄著全球乃至可觀測宇宙部分的動態(tài)信息，供科研人員參考和使用。Mathematica線上資源特別豐富，其中包含了大量的軟件包，每個軟件包以文件的形式存儲在服務(wù)器上。這些軟件包中的函數(shù)稱為包函數(shù)，這里以有限域算術(shù)包為例介紹包函數(shù)的用法。在線上資源庫中，有限域算術(shù)包的名稱為FiniteField，在計算機(jī)聯(lián)網(wǎng)的情況下，可在Notebook中調(diào)入有限域算術(shù)包，即輸入

<<FiniteFields`或Needs["FiniteFields`"]然后，可以使用有限域算術(shù)包中的全部函數(shù)。這里重點(diǎn)介紹有限域算術(shù)中創(chuàng)建有限域?qū)ο蟮暮瘮?shù)，如表7.3所示。7.1.2包函數(shù)7.1函數(shù)函數(shù)是Mathematica程序設(shè)計的基本要素。在Mathematica中，程序由函數(shù)組成，而函數(shù)一般以模塊的形式供用戶調(diào)用。模塊將在第7.3節(jié)中介紹，這里重點(diǎn)介紹自定義函數(shù)及其用法。由于Mathematica系統(tǒng)內(nèi)置函數(shù)和包函數(shù)均以大寫字母開頭，因此，自定義函數(shù)盡可能以小寫字母開頭。自定義函數(shù)的格式為：函數(shù)名[參數(shù)1_,參數(shù)2_,…,參數(shù)n_]=表達(dá)式。需要注意的是，在參數(shù)表中，每個參數(shù)后有一個下劃線“_”，而表達(dá)式中的參數(shù)不再具有下劃線。這里的參數(shù)以所謂的符號或變量的形式出現(xiàn)，其對應(yīng)的實(shí)參可以為數(shù)值、列表或函數(shù)名。7.1.3自定義函數(shù)7.2控制語句程序語句的執(zhí)行只有三種方式，即順序、分支和循環(huán)。在Mathematica的Notebook中輸入的計算語句按照順序執(zhí)行的方式執(zhí)行。這里重點(diǎn)介紹分支和循環(huán)控制語句。7.2.1分支控制（1）If[條件，語句組1，語句組2，語句組3]如果條件為邏輯真，則執(zhí)行語句組1；如果條件為邏輯假，則執(zhí)行語句組2；如果條件在邏輯上非真非假，則執(zhí)行語句組3。每個語句組可以包含多個語句，語句間用分號“;”分隔。If語句的其他形式簡化形式有：(i)If[條件，語句組1]當(dāng)條件為邏輯真時，執(zhí)行語句組1；否則，無操作。(ii)If[條件，，語句組2]當(dāng)條件為邏輯假時，執(zhí)行語句組2；否則，無操作。注意，這里中間的逗號不可缺少。(iii)If[條件，語句組1，語句組2]當(dāng)條件為邏輯真時，執(zhí)行語句組1；當(dāng)條件為邏輯假時，執(zhí)行語句組2。7.2控制語句程序語句的執(zhí)行只有三種方式，即順序、分支和循環(huán)。在Mathematica的Notebook中輸入的計算語句按照順序執(zhí)行的方式執(zhí)行。這里重點(diǎn)介紹分支和循環(huán)控制語句。7.2.1分支控制（2）Switch[表達(dá)式，情況1，語句組1，情況2，語句組2，……]如果表達(dá)式為情況1，則執(zhí)行語句組1；如果表達(dá)式為情況2，則執(zhí)行語句組2。每個語句組可以有多個語句，使用分號“;”分隔。（3）Which[邏輯表達(dá)式1，值1，邏輯表達(dá)式2，值2，……]Which語句依次計算各個邏輯表達(dá)式的值，返回第一個為真的邏輯表達(dá)式i對應(yīng)的值i。（4）條件控制符（“/;”）后面接一個邏輯表達(dá)式test當(dāng)邏輯表達(dá)式test為真時，條件控制起作用。現(xiàn)在，使用分支控制語句實(shí)現(xiàn)下述分段函數(shù)，7.2控制語句Mathematica提供了三種循環(huán)控制語句，即Do語句、While語句和For語句。這三種語句的語法如表7.10所示。7.2.2循環(huán)控制7.2控制語句Wolfram語言是一種高級程序設(shè)計語言，具有高級循環(huán)控制語句，如Table、Array、Nest和NestList等。在大多數(shù)情況下，While、Do和For循環(huán)常被這些高級循環(huán)控制語句替代。這些高級循環(huán)控制語句的語法如表7.12所示。7.2.3高級循環(huán)控制7.2.3高級循環(huán)控制7.2.3高級循環(huán)控制7.2.3高級循環(huán)控制7.3純函數(shù)純函數(shù)又稱純匿名函數(shù)。純函數(shù)的格式有兩種，其一，“Function[函數(shù)體]”；其二，“函數(shù)體&”。這里“函數(shù)體”中，#表示形式參數(shù)；當(dāng)有多個形式參數(shù)時，#1表示第一個形式參數(shù)，#2表示第二個形式參數(shù)，依此類推。在第一種情況下，還可以指定形式參數(shù)，如Function[形式參數(shù)列表，函數(shù)體]。下面列舉幾個典型語句，如表7.14所示。7.3純函數(shù)純函數(shù)一個重要的作用體現(xiàn)在表7.14序號7的應(yīng)用中。Wolfram語言的基本數(shù)據(jù)結(jié)構(gòu)為列表，而借助于純函數(shù)可以把函數(shù)作用施加到列表的各個元素。注意，有些函數(shù)具有Listable屬性（例如，Sin函數(shù)，使用“??Sin”可查看Sin函數(shù)的屬性），這類函數(shù)作用于列表時，將作用于列表的每個元素。例如，Sin[{0,Pi/2,Pi,3Pi/2,2Pi}]將返回{0,1,0,-1,0}。7.4模塊在Wolfram語言中，模塊是程序的基本單位，模塊對應(yīng)著程序設(shè)計實(shí)現(xiàn)的子功能，類似于C語言中的函數(shù)。常用的模塊有四種，即With、Block、Module和Compile。這里借助于Hénon映射生成偽隨機(jī)序列為例，介紹模塊的應(yīng)用方法。Hénon映射是一種二維離散混沌，其吸引子的方程為

其中，a=1.4，b=0.3。給定初始值x0和y0，生成指定長度的偽隨機(jī)序列，序列的每個元素為0至255間取值的整數(shù)。With語句的基本語法為With[初始化變量列表,語句組]7.4.1With模塊7.4.1With模塊其中，“語句組”可包含多個語句，各個語句間用分號“；”分隔?！罢Z句組”中使用的變量可以為全局變量，也可以為“初始化變量列表”中的變量。如果“語句組”中使用了“初始化變量列表”中的變量，這些變量的初始值直接代入語句組中。特別注意：“初始化變量列表”中的變量均為With語句內(nèi)部可見的局部變量。With語句的典型實(shí)例如表7.15所示。With語句的基本語法為With[初始化變量列表,語句組]7.4.1With模塊With語句的處理速度比第7.4.3節(jié)將要介紹的Module語句快。With語句的典型用法如表7.15的序號3所示。注意：在With模塊中出現(xiàn)的非“初始化變量列表”中的變量均為全局變量。借助于With語句使用Hénon映射生成偽隨機(jī)序列的程序如下：henon:={#[[2]],1.0-1.4#[[2]]^2+0.3#[[1]]}&;f[x0_,y0_,m_]:=With[{x1=x0,y1=y0,m1=m},dat1=NestList[henon,{x1,y1},m1];dat2=Flatten[dat1][[3;;-1;;2]];dat3=Mod[Floor[(2+dat2)*10^6],256]]f[0.3,0.23,30]上述代碼中，定義了兩個函數(shù)henon和f，其中，henon使用純函數(shù)定義；函數(shù)f具有三個參數(shù)x0、y0和m，x0和y0作為Hénon映射的迭代初值，m為生成的偽隨機(jī)序列的長度。調(diào)用f[0.3,0.23,30]的執(zhí)行結(jié)果如下：{240,4,204,52,14,106,206,79,232,189,113,50,255,186,203,163,254,170,182,123,147,207,33,81,174,249,115,33,211,113}需要特別注意的是，在函數(shù)f的定義中，dat1、dat2和dat3在函數(shù)f被調(diào)用后，將成為Notebook中的全局變量。但是大多數(shù)情況下，只希望函數(shù)返回值，并不產(chǎn)生任何全局變量。Block和Module模塊可實(shí)現(xiàn)這種變量局部化的要求。7.4模塊類似于With模塊，Block模塊可以直接使用全局變量，并可以生成全局變量。而全局變量的存在破壞了算法的模塊化，因此，Block模塊還實(shí)現(xiàn)了全部變量局部化方法。Block語句的語法為：

Block[局部變量列表,語句組]或Block[局部變量初始化列表,語句組]這里的“語句組”可包含多條語句，各條語句間用分號“；”分隔，語句組的最后一條語句的執(zhí)行結(jié)果為返回值。Block語句的典型用法實(shí)例如表7.16所示。7.4.2Block模塊借助于Block語句使用Hénon映射生成偽隨機(jī)序列的程序如下：henon:={#[[2]],1.0-1.4#[[2]]^2+0.3#[[1]]}&;f[x0_,y0_,m_]:=Block[{x1=x0,y1=y0,m1=m,dat1,dat2,dat3},dat1=NestList[henon,{x1,y1},m1];dat2=Flatten[dat1][[3;;-1;;2]];dat3=Mod[Floor[(2+dat2)*10^6],256]]f[0.3,0.23,30]上述代碼中，先定義了henon函數(shù)和f函數(shù)，在f函數(shù)中使用了Block語句，將Block語句中使用了變量均定義為局部變量，運(yùn)行結(jié)果由“語句組”的最后一條語句返回。執(zhí)行函數(shù)“f[0.3,0.23,30]”返回：{240,4,204,52,14,106,206,79,232,189,113,50,255,186,203,163,254,170,182,123,147,207,33,81,174,249,115,33,211,113}在使用Block模塊時，如果“語句組”中的表達(dá)式包含了局部變量，將使用局部變量定義的初始值，這種方式稱為動態(tài)計算“語句組”中的表達(dá)式（這是和Module模塊唯一的區(qū)別，下文的Module模塊是“靜態(tài)”處理的）。例如：y=x^2+2x+1Block[{x=a},x+y]返回1+3a+a^2。而Module[{x=a},x+y]返回1+a+2x+x^2。7.4.2Block模塊7.4模塊在絕大多數(shù)情況下，使用Module語句實(shí)現(xiàn)模塊的功能，最主要的原因，可能是因?yàn)椤癕odule”英文有“模塊”的含義。由7.4.2節(jié)可知，當(dāng)全部使用局部變量時，Block模塊與Module模塊完全通用，事實(shí)上，Wolfram語言中Block模塊的處理速度比Module模塊更快。Module語句的語法如下：

Module[局部變量列表,語句組]或Module[局部變量初始化列表,語句組]這里的“語句組”可包含多條語句，各條語句間用分號“；”分隔，語句組的最后一條語句的執(zhí)行結(jié)果為返回值。Module語句的典型用法實(shí)例如表7.17所示。7.4.3Module模塊由Module語句使用Hénon映射生成偽隨機(jī)序列的程序如下Clear["`*"]henon:={#[[2]],1.0-1.4#[[2]]^2+0.3#[[1]]}&;f[x0_,y0_,m_]:=Module[{x1=x0,y1=y0,m1=m,dat1,dat2,dat3},dat1=NestList[henon,{x1,y1},m1];dat2=Flatten[dat1][[3;;-1;;2]];dat3=Mod[Floor[(2+dat2)*10^6],256];{dat1[[2;;-1]],dat3}]{ps1,ps2}=f[0.7,0.13,3000]上述程序代碼首先定義了henon函數(shù)和f函數(shù)，f函數(shù)的參數(shù)x0和y0為Hénon映射的迭代初值，m設(shè)定產(chǎn)生的混沌偽隨機(jī)序列的長度。函數(shù)f的返回值為列表{dat1[[2;;-1]],dat3}，列表中有兩個元素dat1[[2;;-1]]和dat3，依次為Hénon映射的狀態(tài)值序列和偽隨機(jī)序列。執(zhí)行“{ps1,ps2}=f[0.7,0.13,3000]”獲得長度為3000的狀態(tài)序列和偽隨機(jī)序列。下面使用ps1借助于ListPlot[Reverse/@ps1,AxesLabel->{"x","y"},AspectRatio->Automatic,ImageSize->Large]繪制Hénon映射的相圖如下：建議在Module模塊中全部使用局部變量，并將Module模塊定義為函數(shù)的實(shí)現(xiàn)部分。由Module語句使用Hénon映射生成偽隨機(jī)序列的程序如下7.4模塊Mathematica軟件使用C語言實(shí)現(xiàn)的，用戶編寫的Mathematica程序是順序調(diào)用Mathematica系統(tǒng)函數(shù)實(shí)現(xiàn)的。一般地，用戶程序的執(zhí)行效率遠(yuǎn)遠(yuǎn)低于C語言可執(zhí)行程序。為了提高M(jìn)athematica用戶程序的執(zhí)行效率，可以使用Compile模塊，設(shè)計經(jīng)過編譯的用戶函數(shù)。這種編譯后的函數(shù)的執(zhí)行效率非常接近于C語言可執(zhí)行程序的效率，但是這類函數(shù)只能使用常規(guī)的數(shù)據(jù)類型，即整型（_Integer）、浮點(diǎn)型（_Real）、復(fù)數(shù)類型（_Complex）和邏輯變量（True或False），使用的列表必須為數(shù)值型數(shù)組?？梢?，雖然Compile模塊可以編譯為機(jī)器代碼，極大地提高了執(zhí)行速度，但是遠(yuǎn)沒有Module模塊靈活。7.4.4Compile模塊Compile模塊的典型語法為Compile[{{變量名1,變量類型1},{變量名2,變量類型2},…},語句組]或Compile[{{變量名1,變量類型1},{變量名2,變量類型2},…},語句組,屬性]或Compile[{{變量名1,變量類型1,變量1維數(shù)},{變量名2,變量類型2,變量2維數(shù)},…},語句組]上述語法中，“語句組”可以包含多條語句，各條語句間用分號“；”隔開。Compile模塊的典型用法實(shí)例如表7.18所示。由Complile語句使用Hénon映射生成偽隨機(jī)序列的程序如下：Clear["`*"]henon:={#[[2]],1.0-1.4#[[2]]^2+0.3#[[1]]}&;f=Compile[{{x0,_Real},{y0,_Real},{m,_Integer}},Module[{x1=x0,y1=y0,m1=m,dat1,dat2,dat3},dat1=NestList[henon,{x1,y1},m1];dat2=Flatten[dat1][[3;;-1;;2]];dat3=Mod[Floor[(2+dat2)*10^6],256]]]ps=f[0.3,0.23,30]返回值為：{240,4,204,52,14,106,206,79,232,189,113,50,255,186,203,163,254,170,182,123,147,207,33,81,174,249,115,33,211,113}。上述程序中在Compile模塊內(nèi)部嵌入了Module模塊。Compile模塊的最大用處在于可以借助于Compile模塊比較不同算法的運(yùn)算速度，Compile模塊本質(zhì)上是機(jī)器語言程序，可以準(zhǔn)確地反映算法的運(yùn)行速度。本節(jié)以兩個常用的對稱密碼算法（RC4和SM4）為例，介紹Mathematica程序設(shè)計方法。這里使用了Module模塊作為程序設(shè)計基本單元。對于加密算法而言，輸入為明文和密鑰，輸出為密文；對于解密算法而言，輸入為密文和密鑰，輸出為明文。RC4是一種常用的流密碼，而SM4是我國的一項(xiàng)文本加密標(biāo)準(zhǔn)。7.5程序設(shè)計實(shí)例7.5.1RC4加密原理與實(shí)現(xiàn)RC4密碼，全稱為“RivestCipher4”，是一種典型的分組密鑰，習(xí)慣上稱之為流密碼，因?yàn)镽C4可用于互聯(lián)網(wǎng)中的實(shí)時數(shù)據(jù)傳輸。RC4的密鑰長度可為1至256個字節(jié)，建議實(shí)際保密通信應(yīng)用中使用128字節(jié)以上的密鑰。這里，設(shè)p表示明文，k表示密鑰，c表示密文，均為基于字節(jié)的向量。結(jié)合表7.6可知，對于RC4加密過程，輸入為密鑰k和長度為n個字節(jié)的明文p，輸出為長度為n個字節(jié)的密文c。具體的加密步驟如下：（1）密碼流初始化第1步：將密鑰k擴(kuò)展為長度為256字節(jié)的key。設(shè)密鑰k的長度為m個字節(jié)，則key[i++]=k[(i++)modm],i=0,1,2,…,255(7.4)第2步：初始化長度為256字節(jié)的數(shù)組sbox，即sbox=[0,1,2,…,255]。第3步：循環(huán)變量i從0至255，循環(huán)執(zhí)行以下兩條語句：(i)j=(j+sbox[i]+key[i])mod256RC4加密過程如圖7-6所示(ii)互換sbox[i]與sbox[j]的值。經(jīng)過上述3步得到的sbox稱為初始密碼流。（2）加密算法已知明文p的長度為n。初始化變量i=0、j=0。變量u從0至n-1，循環(huán)執(zhí)行以