項(xiàng)目八假設(shè)檢驗(yàn)回歸分析與方差分析

1、項(xiàng)目八 假設(shè)檢驗(yàn)、回歸分析與方差分析實(shí)驗(yàn)1 假設(shè)檢驗(yàn)實(shí)驗(yàn)?zāi)康?掌握用Mathematica作單正態(tài)總體均值、方差的假設(shè)檢驗(yàn), 雙正態(tài)總體的均值差、方差比的假設(shè)檢驗(yàn)方法, 了解用Mathematica作分布擬合函數(shù)檢驗(yàn)的方法.單正態(tài)總體均值的假設(shè)檢驗(yàn)(方差已知情形)例1.1 (教材 例1.1) 某車間生產(chǎn)鋼絲, 用表示鋼絲的折斷力, 由經(jīng)驗(yàn)判斷, 其中, 今換了一批材料, 從性能上看, 估計(jì)折斷力的方差不會(huì)有什么變化(即仍有), 但不知折斷力的均值和原先有無差別. 現(xiàn)抽得樣本, 測得其折斷力為578 572 570 568 572 570 570 572 596 584取試檢驗(yàn)折斷力均值有無變化

2、?根據(jù)題意, 要對均值作雙側(cè)假設(shè)檢驗(yàn) 輸入<<StatisticsHypothesisTests.m執(zhí)行后, 再輸入data1=578,572,570,568,572,570,570,572,596,584;MeanTestdata1,570,SignificanceLevel->0.05,KnownVariance->64,TwoSided->True,FullReport->True (*檢驗(yàn)均值, 顯著性水平, 方差已知*)則輸出結(jié)果FullReport->MeanTestStatDistribution575.22.05548 NormalDi

3、stributionTwoSidedPValue->0.0398326,Reject null hypothesis at significance level ->0.05即結(jié)果給出檢驗(yàn)報(bào)告: 樣本均值, 所用的檢驗(yàn)統(tǒng)計(jì)量為統(tǒng)計(jì)量(正態(tài)分布),檢驗(yàn)統(tǒng)計(jì)量的觀測值為2.05548, 雙側(cè)檢驗(yàn)的值為0.0398326, 在顯著性水平下, 拒絕原假設(shè), 即認(rèn)為折斷力的均值發(fā)生了變化.例1.2 (教材 例1.2) 有一工廠生產(chǎn)一種燈管, 已知燈管的壽命X服從正態(tài)分布, 根據(jù)以往的生產(chǎn)經(jīng)驗(yàn), 知道燈管的平均壽命不會(huì)超過1500小時(shí). 為了提高燈管的平均壽命, 工廠采用了新的工藝. 為了弄清

4、楚新工藝是否真的能提高燈管的平均壽命,他們測試了采用新工藝生產(chǎn)的25只燈管的壽命. 其平均值是1575小時(shí), 盡管樣本的平均值大于1500小時(shí), 試問: 可否由此判定這恰是新工藝的效應(yīng), 而非偶然的原因使得抽出的這25只燈管的平均壽命較長呢?根據(jù)題意, 需對均值的作單側(cè)假設(shè)檢驗(yàn) 檢驗(yàn)的統(tǒng)計(jì)量為 , 輸入p1=NormalPValue(1575-1500)/200*Sqrt25ResultOfTestp12,SignificanceLevel ->0.05,FullReport ->True執(zhí)行后的輸出結(jié)果為OneSidedPValue ->0.0303964OneSidedP

5、Value->0.0303964, Fail to reject null hypothesis at significance level ->0.05即輸出結(jié)果拒絕原假設(shè)單正態(tài)總體均值的假設(shè)檢驗(yàn)(方差未知情形)例1.3 (教材 例1.3) 水泥廠用自動(dòng)包裝機(jī)包裝水泥, 每袋額定重量是50kg, 某日開工后隨機(jī)抽查了9袋, 稱得重量如下:49.6 49.3 50.1 50.0 49.2 49.9 49.8 51.0 50.2設(shè)每袋重量服從正態(tài)分布, 問包裝機(jī)工作是否正常()?根據(jù)題意, 要對均值作雙側(cè)假設(shè)檢驗(yàn): 輸入data2=49.6,49.3,50.1,50.0,49.2,4

6、9.9,49.8,51.0,50.2;MeanTestdata2,50.0,SignificanceLevel ->0.05,FullReport ->True(*單邊檢驗(yàn)且未知方差,故選項(xiàng)TwoSided,KnownVariance均采用缺省值*)執(zhí)行后的輸出結(jié)果為FullReport->Mean TestStat Distribution,49.9 -0.559503 StudentTDistribution8OneSidedPValue ->0.295567,Fail to reject null hypothesis at significance level

