

zrpss<-function (ydat) 
{
	obs <- ydat[,1]
	pred <- ydat[,2:4]
	baseline <- c(1/3,1/3,1/3)
    nr <- nrow(pred)
    id <- is.finite(obs) & is.finite(apply(pred, 1, sum))
    obs <- obs[id]
    pred <- matrix(pred[id, ], nrow = nr)
    OBS <- matrix(0, nrow = length(obs), ncol = ncol(pred))
    for (i in 1:nrow(OBS)) {
        OBS[i, obs[i]] <- 1
    }
    OBS2 <- OBS
    for (i in 1:ncol(OBS)) {
        OBS2[, i] <- apply(matrix(OBS[, 1:i], nrow = nr), 1, sum)
    }
    PRED <- OBS
    for (i in 1:ncol(pred)) {
        PRED[, i] <- apply(matrix(pred[, 1:i], nrow = nr), 1, sum)
    }
    RPS <- mean(apply((PRED - OBS2)^2, 1, sum))  #/(ncol(pred) -  1) # As in Wilks ZTS
    if (is.null(baseline)) {
        xxx <- apply(OBS, 2, sum)/nrow(OBS)
        pred.climo <- matrix(xxx, nrow = nrow(OBS), ncol = ncol(OBS), byrow = TRUE)
    }
    else {
        pred.climo <- matrix(baseline, nrow = nrow(OBS), ncol = ncol(OBS), byrow = TRUE)
    }
    PRED.climo <- OBS
    for (i in 1:ncol(pred.climo)) {
        PRED.climo[, i] <- apply(matrix(pred.climo[, 1:i], nrow = nr), 1, sum)
    }
    RPS.climo <- mean(apply((PRED.climo - OBS2)^2, 1, sum))   #/(ncol(PRED.climo) -  1)  #for WILS ZTS
    RPSS <- 1 - RPS/RPS.climo
    RPSS
    #return(list(rps = RPS, rpss = RPSS, rps.clim = RPS.climo))
}


