return(list(

rps = rps,

meanrps = mean(rps),

sumrps = sum(rps)

))

I havent tested your code, but you seem to make an error where you return the values for the last game. This is because you in each iteration overwrite the previous values of rankscore. The cumulative variable should also not be returned, since this too is overwritten in each for each game.

]]>I want to use next (very unrealistic) example: Let’s assume that the attack/defense parameters distributions are equal for both leagues and the winner will promote to PL. In that case the model will assume that the quality of the promoted team is equal to last year’s champion of the PL. Can you eleborate on why the model accuracy improves when we add games from another sub-league (where the level is significant lower)? ]]>

When I use this function, I provide 2 inputs: a predictions outcome matrix with H, D, A probabilities for all matches and a vector with values 1, 2 or 3 depending on the actual outcome of the match.

However, when the rps is returned, I get a value for each match. Is that what it should be? Or should it be one numerical value? Or is it a problem with my inputs to the function?

I reedited the code as follows: Firstly, I returned a list just to see what the function is doing. One thing I notice is that returning “obsvec” only returns one line: [1] 0 0 1.

Then I added one line:

rps[rr] <- (1/(ncat-1))*cumulative

rankscore <- (1/(ncat-1))*cumulative

the rankscore returns one numerical value. In the case of xi = 0.0018 I get rankscore = 0.5849709, using data from Premiership & Championship 2015-16 and 2016-17.

Also, I got an interesting return value for "cumulative" of 1.169942. Not quite sure what that means, but it doesn't sound right.

Below the code:

rankProbScore <- function(predictions, observed){

ncat <- ncol(predictions)

npred <- nrow(predictions)

rps <- numeric(npred)

for (rr in 1:npred){

obsvec <- rep(0, ncat)

obsvec[observed[rr]] <- 1

cumulative <- 0

for (i in 1:ncat){

cumulative <- cumulative + (sum(predictions[rr,1:i]) – sum(obsvec[1:i]))^2

}

rps[rr] <- (1/(ncat-1))*cumulative

rankscore <- (1/(ncat-1))*cumulative

}

return(list(

ncat,

npred,

obsvec,

cumulative,

rps,

rankscore

))

}

I finished the two stage approach and will convert the real DC model to Matlab in a while. ]]>