7、->0.05即結(jié)果給出檢驗(yàn)報(bào)告: 樣本均值, 所用的檢驗(yàn)統(tǒng)計(jì)量為自由度8的分布(檢驗(yàn)),檢驗(yàn)統(tǒng)計(jì)量的觀測值為-0.559503, 雙側(cè)檢驗(yàn)的值為0.295567, 在顯著性水平下, 不拒絕原假設(shè), 即認(rèn)為包裝機(jī)工作正常.例1.4 (教材 例1.4) 從一批零件中任取100件,測其直徑,得平均直徑為5.2,標(biāo)準(zhǔn)差為1.6.在顯著性水平下,判定這批零件的直徑是否符合5的標(biāo)準(zhǔn).根據(jù)題意, 要對均值作假設(shè)檢驗(yàn):檢驗(yàn)的統(tǒng)計(jì)量為, 它服從自由度為的分布. 已知樣本容量 樣本均值, 樣本標(biāo)準(zhǔn)差.輸入StudentTPValue(5.2-5)/1.6*Sqrt100,100-1,TwoSided->

8、;True則輸出TwoSidedPValue->0.214246即值等于0.214246, 大于0.05, 故不拒絕原假設(shè), 認(rèn)為這批零件的直徑符合5的標(biāo)準(zhǔn).單正態(tài)總體的方差的假設(shè)檢驗(yàn)例1.5 (教材 例1.5) 某工廠生產(chǎn)金屬絲, 產(chǎn)品指標(biāo)為折斷力. 折斷力的方差被用作工廠生產(chǎn)精度的表征. 方差越小, 表明精度越高. 以往工廠一直把該方差保持在64(kg)與64以下. 最近從一批產(chǎn)品中抽取10根作折斷力試驗(yàn), 測得的結(jié)果(單位為千克) 如下:578 572 570 568 572 570 572 596 584 570由上述樣本數(shù)據(jù)算得.為此, 廠方懷疑金屬絲折斷力的方差是否變大了.

9、如確實(shí)增大了, 表明生產(chǎn)精度不如以前, 就需對生產(chǎn)流程作一番檢驗(yàn), 以發(fā)現(xiàn)生產(chǎn)環(huán)節(jié)中存在的問題.根據(jù)題意, 要對方差作雙邊假設(shè)檢驗(yàn): 輸入data3=578,572,570,568,572,570,572,596,584,570;VarianceTestdata3,64,SignificanceLevel->0.05,FullReport->True(*方差檢驗(yàn),使用雙邊檢驗(yàn),*)則輸出FullReport->Variance TestStat Distribution75.7333 10.65 ChiSquareDistribution9OneSidedPValue->

10、;0.300464,Fail to reject null hypothesis at significance level->0.05即檢驗(yàn)報(bào)告給出: 樣本方差 所用檢驗(yàn)統(tǒng)計(jì)量為自由度4的分布統(tǒng)計(jì)量(檢驗(yàn)), 檢驗(yàn)統(tǒng)計(jì)量的觀測值為10.65, 雙邊檢驗(yàn)的值為0.300464, 在顯著性水平時(shí), 接受原假設(shè), 即認(rèn)為樣本方差的偏大系偶然因素, 生產(chǎn)流程正常, 故不需再作進(jìn)一步的檢查.例1.6 (教材 例1.6) 某廠生產(chǎn)的某種型號的電池, 其壽命(以小時(shí)計(jì)) 長期以來服從方差的正態(tài)分布, 現(xiàn)有一批這種電池, 從它的生產(chǎn)情況來看, 壽命的波動(dòng)性有所改變. 現(xiàn)隨機(jī)取26只電池, 測出其壽命的

11、樣本方差.問根據(jù)這一數(shù)據(jù)能否推斷這批電池的壽命的波動(dòng)性較以往的有顯著的變化(取)?根據(jù)題意, 要對方差作雙邊假設(shè)檢驗(yàn): 所用的檢驗(yàn)統(tǒng)計(jì)量為它服從自由度為的分布.已知樣本容量 樣本方差輸入ChiSquarePValue(26-1)*9200/5000, 26-1,TwoSided->True則輸出TwoSidedPValue->0.0128357.即值小于0.05, 故拒絕原假設(shè). 認(rèn)為這批電池壽命的波動(dòng)性較以往有顯著的變化.雙正態(tài)總體均值差的檢驗(yàn)(方差未知但相等)例1.7 (教材 例1.7) 某地某年高考后隨機(jī)抽得15名男生、12名女生的物理考試成績?nèi)缦?男生: 49 48 47

