How do I calculate the probability for a given quantile in R?

ecdf returns a function: you need to apply it.

f <- ecdf(x)
f( quantile(x,.91) )
# Equivalently:
ecdf(x)( quantile(x,.91) )

Just for convenience, this function helps:

quantInv <- function(distr, value) ecdf(distr)(value)
set.seed(1)
x <- rnorm(1000, mean=4, sd=2)
quantInv(x, c(4, 5, 6.705755))
[1] 0.518 0.685 0.904

You more or less have the answer yourself. When you want to write

backwards_quantile(x, 5)

just write

ecdf(x)(5)

This corresponds to the inverse of quantile() with type=1. However, if you want other types (I favour the NIST standard, corresponding to Excel's Percentile.exc, which is type=6), you have more work to do.

In these latter cases, consider which use you are going to put it to. If all you want is to plot it, for instance, then consider

yVals<-seq(0,1,0.01)
plot(quantile(x,yVals,type=6))

But if you want the inverse for a single value, like 5, then you need to write a solving function to find the P that makes

quantile(x,P,type=6) = 5

For instance this, which uses binary search between the extreme values of x:

inverse_quantile<-function(x,y,d=0.01,type=1) {
  A<-min(x)
  B<-max(x)
  k<-(log((B-A)/d)/log(2))+1
  P=0.5
  for (i in 1:k) {
    P=P+ifelse((quantile(x,P,type=type)<y),2^{-i-1},-2^{-i-1})
  }
  P
}

So if you wanted the type 4 quantile of your set x for the number 5, with precision 0.00001, then you would write

inverse_quantile<-function(x,5,d=0.00001,type=4)