烏龜性別與溫度的logit模型-數(shù)學(xué)模型

1、摘要本文針對(duì)幼龜性別溫度關(guān)系進(jìn)行建立模型，在建模時(shí)建立溫度與雄龜比例的模型，由于雄龜比例只可能在【0,1】之間取值，所以不能建立一般的統(tǒng)計(jì)回歸模型，所以可以建立Logit模型。利用MATLAB統(tǒng)計(jì)工具箱中的命令glmfit求解。再求出各個(gè)系數(shù)估計(jì)值后用EXCEL表格進(jìn)行處理算出每個(gè)溫度下雄龜比例的估計(jì)值，對(duì)模型進(jìn)行驗(yàn)證，為了提高模型的擬合效果可以在Logit模型中添加t的高次方，并用利用MATLAB統(tǒng)計(jì)工具箱中的命令regress求解,求出各個(gè)系數(shù)估計(jì)值，置信區(qū)間和模型相關(guān)系數(shù)。為了能直觀的比較模型擬合效果可以利用MATLAB中的plot工具會(huì)出（t,logit）的散點(diǎn)圖，并用利用Tools下

2、Basic Fitting工具找出一條擬合效果最好的曲線，并得出其對(duì)應(yīng)的系數(shù)。最后可以通過(guò)logit基本模型得出當(dāng)溫度為27.7329度幼龜比例為1:1，而且溫度每增加一度雄龜與雌龜比例擴(kuò)大9.12484倍。關(guān)鍵詞：logit模型，置信區(qū)間，相關(guān)系數(shù)一、問(wèn)題提出人類的性別是由基因決定的，烏龜?shù)男詣e主要有什么因素決定的呢？科學(xué)研究表明，決定幼龜性別的最關(guān)鍵的因素是烏龜孵化時(shí)的溫度。為了研究溫度是如何影響幼龜?shù)拇菩郾壤?，美?guó)科學(xué)家對(duì)某一類烏龜?shù)姆趸^(guò)程做了實(shí)驗(yàn)，試驗(yàn)在五個(gè)不同恒定溫度下進(jìn)行，每個(gè)溫度下分別觀察3批烏龜?shù)暗姆趸^(guò)程，得到的數(shù)據(jù)如下：溫度烏龜?shù)皞€(gè)數(shù)雄龜個(gè)數(shù)雌龜個(gè)數(shù)雄龜比例27.2101

3、910%8080%91811.1%27.7107370%64266.7%86275%28.313130100%96366.7%87187.5%28.4107370%85362.5%97277.8%29.91110190.9%880100%990100%二、基本假設(shè)假設(shè)1；幼龜性別只與溫度有關(guān)三、符號(hào)說(shuō)明符號(hào)意義單位備注t烏龜孵化時(shí)的溫度P（x）雄龜比例S烏龜?shù)翱倲?shù)個(gè)X雄龜個(gè)數(shù)個(gè)C雌龜個(gè)數(shù)個(gè)四、問(wèn)題分析 在本題由于是求溫度與性別比例的模型，在數(shù)據(jù)表中每個(gè)溫度都記錄了三批烏龜，所以首先要對(duì)數(shù)據(jù)進(jìn)行處理，利用excel算出每個(gè)溫度對(duì)應(yīng)下烏龜?shù)目倲?shù)，雄龜總個(gè)數(shù)，雌龜總個(gè)數(shù)和雄龜比例。溫度烏龜?shù)皞€(gè)數(shù)雄龜

4、個(gè)數(shù)雌龜個(gè)數(shù)雄龜比例27.2272257.41%27.72417770.83%28.33026486.67%28.42719870.37%29.92827196.43%為了使得運(yùn)算更簡(jiǎn)單可以把溫度進(jìn)行預(yù)處理，把27.2看成0，以此類推可得溫度烏龜?shù)皞€(gè)數(shù)雄龜個(gè)數(shù)雌龜個(gè)數(shù)雄龜比例0272257.41%0.52417770.83%1.13026486.67%1.22719870.37%2.72827196.43%為了更直觀觀察其回歸關(guān)系，利用MATLAB繪制出散點(diǎn)圖。從圖中可以看出回歸曲線是一條近于3次樣條的多項(xiàng)式回歸曲線，其回歸模型為 （1）然而在這個(gè)問(wèn)題中（1）是回歸方程中P(x)的取值不一定在