12、53 51 43 39 57 56 46 42 44 55 44 40女生: 46 40 47 51 43 36 43 38 48 54 48 34從這27名學(xué)生的成績能說明這個(gè)地區(qū)男女生的物理考試成績不相上下嗎?(顯著性水平).根據(jù)題意, 要對均值差作單邊假設(shè)檢驗(yàn): 輸入data4=49.0,48,47,53,51,43,39,57,56,46,42,44,55,44,40;data5=46,40,47,51,43,36,43,38,48,54,48,34;MeanDifferenceTestdata4,data5,0,SignificanceLevel->0.05,TwoSided-

13、>True,FullReport->True,EqualVariances->True,FullReport->True (*指定顯著性水平,且方差相等*)則輸出FullReport->MeanDiff TestStat Distribution3.6 1.56528 tudentTDistribution25,OneSidedPValue->0.13009,Fail to reject null hypothesis at significance level->0.05即檢驗(yàn)報(bào)告給出: 兩個(gè)正態(tài)總體的均值差為3.6, 檢驗(yàn)統(tǒng)計(jì)量為自由度25的分布(

14、檢驗(yàn)),檢驗(yàn)統(tǒng)計(jì)量的觀察值為1.56528, 單邊檢驗(yàn)的值為0.13009, 從而沒有充分理由否認(rèn)原假設(shè), 即認(rèn)為這一地區(qū)男女生的物理考試成績不相上下.雙正態(tài)總體方差比的假設(shè)檢驗(yàn)例1.8 (教材 例1.8) 為比較甲、乙兩種安眠藥的療效, 將20名患者分成兩組, 每組10人, 如服藥后延長的睡眠時(shí)間分別服從正態(tài)分布, 其數(shù)據(jù)為(單位:小時(shí)):甲: 5.5 4.6 4.4 3.4 1.9 1.6 1.1 0.8 0.1 -0.1 乙: 3.7 3.4 2.0 2.0 0.8 0.7 0 -0.1 -0.2 -1.6問在顯著性水平下兩重要的療效又無顯著差別.根據(jù)題意, 先在未知的條件下檢驗(yàn)假設(shè):

15、輸入list1=5.5,4.6,4.4,3.4,1.9,1.6,1.1,0.8,0.1,-0.1;list2=3.7,3.4,2.0,2.0,0.8,0.7,0,-0.1,-0.2,-1.6;VarianceRatioTestlist1,list2,1,SignificanceLevel->0.05,TwoSided->True,FullReport->True (*方差比檢驗(yàn),使用雙邊檢驗(yàn),*)則輸出FullReport->Ratio TestStat Distribution1.41267 1.41267 FratioDistribution9,9,TwoSided

16、PValue->0.615073,Fail to reject null hypothesis at significancelevel->0.05即檢驗(yàn)報(bào)告給出: 兩個(gè)正態(tài)總體的樣本方差之比為1.41267, 檢驗(yàn)統(tǒng)計(jì)量的分布為分布(F檢驗(yàn)), 檢驗(yàn)統(tǒng)計(jì)量的觀察值為1.41267, 雙側(cè)檢驗(yàn)的值為0.615073. 由檢驗(yàn)報(bào)告知兩總體方差相等的假設(shè)成立. 其次, 要在方差相等的條件下作均值是否相等的假設(shè)檢驗(yàn): 輸入MeanDifferenceTestlist1,list2,0,EqualVariances->True,SignificanceLevel->0.05,T

17、woSided->True,FullReport->True (*均值差是否為零的檢驗(yàn),已知方差相等,雙邊檢驗(yàn)*)則輸出FullReport->MeanDiff TestStat Distribution1.26 1.52273 StudentTDistribution18,TwoSidedPValue->0.1452,Fail to reject null hypothesis at significance level->0.05根據(jù)輸出的檢驗(yàn)報(bào)告, 應(yīng)接受原假設(shè)因此, 在顯著性水平下可認(rèn)為. 綜合上述討論結(jié)果, 可以認(rèn)為兩種安眠藥療效無顯著差異.例1.9 (

18、教材 例1.9) 甲、乙兩廠生產(chǎn)同一種電阻, 現(xiàn)從甲乙兩廠的產(chǎn)品中分別隨機(jī)抽取12個(gè)和10個(gè)樣品, 測得它們的電阻值后, 計(jì)算出樣本方差分別為 假設(shè)電阻值服從正態(tài)分布, 在顯著性水平下, 我們是否可以認(rèn)為兩廠生產(chǎn)的電阻值的方差相等.根據(jù)題意, 檢驗(yàn)統(tǒng)計(jì)量為它服從自由度()的分布.已知樣本容量, 樣本方差該問題即檢驗(yàn)假設(shè): 輸入 FRatioPValue1.40/4.38,12-1,10-1,TwoSided->True,SignificanceLevel->0.1則輸出TwoSidedPValue->0.0785523,Reject null hypothesis at si

