2017-06-15 56 views
-3

我正在通過自己編寫一個簡單的腳本來計算PDSI,但是當我運行腳本時,我在DT中的輸出中獲得了不同的值。我手動完成了所有這些計算,並在Excel中進行了驗證,並驗證了兩者,並得到了正確的值。在IF中的語句輸出不正確

這是腳本的摘錄,我沒有任何錯誤,所以我不知道我的錯誤在哪裏。

Tb.PDSI<-as.data.table(matrix(ncol = 22, nrow = 190)) 
colnames(Tb.PDSI)<-c("Anho","Mes","Z","Uw.Z","Ud.Z","V","Ze","Q","P","Ms.X1","AX.1","X1","Ms.X2","AX.2","X2","Ms.X3","AX.3","X3","O","IndeX.f","Dur.Mes","Z3") 


Tb.PDSI<-rbind(data.frame(Anho=1999,Mes=12,Z=0,Uw.Z=0,Ud.Z=0,V=0,Ze=0,Q=0,P=0,Ms.X1=0,AX.1=0,X1=0,Ms.X2=0,AX.2=0, 
        X2=0,Ms.X3=0,AX.3=0,X3=0,O=0,IndeX.f=0,Dur.Mes=0, Z3=0),Tb.PDSI, fill=F) 

mb<-36.03 
mmb<-0.07 
b<-33.58 
mb0<-18.01 
m2<-1.22 

