How to interpret/report regression scoring in the form of probabilities? Around a hundred features are measured from a product.
I've successfully modelled many of them based on a features subset thanks to historical data and random forest regression with very high R2 scores (0.95).
How and what evaluation metric should I use in order to be able to make a statement like:
'The physical measurement of this feature can be stopped, replaced by the model's prediction because the chance of the error being greater than XX is only YY%'
Or something similar that could be convincing enough for a quality responsible.
Thanks in advance
 A: 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:
library(randomForest)
library(ggplot2)
data(iris)

set.seed(2020)

df <- iris
test.df <- iris #in your case it will be a real test set
head(df)

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)
head(check.df)
#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)) 
P

Of course you can use any measure you like for the error.
Also prediction intervals may be interesting for you.