19、gnificance level->0.1所以, 我們拒絕原假設(shè), 即認(rèn)為兩廠生產(chǎn)的電阻阻值的方差不同分布擬合檢驗(yàn)檢驗(yàn)法例1.10 (教材 例1.10) 下面列出84個(gè)伊特拉斯坎男子頭顱的最大寬度(單位:mm):141 148 132 138 154 142 150 146 155 158 150 140 147 148 144150 149 145 149 158 143 141 144 144 126 140 144 142 141 140145 135 147 146 141 136 140 146 142 137 148 154 137 139 143140 131 143 14

20、1 149 148 135 148 152 143 144 141 143 147 146150 132 142 142 143 153 149 146 149 138 142 149 142 137 134144 146 147 140 142 140 137 152 145試檢驗(yàn)上述頭顱的最大寬度數(shù)據(jù)是否來自正態(tài)總體()?輸入數(shù)據(jù)data2=141,148,132,138,154,142,150,146,155,158,150,140, 147,148,144,150,149,145,149,158,143,141,144,144,126,140, 144,142,141,140,145,

21、135,147,146,141,136,140,146,142,137, 148,154,137,139,143,140,131,143,141,149,148,135,148,152, 143,144,141,143,147,146,150,132,142,142,143,153,149,146, 149,138,142,149,142,137,134,144,146,147,140,142,140,137,152,145;輸入Mindata2|Maxdata2則輸出126|158即頭顱寬度數(shù)據(jù)的最小值為126, 最大值為158. 考慮區(qū)間124.5,159.5, 它包括了所有的數(shù)據(jù). 以5

22、為間隔, 劃分小區(qū)間. 計(jì)算落入每個(gè)小區(qū)間的頻數(shù), 輸入pshu=BinCountsdata2,124.5,159.5,5則輸出1,4,10,33,24,9,3因?yàn)槌霈F(xiàn)了兩個(gè)區(qū)間內(nèi)的頻數(shù)小于5, 所以要合并小區(qū)間. 現(xiàn)在把頻數(shù)為1, 4的兩個(gè)區(qū)間合并, 再把頻數(shù)為9, 3的兩個(gè)區(qū)間合并. 這樣只有5個(gè)小區(qū)間. 這些區(qū)間為(),為了計(jì)算分布函數(shù)在端點(diǎn)的值, 輸入zu=Table129.5+j*5,j,1,4則輸出134.5,139.5,144.5,149.5以這4個(gè)數(shù)為分點(diǎn),把分成5個(gè)區(qū)間后,落入5個(gè)小區(qū)間的頻數(shù)分別為5, 10, 33, 24, 12.它們除以數(shù)據(jù)的總個(gè)數(shù)就得到頻率. 輸入pl

23、v=5,10,33,24,12/Lengthdata2則輸出下面計(jì)算在成立條件下, 數(shù)據(jù)落入5個(gè)小區(qū)間的概率. 輸入nor=NormalDistributionMeandata2,StandardDeviationMLEdata2; (*Meandata2是總體均值的極大似然估計(jì),StandardDeviationMLEdata2是總體標(biāo)準(zhǔn)差的極大似然估計(jì),NormalDistribution是正態(tài)分布,因此nor是由極大似然估計(jì)得到的正態(tài)分布*)Fhat=CDFnor,zu (*CDF是分布函數(shù)的值*)則輸出0.0590736,0.235726,0.548693,0.832687此即成立條件

24、下分布函數(shù)在分點(diǎn)的值. 再求相鄰兩個(gè)端點(diǎn)的分布函數(shù)值之差, 輸入Fhat2=Join0,Fhat,1;glv=TableFhat2j-Fhat2j-1,j,2,LengthFhat2則輸出0.0590736,0.176652,0.312967,0.283994,0.167313輸入計(jì)算檢驗(yàn)統(tǒng)計(jì)量值的命令chi=ApplyPlus,(plv-glv)2/glv*Lengthdata2則輸出3.59235再輸入求分布的值命令ChiSquarePValuechi,2 (*5-2-1=2為分布的自由度*)則輸出OneSidedPValue->0.165932這個(gè)結(jié)果表明成立條件下, 統(tǒng)計(jì)量取3.

