高教社杯全國大學生數學建模競賽-資料下載頁

　　

【正文】 (:,[4,5])。%低優度處理ra=tiedrank(a)。%對每個指標值分別編秩，即對a的每一列分別編秩[n,m]=size(a)。%計算矩陣的維數b=[,]。%權重w=b/sum(b)。%歸一化處理WRSR=ra*w39。/n。[WRSR i]=sort(WRSR,39。ascend39。)%升序排列WRSRn=32。p=[1:n]/n。%計算累計頻率p(end)=11/(4*n)。%修正最后一個累計頻率，最后一個累積頻率按11/(4n)估計p=[]。 X=[ones(n,1),p]。%構造一元線性回歸分析的數據矩陣[ab,abint,r,rint,stats]=regress(WRSR,X)%一元線性回歸分析WRSRfit=ab(1)+ab(2)*probit%計算WRSR的估計值程序2：aw=[ 1848 1893 1987 2036 ]；%導入數據a=aw(1:end,:)；%提取指標數據a(:,[4,4])=a(:,[4,4])；%低優度處理ra=tiedrank(a)%對每個指標值分別編秩，即對a的每一列分別編秩[n,m]=size(a)%計算矩陣的維數b=[,]%權重w=b/sum(b)%歸一化處理WRSR=ra*w39。/n[WRSR i]=sort(WRSR,39。ascend39。)%升序排列WRSR 程序3：灰色預測的程序x = [11]39。 sizexd2 = size(x,2)。%求數組長度k=0。for y1=x k=k+1。 if k1 x1(k)=x1(k1)+x(k)。 %累加生成 z1(k1)=*(x1(k)+x1(k1))。 %z1維數減1，用于計算B yn1(k1)=x(k)。 else x1(k)=x(k)。 endend%x1,z1,k,yn1sizez1=size(z1,2)。%size(yn1)。z2 = z139。 z3 = ones(1,sizez1)39。YN = yn139。 %轉置%YNB=[z2 z3]。au0=inv(B39。*B)*B39。*YN。au = au039。%B,au0,auafor = au(1)。ufor = au(2)。ua = au(2)./au(1)。%afor,ufor,ua %輸出預測的 a u 和 u/a的值constant1 = x(1)ua。afor1 = afor。x1t1 = 39。x1(t+1)39。estr = 39。exp39。tstr = 39。t39。leftbra = 39。(39。rightbra = 39。)39。%constant1,afor1,x1t1,estr,tstr,leftbra,rightbrastrcat(x1t1,39。=39。,num2str(constant1),estr,leftbra,num2str(afor1),tstr,rightbra,39。+39。,leftbra,num2str(ua),rightbra)%輸出時間響應方程%******************************************************nfinal = sizexd21 + 1。%決定預測的步驟數5 這個步驟可以通過函數傳入%nfinal = sizexd2 1 + 1。%預測的步驟數 1for k3=1:nfinal x3fcast(k3) = constant1*exp(afor1*k3)+ua。end %x3fcast%一次擬合累加值for k31=nfinal:1:0 if k311 x31fcast(k31+1) = x3fcast(k31)x3fcast(k311)。 else if k310 x31fcast(k31+1) = x3fcast(k31)x(1)。 else x31fcast(k31+1) = x(1)。 end end endx31fcastfor i=1:32 x4(i)=x31fcast(i)。end %二次擬合預測值%***精度檢驗p C************//////////////////////////////////k5 = 0。for y5 = x k5 = k5 + 1。 if k5 sizexd2 else err1(k5) = x(k5) x31fcast(k5)。 endend%err1%絕對誤差 xavg = mean(x)。%xavg%x平均值err1avg = mean(err1)。%err1avg%err1平均值k5 = 0。s1total = 0 。for y5 = x k5 = k5 + 1。 if k5 sizexd2 else s1total = s1total + (x(k5) xavg)^2。 endends1suqare = s1total ./ sizexd2。s1sqrt = sqrt(s1suqare)。%s1suqare,s1sqrt%s1suqare 殘差數列x的方差 s1sqrt 為x方差的平方根S1k5 = 0。s2total = 0 。for y5 = x k5 = k5 + 1。 if k5 sizexd2 else s2total = s2total + (err1(k5) err1avg)^2。 endends2suqare = s2total ./ sizexd2。%s2suqare 殘差數列err1的方差S2Cval = sqrt(s2suqare ./ s1suqare)。Cval%nnn = * s1sqrt%Cval C檢驗值k5 = 0。pnum = 0 。for y5 = x k5 = k5 + 1。 if abs( err1(k5) err1avg ) * s1sqrt pnum = pnum + 1。 %ppp = abs( err1(k5) err1avg ) else endendpval = pnum ./ sizexd2。pval%p檢驗值%arr1 = x41fcast(1:6)%預測結果為區間范圍 預測步長和數據長度可調整程序參數進行改進t=[79,80,81,82,83,84,85,96,87,88,89,90,91,92,93,94,95,96,97,98,99,00,01,02,03,04,05,06,07,08,09,10]。plot(t,x,39。go39。,t,x4,39。rd39。)t=33y1=((x(1)ua)*exp((afor1*t))+ua)(((x(1)ua)*exp((afor1*(t1)))+ua))t=34y2=((x(1)ua)*exp((afor1*t))+ua)(((x(1)ua)*exp((afor1*(t1)))+ua))程序4X2=[ 1047 ]。X3=[ ]。X4=[ ]。X5=[ 1633 ]。Y2=[ 1028 ]。Y3=[ 1086 ]。Y4=[ 1112 ]。Y5=[ 1581 318]。P2=ranksum(X2,Y2)P3=ranksum(X3,Y3)P4=ranksum(X4,Y4)P5=ranksum(X5,Y5)