I am not sure if I understood you well. However I think you can construct the empirical probability distribution of the error (using test set), if necessary using bootstrapping. Then you can estimate probability distribution of the error using MLE, use non-parametric estimation, etc.
If you have a distribution, then you can easily estimate the probability of exceeding some given value of error.
Please notice, that you require a lot of observations to make a good estimation of the distribution. In case of this dataset (iris) this approach will not provide very robust outcomes.
If this approach fails, you can always use Chebyshev's inequality
Below you have a R code for a toy example:
df <- iris
test.df <- iris #in your case it will be a real test set
rf.model <- randomForest(Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width ,data=df)
check.df <- data.frame(predicted = predict(rf.model,newdata=test.df), observed = test.df$Sepal.Length)
#uncomment line below to display a scatter plot
#ggplot(data = check.df,aes(x=observed,y=predicted)) + geom_point() + geom_smooth(method="lm")
error.df <- data.frame(error = check.df$observed - check.df$predicted)
#uncomment line below to display a histogram
#ggplot(data = error.df,aes(x=error,y=..density..)) + geom_histogram() + geom_density() # plot the histogram and density
cdf <- ecdf(x=error.df$error) #empirical CDF function
#probability of error exceeding 0.2 on both sides, so including also error exceeding -0.2, P(error<-0.2 or error>0.2)
P = cdf(-0.2) + (1-cdf(0.2))
Of course you can use any measure you like for the error.
Also prediction intervals may be interesting for you.