25、59235及比它更大的概率為0.165932, 因此不拒絕, 即頭顱的最大寬度數(shù)據(jù)服從正態(tài)分布.實(shí)驗(yàn)2 回歸分析實(shí)驗(yàn)?zāi)康?學(xué)習(xí)利用Mathematica求解一元線性回歸問題. 學(xué)會(huì)正確使用命令線性回歸Regress, 并從輸出表中讀懂線性回歸模型中各參數(shù)的估計(jì), 回歸方程, 線性假設(shè)的顯著性檢驗(yàn)結(jié)果, 因變量Y在預(yù)察點(diǎn)的預(yù)測區(qū)間等.例2.1 (教材 例2.1) 某建材實(shí)驗(yàn)室做陶?；炷翆?shí)驗(yàn)室中, 考察每立方米混凝土的水泥用量(kg)對混凝土抗壓強(qiáng)度的影響, 測得下列數(shù)據(jù):(1) 畫出散點(diǎn)圖;(2) 求y關(guān)于x的線性回歸方程并作回歸分析;(3) 設(shè)kg, 求y的預(yù)測值及置信水平為0.95的預(yù)測區(qū)

26、間.先輸入數(shù)據(jù):aa = 150,56.9,160,58.3,170,61.6,180,64.6,190,68.1,200,71.3,210,74.1,220,77.4,230,80.2,240,82.6,250,86.4,260,89.7; (1) 作出數(shù)據(jù)表的散點(diǎn)圖. 輸入ListPlotaa,PlotRange->140,270,50,90則輸出圖2.1.圖2.1(2) 作一元回歸分析, 輸入Regressaa,1,x,x,RegressionReport->BestFit,ParameterCITable,SummaryReport則輸出BestFit->10.282

27、9+0.303986x, ParameterCITable->EstimateSECI110.28290.8503758.388111,12.1776,x0.3039860.004090580.294872,0.3131ParameterTable->EsimateSETstatPValue110.28290.85037512.09222.71852,x0.3039860.0040905874.31374.88498Rsquared->0.998193,AdjustedRSquared->0.998012,EstimatedVariance->0.0407025,

28、ANOVATable->DFSumOfSqMeanSqFratio PValueModel 1 1321.431321.435522.52 4.77396Error 10 2.39280.23928Total 11 1323.82現(xiàn)對上述回歸分析報(bào)告說明如下:BestFit(最優(yōu)擬合)-> 10.2829+0.303986x表示一元回歸方程為;ParameterCITable(參數(shù)置信區(qū)間表)中: Estimate這一列表示回歸函數(shù)中參數(shù)a, b的點(diǎn)估計(jì)為=10.2829 (第一行), = 0.303986 (第二行); SE這一列的第一行表示估計(jì)量的標(biāo)準(zhǔn)差為0.850375, 第

29、二行表示估計(jì)量的標(biāo)準(zhǔn)差為0.00409058; CI這一列分別表示的置信水平為0.95的置信區(qū)間是(8.388111,12.1776), 的置信水平為0.95的置信區(qū)間是(0.294872,0.3131).ParameterTable(參數(shù)表)中前兩列的意義同參數(shù)置信區(qū)間表; Tstat與Pvalue這兩列的第一行表示作假設(shè)檢驗(yàn)(t檢驗(yàn)):時(shí), T統(tǒng)計(jì)量的觀察值為12.0922, 檢驗(yàn)統(tǒng)計(jì)量的P值為2.71852, 這個(gè)P值非常小, 檢驗(yàn)結(jié)果強(qiáng)烈地否定, 接受; 第二行表示作假設(shè)檢驗(yàn)(t檢驗(yàn)): 時(shí)T統(tǒng)計(jì)量的觀察值為74.3137, 檢驗(yàn)統(tǒng)計(jì)量的P值為4.88498, 這個(gè)P值也非常小, 檢驗(yàn)

30、結(jié)果強(qiáng)烈地否定接受.Rsquared->0.998193, 表示 它說明y的變化有99.8%來自x的變化; AdjustedRSquared->0.998012, 表示修正后的0.998012. EstimatedVariance->0.0407025, 表示線性模型中方差的估計(jì)為0.0407025.ANOVATable(回歸方差分析表)中的DF這一列為自由度: Model(一元線性回歸模型)的自由度為1, Error(殘差)的自由度為Total(總的)自由度為SumOfSq這一列為平方和: 回歸平方和1321.43, 殘差平方和2.3928,總的平方和1323.82;Mea

31、nSq這一列是平方和的平均值, 由SumOfSq這一列除以對應(yīng)的DF得到, 即FRatio這一列為統(tǒng)計(jì)量的值, 即最后一列表示統(tǒng)計(jì)量F的P值非常接近于0. 因此在作模型參數(shù)的假設(shè)檢驗(yàn)(F檢驗(yàn)):時(shí), 強(qiáng)烈地否定, 即模型的參數(shù)向量因此回歸效果非常顯著.(3) 在命令RegressionReport的選項(xiàng)中增加RegressionReport->SinglePredictionCITable就可以得到在變量x的觀察點(diǎn)處的y的預(yù)測值和預(yù)測區(qū)間. 雖然不是觀察點(diǎn), 但是可以用線性插值的方法得到近似的置信區(qū)間. 輸入aa=Sortaa; (*對數(shù)據(jù)aa按照水泥用量x的大小進(jìn)行排序*)regres