# ...(3) Index Table 
for(i in 2:6){ 
    for(j in 4:21){ 

    # ... Uw.Z 
    if (j == 4){ 
     Tb.PDSI[[i,j]]<-ifelse(Tb.PDSI$O[i-1]==0,0, 
          ifelse(Tb.PDSI$X3[i-1]<0, 
            ifelse(Tb.PDSI$Uw.Z[i-1]==0, 
              ifelse(Tb.PDSI$Z[i]>-m2,Tb.PDSI$Z[i]+m2,0),Tb.PDSI$Z[i]+m2),0)) 
    } else if(j == 5){ # ... Ud.Z 
     Tb.PDSI[[i,j]]<-ifelse(Tb.PDSI$O[i-1]==0,0, 
          ifelse(Tb.PDSI$X3[i-1]>0, 
            ifelse(Tb.PDSI$Ud.Z[i-1]==0, 
              ifelse(Tb.PDSI$Z[i]<m2,Tb.PDSI$Z[i]-m2,0),Tb.PDSI$Z[i]-m2),0)) 
    } else if(j==6){  # ... V 
     Tb.PDSI[[i,j]]<-ifelse(Tb.PDSI$O[i-1]==0,0, 
          ifelse(Tb.PDSI$O[i-1]==1, 
            ifelse(Tb.PDSI$Ud.Z[i]+Tb.PDSI$V[i-1]>0 && Tb.PDSI$P[i-1]<100,0, 
              ifelse(Tb.PDSI$V[i-1]==0 || Tb.PDSI$P[i-1]==100,Tb.PDSI$Ud.Z[i],Tb.PDSI$V[i-1]+Tb.PDSI$Ud.Z[i])), 
            ifelse(Tb.PDSI$Uw.Z[i]+Tb.PDSI$V[i-1]<0 && Tb.PDSI$P[i-1]<100,0, 
              ifelse(Tb.PDSI$V[i-1]==0 || Tb.PDSI$P[i-1]==100,Tb.PDSI$Uw.Z[i],Tb.PDSI$Uw.Z[i]+Tb.PDSI$V[i-1])))) 
    } else if (j==7){ # ... Ze 
     Tb.PDSI[[i,j]]<-ifelse(Tb.PDSI$O[i-1]==0,0, 
          ifelse(Tb.PDSI$X3[i-1]<0 && Tb.PDSI$Z[i]>0,-par.b_n*Tb.PDSI$X3[i-1]-mbO, 
            ifelse(Tb.PDSI$X3[i-1]>0 && Tb.PDSI$Z[i]<0,-par.b_n*Tb.PDSI$X3[i-1]+mbO, 
              ifelse(Tb.PDSI$Ze[i-1]==0,0, 
                ifelse(Tb.PDSI$O[i-1]==1,-par.b_n*Tb.PDSI$X3[i-1]+mbO,-par.b_n*Tb.PDSI$X3[i-1]-mbO))))) 
    } else if (j==8){ # ... Q 
     Tb.PDSI[[i,j]]<-ifelse(Tb.PDSI$O[i-1]==0,0, 
          ifelse(Tb.PDSI$P[i-1]==100,Tb.PDSI$Ze[i],Tb.PDSI$Ze[i]+Tb.PDSI$V[i-1])) 
    } else if (j==9){ # ... P 
     Tb.PDSI[[i,j]]<-ifelse(Tb.PDSI$Q[i]==0,0, 
          ifelse(100*Tb.PDSI$V[i]/Tb.PDSI$Q[i]>100,100,100*Tb.PDSI$V[i]/Tb.PDSI$Q[i])) 
    } else if (j==10){ #########  WET DATA ######### 
            # ... Ms.1 
     Tb.PDSI[[i,j]]<-ifelse(Tb.PDSI$O[i-1]==1 && Tb.PDSI$P[i]<100,0,-Tb.PDSI$X1[i-1]*mmb) 
    } else if (j==11){    # ... AX1 
     Tb.PDSI[[i,j]]<-ifelse(Tb.PDSI$O[i-1]==1 && Tb.PDSI$P[i]<100,0, 
          ifelse(Tb.PDSI$X1[i-1]==0,0,Tb.PDSI$Z3[i]+Tb.PDSI$Ms.X1[i])) 
    } else if (j==12){    # ... X1 
     Tb.PDSI[[i,j]]<-ifelse(Tb.PDSI$O[i-1]==1,0, 
          ifelse(ifelse(Tb.PDSI$X1[i-1]==0, 
              ifelse(Tb.PDSI$Z3[i]>0,Tb.PDSI$Z3[i],0),Tb.PDSI$AX.1[i]+Tb.PDSI$X1[i-1])>0, 
            ifelse(Tb.PDSI$X1[i-1]==0, 
              ifelse(Tb.PDSI$Z3[i]>0,Tb.PDSI$Z3[i],0),Tb.PDSI$AX.1[i]+Tb.PDSI$X1[i-1]),0)) 
    } else if (j==13){########  drought data  ######### 
           # ... Ms.2 
     Tb.PDSI[[i,j]]<-ifelse(Tb.PDSI$O[i-1]== -1 && Tb.PDSI$P[i]<100,0,Tb.PDSI$X2[i-1]*-mmb) 
    } else if (j==14){   # ... AX2 
     Tb.PDSI[[i,j]]<-ifelse(Tb.PDSI$O[i-1]==-1 && Tb.PDSI$P[i]<100,0,ifelse(Tb.PDSI$X2[i-1]==0,0,Tb.PDSI$Z3[i]+Tb.PDSI$Ms.X2[i])) 
    } else if (j==15){   # ... X2 
     Tb.PDSI[[i,j]]<-ifelse(Tb.PDSI$O[i-1]== -1,0, 
          ifelse(
           ifelse(Tb.PDSI$X2[i-1]==0, 
             ifelse(Tb.PDSI$Z3[i]<0, Tb.PDSI$Z3[i],0),Tb.PDSI$AX.2[i]+Tb.PDSI$X2[i-1])<0, 
           ifelse(Tb.PDSI$X2[i-1]==0, 
             ifelse(Tb.PDSI$Z3[i]<0,Tb.PDSI$Z3[i],0),Tb.PDSI$AX.2[i]+Tb.PDSI$X2[i-1]),0)) 
    } else if (j==16){ ########  Established event ########### 
             # ... Ms.3 
     Tb.PDSI[[i,j]]<-ifelse(Tb.PDSI$O[i-1]==0 || Tb.PDSI$P[i]==100,0,-mmb*Tb.PDSI$X3[i-1]) 
    } else if (j==17){     # ... AX3 
     Tb.PDSI[[i,j]]<-ifelse(Tb.PDSI$O[i-1]==0 || Tb.PDSI$P[i]==100,0,Tb.PDSI$Z3[i]+Tb.PDSI$Ms.X3[i]) 
    } else if (j==18){     # ... X3 
     Tb.PDSI[[i,j]]<-ifelse(Tb.PDSI$O[i-1]==0 || Tb.PDSI$P[i]==100, 
          ifelse(Tb.PDSI$X1[i]>1,Tb.PDSI$X1[i], 
            ifelse(Tb.PDSI$X2[i]<-1,Tb.PDSI$X2[i],0)),Tb.PDSI$X3[i-1]+Tb.PDSI$AX.3[i]) 
    } else if (j==19){ ######## FINAL ############ 
           # ... O 
     Tb.PDSI[[i,j]]<-ifelse(Tb.PDSI$O[i-1]==0, 
          ifelse(Tb.PDSI$X1[i]>1,1, 
            ifelse(Tb.PDSI$X2[i]<-1,-1,0)),ifelse(Tb.PDSI$P[i]==100,ifelse(Tb.PDSI$X1[i]>1,1,ifelse(Tb.PDSI$X2[i]<-1,-1,0)),Tb.PDSI$O[i-1])) 
    }else if (j==20){   # ... Final Index 
    Tb.PDSI[[i,j]]<-ifelse(Tb.PDSI$O[i]==0, ifelse(Tb.PDSI$X1[i]+Tb.PDSI$X2[i]>0,Tb.PDSI$X1[i],Tb.PDSI$X2[i]),Tb.PDSI$X3[i]) 
    }else if (j==21){   # ... Duracion Meses 
    Tb.PDSI[[i,j]]<-ifelse(Tb.PDSI$O[i]==Tb.PDSI$O[i-1] && abs(Tb.PDSI$O[i])>0, 
          ifelse(Tb.PDSI$O[i-1]==1,Tb.PDSI$Dur.Mes[i-1]+1,Tb.PDSI$Dur.Mes[i-1]-1),Tb.PDSI$O[i]) 
    } 
    } 
} 

