# Rao-Stirling diversity index

I need help to calculate the Rao-Stirling diversity index. I tried it several times but cannot achieve equal results to R packages yet (e.g. diverse).

I used the Rao-Stirling diversity index as suggested by several references $$\Delta=\sum_{ij} p_j p_i d_{ij}$$

including the reciprocal cosine similarity for computing the disparity $$d_{ij}=(1-'cosine')$$

When I would have the following simple matrix, please could one show me how to do it?

$$\begin{bmatrix} 0 & 10 & 20 \\ 10 & 0 & 30 \\ 20 & 30 & 0 \end{bmatrix}$$

Edit 02/01/20

delta <- function(matrix){
proportion.vector <- prop.vector(matrix)
distance.matrix <- dist.sim.matrix(matrix)
result <- matrix(rep(NA,nrow(matrix)*ncol(matrix)),nrow=nrow(matrix))
for (i in 1:nrow(matrix)){
for (j in 1:ncol(matrix)){

result[i,j] <- proportion.vector[i]*proportion.vector[j]*distance.matrix[i,j] } }

return(colSums(result))
}


Oh sorry, I missed the point. Here is the code that I have used.

## Computing Rao-Stirling Diversity
#install.packages('diverse')
library(diverse)

rows <- c(1,1,1,2,2,2,3,3,3)
cols <- c(1,2,3,1,2,3,1,2,3)
values <- c(0,10,20,10,0,30,20,30,0)
my.dataframe <- data.frame(rows,cols,values)

RAO = diversity(data=my.dataframe, type="rao-stirling", method="cosine")

• The formula requires the frequencies $p_i$ as well as the similarity matrix $d_{ij}.$ How, then, do you propose to implement it when given only a "simple matrix"?? What does this input matrix represent?
– whuber
Jan 3 '20 at 15:14

I just had a cursory look at the literature. Here comes some code I just wrote:

    # computes cosine similarity between vec1 and vec2
cos.sim <- function(vec1,vec2){
return(sum(vec1*vec2)/sqrt(sum(vec1*vec1)*sum(vec2*vec2)))
}

# returns D_ij matrix
dist.sim.matrix<- function(M){
result <- matrix(rep(NA,nrow(M)^2),nrow=nrow(M))
for (i in 1:nrow(M)){
for (j in 1:nrow(M)){
result[i,j] <- cos.sim(M[,i],M[,j])
}
}
return(1-result)
}

# returns proportion vector
prop.vector <- function(matrix){
result <- numeric(ncol(matrix))
total.sum <- sum(matrix)
for (i in 1:ncol(matrix)){
result[i] <- sum(matrix[,i])/total.sum
}
return(result)
}
# returns delta
delta <- function(matrix){
proportion.vector <- prop.vector(matrix)
distance.matrix <- dist.sim.matrix(matrix)
result <- matrix(rep(NA,nrow(matrix)^2),nrow=nrow(matrix))
for (i in 1:nrow(matrix)){
for (j in 1:nrow(matrix)){
result[i,j] <-            proportion.vector[i]*proportion.vector[j]*distance.matrix[i,j]
}
}
return(sum(result))
}

your.matrix <- matrix(c(0,10,20,10,0,30,20,30,0),nrow=3,byrow=TRUE)