32、s2=Regressaa,1,x,x,RegressionReport->SinglePredictionCITable(*對數(shù)據(jù)aa作線性回歸, 回歸報(bào)告輸出y值的預(yù)測區(qū)間*)執(zhí)行后輸出SinglePredictionCITable->Observed PredictedSECI56.955.88080.5566354.6405,57.12158.358.92060.54139157.7143,60.1269 61.661.96050.52888360.7821,63.1389 64.665.00030.51930563.8433,66.157468.168.04020.5128

33、266.8976,69.182871.371.08010.50954769.9447,72.215474.174.11990.50954772.9846,75.255377.477.15980.5128276.0172,78.302480.280.19970.51930579.0426,81.356782.683.23950.52888382.0611,84.417986.486.27940.54139185.0731,87.485789.789.31920.5566388.079,90.5595上表中第一列是觀察到的y的值, 第二列是y的預(yù)測值, 第三列是標(biāo)準(zhǔn)差, 第四列是相應(yīng)的預(yù)測區(qū)間(置

34、信度為0.95). 從上表可見在時(shí), y的預(yù)測值為77.1598, 置信度為0.95的預(yù)測區(qū)間為(76.0172,75.2553), 在時(shí), y的預(yù)測值為80.1997, 置信度為0.95的預(yù)測區(qū)間為79.0426,81.3567. 利用線性回歸方程, 可算得時(shí), y的預(yù)測值為78.68, 置信度為0.95的預(yù)測區(qū)間為(77.546, 79.814).利用上述插值思想, 可以進(jìn)一步作出預(yù)測區(qū)間的圖形.先輸入調(diào)用圖軟件包命令<<Graphics執(zhí)行后再輸入observed2,predicted2,se2,ci2=Transpose(SinglePredictionCITable/.r

35、egress2)1;(*取出上面輸出表中的四組數(shù)據(jù), 分別記作observed2,predicted2,se2,ci2*)xva12=MapFirst,aa;(*取出數(shù)據(jù)aa中的第一列, 即數(shù)據(jù)中x的值, 記作xva12*)Predicted3=Transposexva12,predicted2;(*把x的值xva12與相應(yīng)的預(yù)測值predicted2配成數(shù)對, 它們應(yīng)該在一條回歸直線上*)lowerCI2=Transposexva12,MapFirst,ci2;(*MapFirst,ci2取出預(yù)測區(qū)間的第一個(gè)值, 即置信下限. x的值xva12與相應(yīng)的置信下限配成數(shù)對*)upperCI2=T

36、ransposexva12,MapLast,ci2;(*MapLast,ci2取出預(yù)測區(qū)間的第二個(gè)值, 即置信上限. x的值xva12與相應(yīng)的置信上限配成數(shù)對*)MultipleListPlotaa,Predicted3,lowerCI2,upperCI2,PlotJoined->False,True,True,True,SymbolShape->PlotSymbolDiamond,None,None, None,PlotStyle->Automatic,Automatic,Dashing0.04,0.04,Dashing0.04,0.04(*把原始數(shù)據(jù)aa和上面命令得到的三

37、組數(shù)對predicted3,lowerCI2,upperCI2用多重散點(diǎn)圖命令MultipleListPlot在同一個(gè)坐標(biāo)中畫出來. 圖形中數(shù)據(jù)aa的散點(diǎn)圖不用線段連接起來, 其余的三組散點(diǎn)圖用線段連接起來, 而且最后兩組數(shù)據(jù)的散點(diǎn)圖用虛線連接.*)則輸出圖2.2.圖2.2從圖形中可以看到, 由Y的預(yù)測值連接起來的實(shí)線就是回歸直線. 鉆石形的點(diǎn)是原始數(shù)據(jù). 虛線構(gòu)成預(yù)測區(qū)間.多元線性回歸例2.2 (教材 例2.2) 一種合金在某種添加劑的不同濃度下, 各做三次試驗(yàn), 得到數(shù)據(jù)如下表:(1) 作散點(diǎn)圖;(2) 以模型擬合數(shù)據(jù), 其中與x無關(guān); (3) 求回歸方程并作回歸分析.先輸入數(shù)據(jù)bb=10