5、0,1中，即使P（x）取值在0,1中，有意在給定t是，誤差項(xiàng)也只能取0,1兩個(gè)值，顯然不具有正態(tài)性，而且的方差依賴于與t，具有異方差性，這些都違反了普通回歸分析的前提條件，因此，該題不能用用普通回歸分析。 由于P（x）在0,1之間取值，可以使用Logit模型。五、模型的建立與求解5.1模型的建立5.2模型的求解Logit模型是一種廣義線性模型，可利用MATLAB統(tǒng)計(jì)工具箱總的命令glmfit求解。參數(shù)參數(shù)估計(jì)值標(biāo)準(zhǔn)差B0-1017830.3739B12.21100.4309所以估計(jì)值為即所以當(dāng)幼龜比例為1:1時(shí)計(jì)算出溫度為27.7329度令odd(t)為雄龜與雌龜比例，故有當(dāng)溫度增加1度時(shí)od

6、d比為于是由于b1=2.2110,所以溫度每增加一度時(shí)雄龜與雌龜比例增加到原來(lái)的9.12484倍5.3結(jié)果的分析及驗(yàn)證對(duì)模型各個(gè)溫度進(jìn)行驗(yàn)證溫度烏龜?shù)皞€(gè)數(shù)雄龜個(gè)數(shù)雌龜個(gè)數(shù)雄龜比例雄龜估計(jì)值0272257.41%0.2353570.52417770.83%0.4818121.13026486.67%0.7779551.22719870.37%0.8138012.72827196.43%0.9917685.3模型改進(jìn)從以上結(jié)果可知擬合偏差太大，不適合于做為最終結(jié)果。由于模型的右端是溫度t的線性函數(shù)，可以考慮加入t的二次項(xiàng)后，看是否能提高模型的擬合程度。即考慮模型為執(zhí)行以下程序 b2,dev2=gl

7、mfit(t t.2,x s,binomial,logit);b2,pval=1-chi2cdf(dev-dev2,1)b2 = -1.6582 3.7840 -0.7745pval =0.0304計(jì)算出b0,b1,b2的估計(jì)值為-1.6582、 3.7840、 -0.7745所以模型為溫度烏龜?shù)皞€(gè)數(shù)雄龜個(gè)數(shù)雌龜個(gè)數(shù)雄龜比例雄龜估計(jì)值0272257.41%0.1600030.52417770.83%0.5121781.13026486.67%0.8273631.22719870.37%0.8541242.72827196.43%0.948461由以上表可知擬合偏差減小，由此可知，加入高次方后可

8、以提高擬合偏差，所以為了進(jìn)一步提高擬合效果，可以先計(jì)算出logit對(duì)應(yīng)的值并畫出并繪出（t,logit）的散點(diǎn)圖溫度烏龜?shù)皞€(gè)數(shù)雄龜個(gè)數(shù)雌龜個(gè)數(shù)雄龜比例Logit(t)0272257.41%-2.525350.52417770.83%0.88714181.13026486.67%1.872090671.22719870.37%0.86497962.72827196.43%3.296252再利用Tools下Basic Fitting工具找出一條擬合效果最好的曲線，并得出其對(duì)應(yīng)的系數(shù)可知模型為個(gè)系數(shù)為所以模型為當(dāng)所以幼龜比例為1:1時(shí)P(x)=0.5.所以logit=0;即可以得出當(dāng)溫度為27.57

9、932852、28.4627667和29.854698時(shí)幼龜比例為1:1；再次各個(gè)數(shù)據(jù)進(jìn)行驗(yàn)證可得溫度烏龜?shù)皞€(gè)數(shù)雄龜個(gè)數(shù)雌龜個(gè)數(shù)雄龜比例雄龜估計(jì)值0272257.41%0.07410.52417770.83%0.70831.13026486.67%0.866661.22719870.37%0.70352.72827196.43%0.9633由以上結(jié)果可以看出擬合程度很好。在用MATLAB統(tǒng)計(jì)工具中的命令regress求解得知相關(guān)系數(shù)R為1；擬合效果很好。六、參考文獻(xiàn)數(shù)學(xué)模型 第四版 高等教育出版社七、附錄模型1求解程序 t=0 0.5 1.1 1.2 2.7; x=2 17 26 19 27;

