【建模教程】-建模-數(shù)學建模夏令營

----宋停云與您分享--------宋停云與您分享----2012年西安交通大學數(shù)學建模夏令營參 賽 承 諾 書我們仔細閱讀了中國大學生數(shù)學建模競賽的競賽規(guī)則.我們完全明白，在競賽開始后參賽隊員不能以任何方式（包括電話、電子郵件、網(wǎng)上咨詢等）與隊外的任何人（包括指導教師）研究、討論與賽題有關的問題。我們知道，抄襲別人的成果是違反競賽規(guī)則的,如果引用別人的成果或其他公開的資料（包括網(wǎng)上查到的資料，必須按照規(guī)定的的表述方式在正文引用處和中明確列出。我們鄭重承諾，嚴格遵守競賽規(guī)則，以保證競賽的公正、公平性。如有違反競賽規(guī)則的行為，我們將受到嚴肅處理。我們參賽選擇的題號是（A/B/C/D中選擇一項填寫：我們的參賽報名號為（如果賽區(qū)設置報名號的話：所屬學校（請?zhí)顚懲暾娜?打印并簽名)：1.2.3.指導教師或指導教師組負責人 (打印并簽名)：日期： 年 月 日賽區(qū)評閱編號（由賽區(qū)組委會評閱前進行編號）：----宋停云與您分享--------宋停云與您分享----2012西安交通大學數(shù)學建模夏令營編號專用頁賽區(qū)評閱編號（由賽區(qū)組委會評閱前進行編號）：賽區(qū)評閱記錄（可供賽區(qū)評閱時使用）：評評閱人評分備注全國統(tǒng)一編號（由賽區(qū)組委會送交全國前編號）：全國評閱編號（由全國組委會評閱前進行編號）：----宋停云與您分享--------宋停云與您分享----基于微分方程模型的長江江豚種群數(shù)量變化預測摘要作為目前長江里唯一的淡水哺乳動物，江豚屬于我國二級保護動物 ,目前由于長江流域水體污染加劇、人類肆意采挖江砂、非法使用漁具等原因，長江江豚的生存和繁育受到嚴重影響，其數(shù)量正在銳減。故準確預測長江江豚的數(shù)量變化對我們保護江豚具有深遠意義。本文從長江江豚近二十年的實際生存情況和長江的生態(tài)環(huán)境出發(fā)，建立了基a+be針對題目所給的問題，根據(jù)問題一和問題二的要求，我們先對功能性滅絕的ct=1t13微分方程功能性滅絕徹底滅絕 生殖率 死亡率 保護措施----宋停云與您分享--------宋停云與您分享----一．問題重述根據(jù)題給信息，江豚是目前長江里唯一的淡水哺乳動物，分布于長江根據(jù)題給信息，江豚是目前長江里唯一的淡水哺乳動物，分布于長江中下游干流以及洞庭湖和鄱陽湖等區(qū)域，近 20年來種群數(shù)量快速衰減。根據(jù)題述的資料顯示，1991年長江江豚數(shù)量是 2700多頭；2006年，國際聯(lián)合科學考察組經(jīng)一個多月調查發(fā)現(xiàn)，江豚數(shù)量已不足18002011們的數(shù)量可能僅為 1000余頭，洞庭湖僅剩 85頭江豚。資料顯示，長江江豚的生存面臨巨大威脅。截止 2012年，由于長江流域水體污染加劇、人類肆意采挖江砂、非法使用漁具等原因，長江江豚的生存和繁育受到嚴重影響。長江江豚的數(shù)量只剩下12006.4%10～15本題給出的數(shù)據(jù)和相關資料都比較少，故要求參賽者搜集大量相關資料，建立合理的數(shù)學模型研究長江江豚種群數(shù)量。本文主要研究以下幾個問題：問題一：長江江豚出現(xiàn)功能性滅絕的時間和數(shù)量；問題二：長江江豚出現(xiàn)功能性滅絕至徹底滅絕的時間；問題三：今年3月份以來，洞庭湖連續(xù)發(fā)現(xiàn)江豚死亡,這說明什么？問題四：根據(jù)你的分析，給出保護這一珍稀物種的建議。二．符號說明x(t) t1991（單位：頭；長江江豚的年生殖率（20%[1]a 外界環(huán)境因素（長江流域水體污染加劇、人類肆意采挖江砂、非法使用漁具等）導致的長江江豚年死亡率，即直接死亡率，這里初步設為正常量；b、c用于確定由長江江豚本身數(shù)量的減少導致的死亡率的參系數(shù)（存在長江流域水體污染加劇、人類肆意采挖江砂、非法使用漁具等外界干擾）；t1 長江江豚開始出現(xiàn)功能性滅絕的時間，即第 t1年開始出現(xiàn)功能性滅絕；m1 長江江豚開始出現(xiàn)功能性滅絕時其數(shù)量（單位：頭）；----宋停云與您分享--------宋停云與您分享--------宋停云與您分享--------宋停云與您分享----t 長江江豚徹底滅絕的時間，即第 t2 2