38、.0,25.2,10.0,27.3,10.0,28.7,15.0,29.8,15.0,31.1,15.0,27.8,20.0,31.2,20.0,32.6,20.0,29.7,25.0,31.7,25.0,30.1,25.0,32.3,30.0,29.4,30.0,30.8,30.0,32.8; (1) 作散點(diǎn)圖, 輸入ListPlotbb,PlotRange->5,32,23,33,AxesOrigin->8,24則輸出圖2.3.圖2.3(2) 作二元線性回歸, 輸入Regressbb,1,x,x2,x,RegressionReport->BestFit,Parameter

39、CITable,SummaryReport(*對數(shù)據(jù)bb作回歸分析, 回歸函數(shù)為用1,x,x2表示, 自變量為x, 參數(shù),的置信水平為0.95的置信區(qū)間)執(zhí)行后得到輸出的結(jié)果:2,ParameterCITable->EstimateSECI1 19.03333.2775511.8922,26.1745x 1.008570.3564310.231975,1.78517x2 -0.0203810.00881488-0.0395869,-0.00117497ParameterTable->EstimateSETstatPValue1 19.03333.277555.807180.0000

40、837856x 1.008570.3564312.829640.0151859x2 -0.0203810.00881488-2.312110.0393258Rsquared->0.614021,AdjustedRSquared->0.549692,EstimatedVariance->2.03968,ANOVATable->DF SumOfSqMeanSqFratioPValueMode1 2 38.937119.46869.54490.00330658Error 12 24.47622.03968Total 14 63.4133從輸出結(jié)果可見: 回歸方程為它們的置信水

41、平為0.95的置信區(qū)間分別是(11.8922,26.1745),(0.231975,1.78517),(-0.0395869,-0.00117497).假設(shè)檢驗(yàn)的結(jié)果是: 在顯著性水平為0.95時(shí)它們都不等于零. 模型中,的估計(jì)為2.03968. 對模型參數(shù)是否等于零的檢驗(yàn)結(jié)果是: 因此回歸效果顯著.非線性回歸例2.3 下面的數(shù)據(jù)來自對某種遺傳特征的研究結(jié)果, 一共有2723對數(shù)據(jù), 把它們分成8類后歸納為下表.研究者通過散點(diǎn)圖認(rèn)為y和x符合指數(shù)關(guān)系: 其中是參數(shù). 求參數(shù)的最小二乘估計(jì).因?yàn)閥和x的關(guān)系不是能用Fit命令擬合的線性關(guān)系, 也不能轉(zhuǎn)換為線性回歸模型. 因此考慮用(1)多元微積分

42、的方法求的最小二乘估計(jì); (2)非線性擬合命令NonlinearFit求的最小二乘估計(jì).(1) 微積分方法輸入OffGenera1:spe11OffGenera1:spe111Clearx,y,a,b,cdataset=579,1,38.08,1021,2,29.70,607,3,25.42,324,4,23.15,120,5,21.79,46,6,20.91,17,7,19.37,9,8,19.36; (*輸入數(shù)據(jù)集*)yx_:=a Expb x+c (*定義函數(shù)關(guān)系*)下面一組命令先定義了曲線與2723個(gè)數(shù)據(jù)點(diǎn)的垂直方向的距離平方和, 記為再求對的偏導(dǎo)數(shù)分別記為用FindRoot命令解三個(gè)

43、偏導(dǎo)數(shù)等于零組成的方程組(求解). 其結(jié)果就是所要求的的最小二乘估計(jì). 輸入Cleara,b,c,f,fa,fb,fcga_,b_,c_:=Sumdataseti,1*(dataseti,3-a*Expdataseti,2*b-c)2,i,1,Lengthdatasetgaa_,b_,c_=Dga,b,c,a;gba_,b_,c_=Dga,b,c,b;gca_,b_,c_=Dga,b,c,c;Cleara,b,coursolution=FindRootgaa,b,c=0,gba,b,c=0,gca,b,c=0,a,40.,b,-1.,c,20.(* 40是a的初值, -1是b的初值, 20是c

44、的初值*)則輸出a->33.2221,b->-0.626855,c->20.2913再輸入yhatx_=yx/.oursolution則輸出20.2913+33.2221這就是y和x的最佳擬合關(guān)系. 輸入以下命令可以得到擬合函數(shù)和數(shù)據(jù)點(diǎn)的圖形:p1=Plotyhatx,x,0,12,PlotRange->15,55,DisplayFunction->Identity;pts=Tabledataseti,2,dataseti,3,i,1,Lengthdataset;p2=ListPlotpts,PlotStyle->PointSize.01,DisplayFu

45、nction->Identity;Showp1,p2,DisplayFunction->$DisplayFunction;則輸出圖2.4.圖2.4(2) 直接用非線性擬合命令NonlinearFit方法輸入data2=FlattenTableTabledatasetj,2,datasetj, 3,i,datasetj,1,j,1,Lengthdataset,1;(*把數(shù)據(jù)集恢復(fù)成2723個(gè)數(shù)對的形式*)<<Statisticsw=NonlinearFitdata2,a*Expb*x+c,x,a,40,b,-1,c,20則輸出這個(gè)結(jié)果與(1)的結(jié)果完全相同. 這里同樣要注

