-3

I'm working on the calculation of PDSI with a simple script build by myself, but when i ran the script i get a different value in the outputs in the DT. I did all this calculation manually in and sheet of paper and also in excel and verified both and i get the correct values.

Here is an extract of the script and i didnt have any error so i don´t know where is my error.

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])
   }
  }
}

This is my result in 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

This is the correct one i got in excel and verified manually

    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

I don´t know why R gave me those incorrect values.

Guisseppe
  • 3
  • 4
  • 1
    You should work on making a minimal example. I doubt you need your full script or this many columns to illustrate the issue. See [mcve] and https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example/28481250#28481250 – Frank Jun 15 '17 at 23:23
  • 1
    You've given us no reason to think that your coding within Excel should be a "gold standard". You clearly have done something that is different on two different platforms and there is no way we can possbly know which one is correct (since you have not explained in natural language what the algorithm might be and we also don't have the Excel worksheet. – IRTFM Jun 16 '17 at 00:13

1 Answers1

0

Thanks for the critics finally I solved the problem for my own, it was because of the spaces ifelse(Tb.PDSI$X2[i]<-1 , when the correct was ifelse(Tb.PDSI$X2[i]< (-1)

Guisseppe
  • 3
  • 4