這是我的結果R

Anho Mes   Z Uw.Z  Ud.Z V  Ze   Q P  Ms.X1  AX.1  X1 Ms.X2 AX.2 X2  Ms.X3  AX.3 
1: 1999 12 0.00000 0 0.000000 0 0.00000 0.00000 0 0.00000000 0.0000000 0.0000000  0 0 0 0.00000000 0.0000000 
2: 2000 1 22.82134 0 0.000000 0 0.00000 0.00000 0 0.00000000 0.0000000 0.6334310  0 0 1 0.00000000 0.0000000 
3: 2000 2 39.40637 0 0.000000 0 0.00000 0.00000 0 -0.04303857 1.0507278 1.6841588  0 0 0 -0.06794515 1.0258212 
4: 2000 3 -16.39223 0 -17.616196 0 -50.01341 -50.01341 0 -0.11443043 -0.5694145 1.1147443  0 0 0 -0.13764473 -0.5926288 
5: 2000 4 -8.35648 0 -9.580449 0 -30.11282 -30.11282 0 -0.07574147 -0.3076846 0.8070597  0 0 0 -0.09737848 -0.3293216 
6: 2000 5 -10.62170 0 -11.845673 0 -19.05413 -19.05413 0 -0.05483579 -0.3496527 0.4574070  0 0 0 -0.07500267 -0.3698196 
      X3 O IndeX.f Dur.Mes   Z3 
1: 0.0000000 0 0.0000000  0 0.0000000 
2: 1.0000000 -1 1.0000000  -1 0.6334310 
3: 2.0258212 -1 2.0258212  -2 1.0937664 
4: 1.4331925 -1 1.4331925  -3 -0.4549841 
5: 1.1038708 -1 1.1038708  -4 -0.2319432 
6: 0.7340512 -1 0.7340512  -5 -0.2948169 

這是正確的我在Excel中得到了和手工驗證

Anho Mes   Z Uw.z  Ud.z   V   Ze   Q   P  Ms.X1  AX1 
1: 1999 12 0.000000 0 0.000000 0.00000 0.000000 0.00000 0.00000 0.00000000 0.000000 
2: 2000 1 22.818904 0 0.000000 0.00000 0.000000 0.00000 0.00000 0.00000000 0.000000 
3: 2000 2 39.402163 0 0.000000 0.00000 0.000000 0.00000 0.00000 -0.04303857 1.050728 
4: 2000 3 -16.390479 0 -17.614317 -17.61432 -38.536210 -38.53621 45.70848 0.00000000 0.000000 
5: 2000 4 -8.355589 0 -9.579427 -27.19374 -19.417198 -37.03152 73.43406 0.00000000 0.000000 
6: 2000 5 -10.620571 0 -11.844410 -39.03815 -9.086188 -36.27993 100.00000 0.00000000 0.000000 
     X1  Ms.X2  AX2   X2  Ms.X3  AX3  X3 O IndeX.f Dur.Mes 
1: 0.000000 0.00000000 0.0000000 0.0000000 0.00000000 0.0000000 0.0000000 0 0.0000000  0 
2: 0.633431 0.00000000 0.0000000 0.0000000 0.00000000 0.0000000 0.0000000 0 0.6334310  0 
3: 1.684159 0.00000000 0.0000000 0.0000000 0.00000000 0.0000000 1.6841588 1 1.6841588  1 
4: 0.000000 0.00000000 0.0000000 -0.4549840 -0.11443043 -0.5694145 1.1147443 1 1.1147443  2 
5: 0.000000 0.03091396 -0.2010292 -0.6560133 -0.07574148 -0.3076846 0.8070597 1 0.8070597  3 
6: 0.000000 0.04457292 -0.2502440 -0.9062573 0.00000000 0.0000000 0.0000000 0 -0.9062573  0 
      Z3 
1: 0.0000000 
2: 0.6334310 
3: 1.0937664 
4: -0.4549840 
5: -0.2319432 
6: -0.2948169 

我不知道爲什麼R給我的那些不正確值。

+1

你應該做一個小例子的工作。我懷疑你需要你的完整腳本或這麼多列來說明這個問題。請參閱[mcve]和https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example/28481250#28481250 – Frank

+1

您已經給我們沒有理由認爲您的代碼在Excel應該是一個「黃金標準」。你明顯在兩個不同的平臺上做了一些不同的事情,我們不可能知道哪一個是正確的(因爲你沒有用自然語言解釋算法可能是什麼,而且我們也沒有Excel工作表。 –

回答

0

謝謝你的批評我終於解決了這個問題對我自己的,這是因爲空間ifelse(Tb.PDSI$X2[i]<-1,當正確的是ifelse(Tb.PDSI$X2[i]< (-1)