10、 s=27 24 30 27 28; proport=x./s; b,dev,stats=glmfit(t,x s,binomial,logit); logitfit=glmval(b,t,logit); plot(t,proport,o,t,logitfit,r-); xlabel(t);ylabel(proportion of x) b,bi=stats.se,devb = -1.1783 2.2110bi= 0.3739 0.4309dev = 14.8629模型1驗(yàn)證程序 x=0 1;0.5 1;1.1 1;1.2 1;2.7 1x = 0 1.0000 0.5000 1.0000 1

11、.1000 1.0000 1.2000 1.0000 2.7000 1.0000 y=-2.52535 0.8871418 1.87209067 0.8649796 3.296252;alpha=0.05;b,bint,r,rint,stats=regress(y,x,alpha)b = 1.8372 -1.1419bint = -0.0667 3.7411 -3.8599 1.5760r = -1.3834 1.1105 0.9931 -0.1978 -0.5223rint = -2.7014 -0.0654 -2.0210 4.2419 -2.6144 4.6005 -4.4115 4.0

12、160 -2.1677 1.1230stats =0.7587 9.4311 0.0545 1.4817模型2的驗(yàn)證程序x=0 0 1;0.5 0.52 1;1.12 1.1 1;1.22 1.2 1;2.72 2.7 1x = 0 0 1.0000 0.5000 0.2500 1.0000 1.2100 1.1000 1.0000 1.4400 1.2000 1.0000 7.2900 2.7000 1.0000 y=-2.52535 0.8871418 1.87209067 0.8649796 3.296252;alpha=0.05;b,bint,r,rint,stats=regress(

13、y,x,alpha)b = -0.3797 2.7004 -1.1635bint = -3.9912 3.2317 -7.4018 12.8026 -6.3645 4.0374r = -1.3618 1.5654 0.5246 -0.6652 -0.0631rint = -3.5130 0.7893 -1.3332 4.4641 -6.0070 7.0562 -6.7698 5.4395 -0.3836 0.2575stats = 0.7271 2.6642 0.2729 2.5134繪制（t,P(x)）散點(diǎn)圖程序 t=0 0.5 1.1 1.2 2.7;p=0.0741 0.7083 0.8

14、667 0.7037 0.9643; plot(t,p,r*)繪制（t,logit）散點(diǎn)圖程序 t=0 0.5 1.1 1.2 2.7;p=-2.52535 0.8871418 1.87209067 0.8649796 3.296252; plot(t,p,r*)用第三模型計(jì)算幼龜比例為1:1時(shí)溫度程序 y=4.1581 -15.113 9.9691 0 5.9619 -2.5254y = 4.1581 -15.1130 9.9691 0 5.9619 -2.5254 roots(y)ans = 2.6755 1.2628 -0.3415 + 0.5977i -0.3415 - 0.5977i

15、 0.3793 vpa( roots(y),8) ans = 2.6754698 1.2627667 0.59774412*i-0.34148611 -0.59774412*i-0.34148611 0.37932852第三模型各個(gè)溫度對(duì)應(yīng)雄龜比例結(jié)果驗(yàn)證程序 Y=(t)-2.5254+5.9619*t+9.9691*t3-15.113*t4+4.1581*t5Y = (t)-2.5254+5.9619*t+9.9691*t3-15.113*t4+4.1581*t5 f=(t)exp(Y(t)/(1+exp(Y(t)f = (t)exp(Y(t)/(1+exp(Y(t) f(0)ans = 0

16、.0741 f(0.5)ans = 0.7083 f(1.1)ans = 0.8666 f(1.2)ans = 0.7035 f(2.7)ans =0.9633用regress驗(yàn)證模型程序x=0 0 0 0 0 1;0.55 0.54 0.53 0.52 0.5 1;1.15 1.14 1.13 1.12 1.1 1;1.25 1.24 1.23 1.22 1.2 1;2.75 2.74 2.73 2.72 2.7 1x = 0 0 0 0 0 1.0000 0.0313 0.0625 0.1250 0.2500 0.5000 1.0000 1.6105 1.4641 1.3310 1.2100 1.1000 1.0000 2.4883 2.0736 1.7280 1.4400 1.2000 1.0000 143.4891 53.1441 19.6830 7.2900 2.7000 1.0000 y=-2.52535 0.8871418 1.87209067 0.8649796 3.296252;alpha=0.05;b,bint,r,rint,stats=regress(y,x,alpha)Warning: X is rank d