46、意的是參數(shù)必須選擇合適的初值.如果要評價(jià)回歸效果, 則只要求出2723個(gè)數(shù)據(jù)的殘差平方和 輸入yest=Tableyhatdataseti,2,i,1,Lengthdataset;yact=Tabledataseti,3,i,1,Lengthdataset;wts=Tabledataseti,1,i,1,Lengthdataset;sse=wts.(yact-yest)2 (*作點(diǎn)乘運(yùn)算*)則輸出59.9664即2723個(gè)數(shù)據(jù)的殘差平方和是59.9664. 再求出2723個(gè)數(shù)據(jù)的總的相對誤差的平方和 輸入sse2=wts.(yact-yest)2/yest) (*作點(diǎn)乘運(yùn)算)則輸出2.7407

47、5由此可見, 回歸效果是顯著的.實(shí)驗(yàn)3 方差分析實(shí)驗(yàn)?zāi)康?學(xué)習(xí)利用Mathematica求單因素方差分析的方法.例3.1 (教材 例3.1) 今有某種型號的電池三批, 它們分別是A,B,C三個(gè)工廠所生產(chǎn)的. 為評比起質(zhì)量, 各隨機(jī)抽取5只電池為樣品, 經(jīng)試驗(yàn)得其壽命(單位:h)如下表:A4042484538B2628343230C3950405043試在顯著性水平0.05下檢驗(yàn)電池的平均壽命有無顯著的差異. 若差異是顯著的, 試求均值差及的置信水平為95%的置信區(qū)間.這是方差分析問題, 先把它轉(zhuǎn)化為線性模型:Y = Xb +e. 令 則線性模型(3.3)與方差分析模型(3.1)完全等價(jià). 模型

48、(3.3)完全可以用DesignedRegress命令作設(shè)計(jì)回歸, 得到所要的方差分析表.我們面臨的任務(wù)是: (1) 檢驗(yàn)3個(gè)總體的均值是否相等,即作假設(shè)檢驗(yàn)(2) 求均值差及的置信水平為95%的置信區(qū)間.任務(wù)(1)等價(jià)于對模型(3.3)作檢驗(yàn):而任務(wù)(2)等價(jià)于求的置信區(qū)間. 在DesignedRegress命令中加入選項(xiàng)RegressionReport->ParameterCITable,MeanPredictionCITable,SummaryReport后便能完成上述任務(wù). 用回歸分析作單因素方差分析 完成對模型的假設(shè)檢驗(yàn)和對模型參數(shù)的區(qū)間估計(jì)任務(wù).輸入設(shè)計(jì)矩陣和數(shù)據(jù)X1=1.0

49、,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,1,0,1,1,0,1,1,0,1,1,0,1,1,0,1,0,1,1,0,1,1,0,1,1,0,1,1,0,1; Y1=40,42,48,45,38,26,28,34,32,30,39,50,40,50,43;再輸入設(shè)計(jì)回歸命令DesignedRegressX1,Y1,RegressionReport-> ParameterCITable,MeanPredictionCITable,SummaryReport(*回歸報(bào)告輸出參數(shù)的置信區(qū)間,均值的置信區(qū)間和總結(jié)報(bào)告*)執(zhí)行后得到輸出Estimate SE CI142.61.

50、8991238.4622,46.7378ParameterCITable->2-12.62.68576-18.4518,-6.74822 31.82.68576-4.05178,7.65178MeanPredictionCITable->ObservedPredictedSECI40.42.61.8991238.4622,46.737842.42.61.8991238.4622,46.737848.42.61.8991238.4622,46.737845.42.61.8991238.4622,46.737838.42.61.8991238.4622,46.737826.30.1.8

51、991225.8622,34.137828.30.1.8991225.8622,34.137834.30.1.8991225.8622,34.137832.30.1.8991225.8622,34.137830.30.1.8991225.8622,34.137839.44.41.8991240.2622,48.537850.44.41.8991240.2622,48.537840.44.41.8991240.2622,48.537850.44.41.8991240.2622,48.537843.44.41.8991240.2622,48.5378Estimate SE TStatPValue142.61.8991222.43143.63987×10-11ParameterCITable->2-12.62.68576-4.69140.00052196 31.82.685760.67020.515421Rsquared->0.739904,AdjustedRSquared->0.696554,EstimatedVariance->18.0333,ANOVATable-> DFSumOfsqMeanSqFratioPvalueModel2615.6307.8