年徹底滅絕；----宋停云與您分享--------宋停云與您分享----△t長江江豚從開始功能性滅絕至徹底滅絕的時間（單位 :年）。三．模型假設1.假設雌性長江江豚都正常分娩，江豚的年生殖率變化不大，即為目前估20%2．理想情況下假設長江江豚雌雄性別比其例為 1:1，且雌性江豚每年10月份生產，每年一胎，每胎只產一仔； [2]3.假設長江江豚的死亡率不受年齡和性別影響，只考慮種群整體死亡率；4．根據(jù)目前長江江豚數(shù)量急劇下降的現(xiàn)象，不考慮江豚因為數(shù)量達到環(huán)境容量后，環(huán)境對其數(shù)量的增加產生的阻力；5．假設政府不出臺保護政策，環(huán)境也沒有進一步惡化。四：問題分析問題一的分析對于問題一，由于目前研究長江江豚種群數(shù)量變化的相關數(shù)據(jù)統(tǒng)計偏少，我們直接構造最原始的針對一般生物種群的數(shù)量變化的常微分方程模型，即：種群數(shù)量增長率（以年為微分元） =種群年生殖率-種群年死亡率。并且根據(jù)題中提供的數(shù)據(jù)，我們假定 1991年為第一年，即t=1,其余年份依次類推。顯然根據(jù)事實長江江豚種群的增長率是負的，即其目前的數(shù)量變化呈銳減形式。首先，這里種群生殖率根據(jù)假設為 0.2，死亡率由內外因共同決定。而我們在構造內外、因導致的江豚的死亡率函數(shù)時考慮了以下兩個事實：（1）當江豚的數(shù)量比較多時，由環(huán)境造成的江豚的死亡數(shù)也會多一些，當江豚的數(shù)量比較少時，由環(huán)境造成的江豚的死亡數(shù)也會少一些；（2）當江豚的數(shù)量越來越少，近親交配的幾率越來越大，由內因造成的死亡率會不斷上升。----宋停云與您分享--------宋停云與您分享----結合模型假設5和事實（1）我們可以認為由環(huán)境造成的江豚的死亡率每年都是相同的，即外因引起的死亡率部分可根據(jù)假設設為 a，a是待求解的常參數(shù)。而對于內因引起的死亡率部分，由（ 2）這條規(guī)律，再聯(lián)系到生物學當中的s型和j型曲線，我們用指數(shù)函數(shù)來模擬長江江豚種群內因引起的死亡率，設為 bect（這里b,c均為帶求解的正常量）。這樣就初步確立了長江江豚數(shù)量變化關于時間的微分方程。其次，問題一要求我們求解長江江豚出現(xiàn)功能性滅絕的時間和數(shù)量。經(jīng)過查找相關資料我們知道，生物種群出現(xiàn)功能性滅絕是指該物種因其生存環(huán)境被破壞，數(shù)量非常稀少，以致在自然狀態(tài)下基本喪失了維持繁殖的能力，甚至喪失了維持生存的能力。功能性滅絕是對生物生存狀況的一個客觀準確描述，主要包含兩個層面的意思，一是指一個物種由于其種群數(shù)量非常稀少，在其生存的生態(tài)系統(tǒng)中不再扮演重要角色，即指其生態(tài)功能的喪失。第二個層面的意思是指殘存種群中已經(jīng)沒有能夠繁殖的個體，或者由于種群數(shù)量稀少，受人類活動干擾及遺傳因素的影響，該物種在自然環(huán).[3]x(t)曲線的走勢來看，ct=1ct<1時江豚的死亡率整體上遠遠小于ct>1時江豚的死亡率，即ct>1以后江豚死亡率突然急劇增長。故我們定義江豚開始出現(xiàn)功能性滅絕的時間點即為1ct=1t1c年長江江豚開始出現(xiàn)功能性滅亡。結合題給數(shù)據(jù)，我們用matlab軟件求解微分方程的解，從而確定長江江豚開始出現(xiàn)功能性滅絕的時間和與此時間點相應的江豚數(shù)量。問題二的分析問題二要求解長江江豚出現(xiàn)功能性滅絕至徹底滅絕的時間。 根據(jù)本題第一問的分析可以求出長江江豚出現(xiàn)功能性滅絕的時間為 t1。我們只需要解出問題一中微分方程解的零點對應的時間 t2（或者接近零點的點對應的時間）。此時，長江江豚從開始出現(xiàn)功能性滅絕至全部滅絕的時間△ t得到求解，即△t=t2-t1。----宋停云與您分享--------宋停云與您分享----問題三的分析問題三中指明：今年 3月份以來，洞庭湖連續(xù)發(fā)現(xiàn)江豚死亡，我們要結合我們的研究模型來解釋這一現(xiàn)象說明什么在我們的模型中江豚數(shù)量增長率由內因和外因兩方面決定，故我們直接從這兩方面著手分析這一現(xiàn)象說明了什么，從而揭示這一現(xiàn)象發(fā)生的主要原因。問題四的分析五．模型的建立與求解常微分方程模型預測長江江豚出現(xiàn)功能性滅絕的時間和數(shù)量下：=0.2（20%；a（長導致的長江江豚年b、c用于確定由長江江豚本身數(shù)量的減少導致的死亡率的參系數(shù)（存在長江流域水體污染加劇、人類肆意采挖江砂、非法使用漁具等外界干擾），則得帶參數(shù)的微分方程為：dx(t)=x(t)[-（a+bect）]dtx(1)=2700(1991年江豚數(shù)量)x(16)=1800（2006）x(22)=1200（2012）此外題中還給出 2006年的江豚數(shù)據(jù)可能不足 1000余頭，微分方程模型中的三個初值參數(shù)數(shù)據(jù)也是估計值，這樣解出的方程可能有較大誤差，因此我們選用下面一種方法誤差相對較小的方法來求解三個參數(shù) abc，然后回代入上面模型中的微分方程，從而求解 x(t).----宋停云與您分享--------宋停云與您分享----用原始的定義法和 MATLAB軟件求解三個參數(shù) a、b、c首先，按照符號約定，x(t+1)表示的是從1991年開始，過了 t年后江豚數(shù)量,則我們可以得到如下所示的一種連乘積形式的遞推關系式：x(t+1)=2700(a-bec)(a-be2c)...(a-betc)，再由題給定的已知條件，我們知道x(16)=1800...(1)x(22)=1200...(2)聯(lián)立方程(1)(2)MATLAB軟件求解其次，因為該方程組中未知數(shù)大于方程數(shù)，所以 MATLAB給出了多組解，但是在這多組解當中有一個最優(yōu)解，可以求得當a=0.1202,b=0.0828,c=0.0303 時，誤差最小，方程（1）和方程（2）的最小誤差數(shù)量級約為 10。再利用 a,b,c 的值我們可以算出x(11)=2242...(3) ，再聯(lián)立方程（1）、（2）、（3）重新解得a=0.1237,b=0.0794,c=0.0312 ，并且此時方程（1）、方程（2）和方程（3）的誤差數(shù)量級縮小至約為 1011。至此，在較小的誤差情況下，我們求出了原模型中三個待求解的常數(shù)參數(shù)，分別為 a=0.1237,b=0.0794,c=0.0312 。參數(shù)回代求解微分方程模型由上一步求解得三個參數(shù)值 a=0.1237,b=0.0794,c=0.0312 此時微分方程模型成為不帶參數(shù)的微分方程 ,且只需一個初值即可求解 x(t),這里取x(1)=2700,則得到:dx(t)=x(t)[0.2-（0.1237+0.0794e0.0312t）]dtx(1)=2700(1991年江豚數(shù)量)----宋停云與您分享--------宋停云與您分享----求解該微分方程得ln(x(t))= 0.0763t-0.07940.0794

