《乌龟性别与温度的logit模型.docx》由会员分享，可在线阅读，更多相关《乌龟性别与温度的logit模型.docx（15页珍藏版）》请在三一办公上搜索。
1、数学建模第四次作业第一小组乌龟性别的决定因素问题（习题1010）钟杰 20107122赵新 20107134卢昌 20107127摘要本文针对幼龟性别温度关系进行建立模型，在建模时建立温度与雄龟比例的 模型，由于雄龟比例只可能在【0,1】之间取值，所以不能建立一般的统计回归 模型，所以可以建立Logit模型。利用MATLAB统计工具箱中的命令glmfit求解。 再求出各个系数估计值后用EXCEL表格进行处理算出每个温度下雄龟比例的估 计值，对模型进行验证，为了提高模型的拟合效果可以在Logit模型中添加t 的高次方，并用利用MATLAB统计工具箱中的命令regress求解，求出各个系数估 计值
2、，置信区间和模型相关系数。为了能直观的比较模型拟合效果可以利用 MATLAB中的plot工具会出（t,logit）的散点图，并用利用Tools下Basic Fitting 工具找出一条拟合效果最好的曲线，并得出其对应的系数。最后可以通过logit基本模型得出当温度为27.7329度幼龟比例为1:1，而且温度每增加 一度雄龟与雌龟比例扩大9.12484倍。关键词：logit模型，置信区间，相关系数、问题提出人类的性别是由基因决定的，乌龟的性别主要有什么因素决定的呢？科学研究表 明，决定幼龟性别的最关键的因素是乌龟孵化时的温度。为了研究温度是如何影 响幼龟的雌雄比例，美国科学家对某一类乌龟的孵化过
3、程做了实验，试验在五个 不同恒定温度下进行，每个温度下分别观察3批乌龟蛋的孵化过程，得到的数据 如下：温度乌龟蛋个数雄龟个数雌龟个数雄龟比例27.2101910%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%二、基本假设假设1；幼龟性别只与温度有关三、符号说明符号意义单位备注t乌龟孵化时的温度CP（X）雄龟比例S乌龟蛋总数个X雄龟个数个C雌龟个数个四、问题分析在本题由于是求温度与性别比例的模型，在数
4、据表中每个温度都记录了三批乌 龟，所以首先要对数据进行处理，利用excel算出每个温度对应下乌龟的总数， 雄龟总个数，雌龟总个数和雄龟比例。温度乌龟蛋个数雄龟个数雌龟个数雄龟比例27.2272257.41%27.72417770.83%28.33026486.67%28.42719870.37%29.92827196.43%为了使得运算更简单可以把温度进行预处理，把27.2看成0，以此类推可得温度乌龟蛋个数雄龟个数雌龟个数雄龟比例0272257.41%0.52417770.83%1.13026486.67%1.22719870.37%2.72827196.43%为了更直观观察其回归关系，利用M
5、ATLAB绘制出散点图。从图中可以看出回归曲线是一条近于3次样条的多项式回归曲线，其回归模型为P（x） - b0 + b1*t + b2*tA 2 + b3*t A3 + s（1）然而在这个问题中（1）是回归方程中P（x）的取值不一定在0,1中，即使P（x）取值在0,1 中，有意在给定t是，误差项也只能取0,1两个值，显然不具有正态性，而且的方差依 赖于与t，具有异方差性，这些都违反了普通回归分析的前提条件，因此，该题不能用用普 通回归分析。由于P (x)在0,1之间取值，可以使用Logit模型。五、模型的建立与求解5.1模型的建立T、P 、，Logit(P(t)=Ln(二)=b 0 + b1
6、* t 1-P(t)P,v、()1 + eb0+bi*t5.2模型的求解File Edit View Insert Tools Desktop Window HelpLhFigure 1Logit模型是一种广义线性模型，可利用MATLAB统计工具箱总的命令glmfit求解。Q a d | M、上档1闵l业，也口田|参数参数估计值标准差B0B1-1017832.21100.37390.4309b 0 = 1.1783 所以估计值为b = 2.2110,、, P e1.1783+2.2110气Logit(P ) = Ln(T) ) = 1.1783 + 2.2110* tPgf(X)1 + e1.
7、1783+2.2110*t所以当幼龟比例为1:1时计算出温度为27.7329度 令odd(t)为雄龟与雌龟比例，故有odd (t) =eb 0+b1t当温度增加1度时odd比为odd (t +1) eb 0+b1( t+1)odd (t)eb 0+b1t=eb1于是odd (t +1) = eb1odd (t)由于b1=2.2110，所以温度每增加一度时雄龟与雌龟比例增加到原来的9.12484倍5.3结果的分析及验证对模型各个温度进行验证温度乌龟蛋个数雄龟个数雌龟个数雄龟比例雄龟估计值0272257.41%0.2353570.52417770.83%0.4818121.13026486.67%
8、0.7779551.22719870.37%0.8138012.72827196.43%0.9917685.3模型改进从以上结果可知拟合偏差太大，不适合于做为最终结果。由于模型的右端是温度 t的线性函数，可以考虑加入t的二次项后，看是否能提高模型的拟合程度。即 考虑模型为D, log it(p) = ln( D = b0 + b1t + b2t2 1 - p (t)执行以下程序 b2,dev2=glmfit(t t.A2,x s,binomial,logit);b2,pval=1-chi2cdf(dev-dev2,1)-1.65823.7840-0.7745pval =0.0304计算出 b0
9、,b1,b2 的估计值为-1.6582、3.7840、-0.7745所以模型为log it（p ） = ln（ P（t） = -1.6582 + 3.7840t - 0.7745t2（t）1 - p（ t）温度00.5U 12 2.7乌龟蛋个雄龟个数数27224302728雌龟个数雄龟比例12261227257L817.41%70.83%86.67%70.37%96.43%雄龟估计值0.1600030.5121780.8273630.8541240.948461由以上表可知拟合偏差减小，由此可知，加入高次方后可以提高拟合偏差，所以为了进步提高拟合效果，可以先计算出logit对应的值并画出并绘出
10、（t,logit）的散点图温度乌龟蛋个数雄龟个数雌龟个数雄龟比例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工具找出一条拟合效果最好的曲线，并得出其对应的系数可知模型为/、 p,、log it(p(t) = ln(J = bo + b1t + b2t 2 + b3t 3 + b 4t 4 + b 5t 5(t)1 -p (t)个系数为b0 = -2.5254外
11、=5.9619b 2 = 0b3 = 9.9691b 4 = -15.113b5 = 4.1581所以模型为log it(p ) = ln(-) = -2.5254 + 5.9619t + 9.9691t3 -15.113t4 + 4.1581t5(t)1 - P (t)当所以幼龟比例为1:1时P(x)=0.5.所以logit=0;即log it(p )=血(一) = -2.5254 + 5.9619t + 9.9691t3 - 15.113t4 + 4.1581t5 = 0(t)1 - p (t)可以得出当温度为27.57932852、28.4627667和29.854698时幼龟比例为1:
12、1；再次各个数据进行验证可得温度乌龟蛋个数雄龟个数雌龟个数雄龟比例雄龟估计值0272257.41%0.07410.52417770.83%0.70831.13026486.67%0.866661.22719870.37%0.70352.72827196.43%0.9633由以上结果可以看出拟合程度很好。在用MATLAB统计工具中的命令regress求解得知相关系数R为1；拟合效果很好。六、参考文献数学模型第四版高等教育出版社七、附录模型1求解程序 t=0 0.5 1.1 1.2 2.7; x=2 17 26 19 27; s=27 24 30 27 28; proport=x./s; b,de
13、v,stats=glmfit(t,x s,binomial,Togit); logitfit二glmval(b,t,logit); plot(t,proport,o,t,logitfit,r-); xlabel(t);ylabel(proportion of x) b,bi=stats.se,devb =-1.17832.2110bi=0.37390.4309dev =14.8629模型1验证程序 x=0 1;0.5 1;1.1 1;1.2 1;2.7 101.00000.50001.00001.10001.00001.20001.00002.70001.0000 y=-2.52535 0.8
14、871418 1.87209067 0.8649796 3.296252;alpha=0.05;b,bint,r,rint,stats=regress(y,x,alpha)b =1.8372-1.1419bint =-0.06673.7411-3.85991.5760r =-1.38341.11050.9931-0.1978-0.5223rint =-2.7014-0.0654-2.02104.2419-2.61444.6005-4.41154.0160-2.16771.1230stats =0.75879.43110.05451.4817模型2的验证程序x=0 0 1;0.5 0.5A2 1
15、;1.1A2 1.1 1;1.2A2 1.2 1;2.7A2 2.7 1 x =001.00000.50000.25001.00001.21001.10001.00001.44001.20001.00007.29002.70001.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 =-0.37972.7004-1.1635bint =-3.99123.2317-7.401812.8026-6.36454.0374r =-1.3618
16、1.56540.5246-0.6652-0.0631rint =-3.51300.7893-1.33324.4641-6.00707.0562-6.76985.4395-0.38360.2575stats =0.72712.66420.27292.5134绘制(t,P(x)散点图程序 t=0 0.5 1.1 1.2 2.7;p=0.0741 0.7083 0.8667 0.7037 0.9643; plot(t,p,r*)绘制(t,logit)散点图程序 t=0 0.5 1.1 1.2 2.7;p=-2.52535 0.8871418 1.87209067 0.8649796 3.296252
17、; plot(t,p,r*)用第三模型计算幼龟比例为1:1时温度程序 y=4.1581 -15.113 9.9691 0 5.9619 -2.5254y =4.1581 -15.11309.969105.9619-2.5254 roots(y)ans =2.67551.2628-0.3415 + 0.5977i-0.3415 - 0.5977i0.3793vpa( roots(y),8)ans =2.67546981.26276670.59774412*i-0.34148611-0.59774412*i-0.341486110.37932852 第三模型各个温度对应雄龟比例结果验证程序 Y=(
18、t)-2.5254+5.9619*t+9.9691*tA3-15.113*tA4+4.1581*tA5Y =(。-2.5254+5.9619气+9.9691气人3-15.113气人4+4.1581气人5 f=(t)exp(Y(t)/(1+exp(Y(t)f =(t)exp(Y(t)/(1+exp(Y(t) f(0)ans =0.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验证模型程序x=0 0 0 0 0 1;0.5人5 0.5人4 0.5人3 0.5人2 0.5 1;1
19、.1人5 1.1人4 1.1人3 1.1人2 1.1 1；1.2人5 1.2人4 1.2人31.2A2 1.2 1;2.7A5 2.7A4 2.7A3 2.7A2 2.7 1000001.00000.03130.06250.12500.25000.50001.00001.61051.46411.33101.21001.10001.00002.48832.07361.72801.44001.20001.0000143.489153.144119.68307.29002.70001.0000 y=-2.52535 0.8871418 1.87209067 0.8649796 3.296252;al
20、pha=0.05;b,bint,r,rint,stats=regress(y,x,alpha)Warning: X is rank deficient to within machine precision. In regress at 82b =4.1581-15.11259.969105.9619-2.5254 bint =NaN NaNNaN NaNNaN NaN00NaN NaNNaN NaN r =1.0e-013 *0.3775-0.05110.03110.01000.4263 rint =NaN NaNNaN NaNNaN NaNNaN NaNNaN NaN stats =1 NaN NaN NaN