e0.0312t+10.4502,進一步----宋停云與您分享--------宋停云與您分享----求解該方程可得到x(t)=e(0.0763t e0.0312t)+e10.4502可以看出長江江0.0312其結果的圖示如下：----宋停云與您分享--------宋停云與您分享--------宋停云與您分享--------宋停云與您分享----300025002000頭數(shù) 1500豚江( 1000500

江豚數(shù)量隨時間變化的模擬曲線x(t)----宋停云與您分享--------宋停云與您分享--------宋停云與您分享--------宋停云與您分享----（（年0 10 20 30 40t-時間/年

50 60 70 80----宋停云與您分享--------宋停云與您分享----圖1 微分方程預測江豚數(shù)量變化曲線----宋停云與您分享--------宋停云與您分享----最后，問題一可求解得到t

1 1= = =32.0513t1 c 0.0312

=32，1----宋停云與您分享--------宋停云與您分享----即在以1991為第一年的假設前提下第32年長江江豚開始出現(xiàn)功能性滅絕，亦即在2022年長江江豚開始出現(xiàn)功能性滅絕。此時江豚的數(shù)量m=x(32)=397頭。換句話說，若以現(xiàn)在的 2012年為第一年的話，再過10年1江豚開始出現(xiàn)功能性滅絕，和題述預測 10-15年基本一致。該模型相對較準----宋停云與您分享--------宋停云與您分享----確的預測了江豚開始出現(xiàn)功能性滅絕的時間。即 t解答。圖示為：

=32,m1

=397為問題一的1----宋停云與您分享--------宋停云與您分享--------宋停云與您分享--------宋停云與您分享----300025002000頭數(shù) 1500豚江( 1000

x(t)----宋停云與您分享--------宋停云與您分享--------宋停云與您分享--------宋停云與您分享----500

X:32Y:397.5(32,397)功能滅絕點----宋停云與您分享--------宋停云與您分享--------宋停云與您分享--------宋停云與您分享----0(1,0)1991年0 10 20 30 40t-時間/年

50 60 70 80----宋停云與您分享--------宋停云與您分享----圖2 求解江豚開始功能滅絕的時間和此時江豚數(shù)量．1借助問題一中的微分方程求解江豚徹底滅絕時間0.0794已由問題一求得ln(x(t))= 0.0763t-0.0312 e0.0312t+10.4502。各年份長江江豚數(shù)量的預測值（各個數(shù)據(jù)均是由軟件求解值經(jīng)取整處理后所得）由x(t)函數(shù)表達式我們可以預測每一年長江江豚的種群數(shù)量如下面表格所示：表1 年-2050年長江江豚種群數(shù)量預測值各年份長江江豚數(shù)量的預測值（各個數(shù)據(jù)均是由軟件求解值經(jīng)取整處理后所得）年份（年）1991199219931994199519961997199819992000預測值（頭）2700268126562623258425382486242723612290年份（年）2001200220032004200520062007200820092010預測值（頭）2214213320471958186517701672157414751375年份（年）2011201220132014201520162017201820192020預測值（頭）127711801085993903818736659586518----宋停云與您分享--------宋停云與您分享----年份（年）2021202220232024202520262027202820292030預測值（頭）455397345297254215181152126103年份（年）2031203220332034203520362037203820392040預測值（頭）8468544334262015129年份（年）2041204220432044204520462047204820492050預測值（頭）2021202220232024202520262027202820292030----宋停云與您分享--------宋停云與您分享----其次,求出lnx(t)= 0.0763t-0.07940.0312

e0.0312t+10.4502=0的根即為t2。有----宋停云與您分享--------宋停云與您分享--------宋停云與您分享--------宋停云與您分享----MATLAB軟件可求得t2

=56.3125，四舍五入取整后t2

=56，即在以1991為第一----宋停云與您分享--------宋停云與您分享----年的假設前提下，再經(jīng)過55年長江江豚只剩下1頭，結合徹底滅絕原始定義和各年的江豚數(shù)量預測，估計 2047年江豚會徹底滅絕。換句話說，若以現(xiàn)在的2012年為第一年的話，再過34年江豚只剩下1頭，再過25年江豚徹底滅絕。問題二的解即為：△t=t -t=56-32=24，故根據(jù)我們的模型預測所2 1得數(shù)據(jù)我們說長江江豚出現(xiàn)功能性滅絕至徹底滅絕的時間為 24年。5.2.2江豚從功能性滅絕至徹底滅絕的時間另解江豚的徹底滅絕時間也可通過 5.1.2中的連乘積形式的遞推關系式來求解。我們直接利用數(shù)學軟件求解下面的不等式：x(t+1)=2700(a-bec)(a-be2c)...(a-betc)<1,t>54,t=55,可知x(56)<1, 即在以1991為第一年的假設前提下，再經(jīng)過 55年長江江豚徹底滅絕，這和5.2.1的結果一致。具體時間圖示如下：----宋停云與您分享--------宋停云與您分享----藍色曲線代表的是根據(jù)我們的微分方程模型求解出來的 x(t)x(t)300025002000頭數(shù) 1500豚江( 1000----宋停云與您分享--------宋停云與您分享----5000(1,0)1991年0 10 20 30 40 50t-時間/年

X:56Y:1.127(56,1)徹底滅絕點60 70 80----宋停云與您分享--------宋停云與您分享----圖3求解長江江豚徹底滅絕的時間至此，我們已經(jīng)全部解決了問題一和問題二，其整體結果如下圖所示：圖示中，藍色曲線代表的是根據(jù)我們的微分方程模型求解出來的 x(t),即長江江豚種群數(shù)量隨時間變化的函數(shù)曲線圖；而紅星點代表長江江豚的數(shù)量的實際值。結合整個圖像我們可以看出，我們的微分方程模型在誤差較小的范圍內比較準確的預測出了各個年份長江江豚種群數(shù)量的變化。----宋停云與您分享--------宋停云與您分享----FirstSecond為實際江豚數(shù)量FirstSecond為實際江豚數(shù)量25002000頭數(shù) 1500豚江( 1000500 ----宋停云與您分享--------宋停云與您分享----0年

)

底滅絕點----宋停云與您分享--------宋停云與您分享----0 20 40 60 80t-時間/年圖示4問題一和問題二求解結果的綜合結果根據(jù)模型解釋洞庭湖連續(xù)發(fā)現(xiàn)江豚死亡現(xiàn)象的啟示----宋停云與您分享--------宋停云與您分享----310號開始的，經(jīng)過禁魚季之前長時再次，人類肆意采挖江砂導致江豚生活環(huán)境的變化也可能是此次死亡的原因最后，造紙廠等工廠排放的廢水導致水體污染使江豚中毒也可能是原因之一。保護江豚這一珍稀物種的建議從上面的分析可以看出，環(huán)境的惡化直接導致江豚死亡，江豚的死亡又產生具體保護措施可有以下幾方面：就地保護或者遷地保護；加強對湖污水進行嚴格監(jiān)測監(jiān)管；依法嚴厲打擊非法捕撈行為；繼續(xù)依法嚴厲整治河道非法采砂行為；成立江豚救治小組；建立和完善江豚保護的管理機制體制；各個保護區(qū)抓緊啟動江豚館建設的前期工作；加強對群眾保護江豚的宣傳教育工作等；（10）----宋停云與您分享--------宋停云與您分享----六：模型優(yōu)缺點評價對于問題一、二、三、四我們使用了同一個微分方程，微分方程的優(yōu)點在于最后綜合考慮模型的優(yōu)缺點如下：本模型用常量來表示環(huán)境導致的死亡率，七．模型的改進對于問題一、二、三、四的微分方程模型，由于我們考慮的影響因素有限導----宋停云與您分享--------宋停云與您分享----八．張先鋒，武漢。江豚的年齡鑒定、生長和生殖的研究。水生生物學報。第十六卷第四期，199212月/view/68146.htm/view/1290297.htm/view/445269.htm----宋停云與您分享--------宋停云與您分享--------宋停云與您分享--------宋停云與您分享----統(tǒng)計建模與R軟件第七章習題答案(方差分析)Ex7.1(1)>lamp<-data.frame(X=c(115,116,98,83,103,107,118,116,73,89,85,97),A=factor(rep(1:3,c(4,4,4))))>lamp.aov<-aov(X~A,data=lamp);summary(lamp.aov)DfSumSqMeanSqFvaluePr(>F)A 2 1304 652.0 4.9230.0359*Residuals 9 1192 P0.05，有顯著差異。(2)對甲的區(qū)間估計：>a<-c(115,116,98,83)>t.test(a)OneSamplet-testdata: at=13.1341,df=3,p-value=0.0009534alternativehypothesis:truemeanisnotequalto095percentconfidenceinterval:78.04264127.95736sampleestimates:meanofx103或者用這個命令更簡單：>attach(lamp)>t.test(X[A==1])乙的均值估計為111，95%置信區(qū)間為99.59932,122.40068。丙的均值估計為86，95%置信區(qū)間為70.08777,101.91223。(3)多重檢驗：>attach(lamp)P值不做調整：>pairwise.t.test(X,A,p.adjust.method="none")PairwisecomparisonsusingttestswithpooledSDdata: XandA1 220.351-30.0660.013----宋停云與您分享--------宋停云與您分享----PHolm調整：Pvalueadjustmentmethod:none>pairwise.t.test(X,A,p.adjust.method="holm",data)PairwisecomparisonsusingttestswithpooledSDdata: XandA1 220.35-30.130.04Pvalueadjustmentmethod:holm不論采取哪種方法，都可看出乙和丙有顯著差異。Ex7.2(1)>lamp<-data.frame(X=c(20,18,18,17,15,16,13,18,22,17,26,19,26,28,23,25,24,25,18,22,27,24,12,14),A=factor(rep(1:4,c(10,6,6,2))))>lamp.aov<-aov(X~A,data=lamp);summary(lamp.aov)DfSumSqMeanSqFvalue Pr(>F)A 3 351.7 117.24 15.112.28e-05***Residuals 20 155.2 7.76P值小于0.05，可認為四個廠生產的產品的變化率有顯著差異。(2)>attach(lamp)P值不做調整：>pairwise.t.test(X,A,p.adjust.method="none")PairwisecomparisonsusingttestswithpooledSDdata: XandA1 2 328.0e-05- -30.000530.47666-40.054906.1e-050.00020Pvalueadjustmentmethod:nonePHolm調整：>pairwise.t.test(X,A,p.adjust.method="holm")PairwisecomparisonsusingttestswithpooledSD----宋停云與您分享--------宋停云與您分享----data: XandA1 2 320.00040- -30.001580.47666-40.109790.000360.00079Pvalueadjustmentmethod:holm由此可得，除了A1和A4，A2和A3這兩組的差異不顯著外，其他組合的差異都很顯著。Ex7.3>lamp1<-data.frame(X=c(30,27,35,35,29,33,32,36,26,41,33,31,43,45,53,44,51,53,54,37,47,57,48,42,82,66,66,86,56,52 ep(1:3,c(12,12,12))))>attach(lamp1)正態(tài)性檢驗：>shapiro.test(X[A==1])Shapiro-Wilknormalitytestdata: X[A==1]W=0.9731,p-value=0.9407>shapiro.test(X[A==2])Shapiro-Wilknormalitytestdata: X[A==2]W=0.9708,p-value=0.9193>shapiro.test(X[A==3])Shapiro-Wilknormalitytestdata: X[A==3]W=0.9371,p-value=0.4613數(shù)據(jù)在三種水平下均是正態(tài)的。方差齊性檢驗：>bartlett.test(X~A,data=lamp1)Bartletttestofhomogeneityofvariancesdata: XbyABartlett'sK-squared=12.139,df=2,p-value=0.002312----宋停云與您分享--------宋停云與您分享----P值小于0.05，認為各組方差不等。Ex7.4>lamp<-data.frame(X=c(2.79,2.69,3.11,3.47,1.77,2.44,2.83,2.52,3.83,3.15,4.70,3.97,2.03,2.87,3.65,5.09,5.41,3.47,4.92,4.07,2.18,3.13,3.77,4.26),g=factor(rep(1:3,c(8,8,8))))先進行正態(tài)性和方差齊性檢驗以選擇使用方差分析aov()還是KW檢驗kruskal.test()。正態(tài)性檢驗：>attach(lamp)>shapiro.test(X[g==1])Shapiro-Wilknormalitytestdata: X[g==1]W=0.9659,p-value=0.8638>shapiro.test(X[g==2])Shapiro-Wilknormalitytestdata: X[g==2]W=0.983,p-value=0.9763>shapiro.test(X[g==3])Shapiro-Wilknormalitytestdata: X[g==3]W=0.99,p-value=0.9951三組數(shù)據(jù)都服從正態(tài)分布。方差齊性檢驗：>bartlett.test(X~g,data=lamp)Bartletttestofhomogeneityofvariancesdata: XbygBartlett'sK-squared=3.4559,df=2,p-value=0.1776p值大于5)或者W檢驗均可。>kruskal.test(X~g,data=lamp)Kruskal-Wallisranksumtestdata: Xbyg----宋停云與您分享--------宋停云與您分享----Kruskal-Wallischi-squared=7.9322,df=2,p-value=0.01895>lamp.aov<-aov(X~g,data=lamp);summary(lamp.aov)DfSumSqMeanSqFvaluePr(>F)g 2 6.437 3.218 4.2840.0275*Residuals 2115.776 0.751兩種方法得到的p值均小于0.05，可認為3種不同處理的誘導作用不同。Ex7.5根據(jù)題意，適用配伍組設計的Friedman秩和檢驗。>lamp<-data.frame(X=c(23.1,57.6,10.5,23.6,11.9,54.6,21.0,20.3,22.7,53.2,9.7,19.6,13.8,47.1,13.6,23.6,22.5,53.7,10.8,21.1,13.7,39.2,13.7,16.3,22.6,53.1,8.3,21.6,13.3,37.0,14.8,14.8),g=gl(4,8),b=gl(8,1,32))#其中g表示group，b表示block。>friedman.test(X~g|b,data=lamp)Friedmanranksumtestdata: XandgandbFriedmanchi-squared=6.45,df=3,p-value=0.09166P值大于0.05，尚不能拒絕原假設。Ex7.6(1)>qua<-data.frame(x=c(4.6,4.3,6.1,6.5,6.8,6.4,6.3,6.7,3.4,3.8,4.0,3.8,4.7,4.3,3.9,3.5,6.5,7.0),a=gl(3,6,18),b=gl(3,2,18))>qua.aov<-aov(x~a+b+a:b,data=qua);summary(qua.aov)DfSumSqMeanSqFvalue Pr(>F)a2 3.9741.987 26.690.000164***b2 4.4412.221 29.830.000107***a:b421.1595.290 71.068.34e-07***Residuals9 0.074兩種因素以及其交互作用對產品質量的影響都很顯著。(2)最優(yōu)條件為A3和B3組合。>t.test(c(6.5,7.0))OneSamplet-testdata: c(6.5,7)t=27,df=1,p-value=0.02357alternativehypothesis:truemeanisnotequalto095percentconfidenceinterval:3.5734499.926551sampleestimates:meanofx6.75----宋停云與您分享--------宋停云與您分享----點估計為6.75，區(qū)間估計為3.573449，9.926551(3)雙因素方差分析的多重比較？不會...Ex7.7正交試驗的方差分析。L9(3^4)正交表。>pro<-data.frame(Y=c(62.925,57.075,51.6,55.05,58.05,56.55,63.225,50.7,54.45),A=gl(3,3),