knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
We illustrate the process of fitting a probabilty distribution to an expert's judgements, and obtaining feedback.
Various elicitation methods are suggested in SHELF, but they all involve either asking for probability or quantile judgements: the expert is either asked questions of the form, 'What is your probability $P(X\le10)$?' or, 'Provide a value $x$ such that you think $P(X\le x)=0.5.$'' Either way, we think of the expert as providing points on his/her cumulative distribution function (CDF).
As an example, suppose the uncertain quantity is the percentage $X$ of voters who will vote for a particular candidate in an upcoming leadership election.
Suppose we have two experts, who first consider their judgements independently. Expert A states $$P(X\le25)=0.25, P(X\le30)=0.5, P(X\le35)=0.75,$$ and Expert B states $$P(X\le30)=0.25, P(X\le35)=0.5, P(X<50)=0.75.$$
The values of the quantiles are arranged in a matrix, with one column per expert
v <- matrix(c(25, 30, 35, 30, 35, 50), nrow = 3, ncol = 2)
v
The probabilities should be arranged in matrix if they are different for the two experts (for example, if one expert had specified tertiles instead of lower and upper quartiles), otherwise they are stored in a vector
p <- c(0.25, 0.5, 0.75)
We now fit probability distributions to these judgements with the fitdist
command. Without specifying any limits (lower
and upper
), only normal and Student-$t$ distributions (with 3 degrees of freedom used as a default) will be fitted. Including a lower
limit will result in log-normal, log Student-$t$ and gamma distributions also being fitted, and including both a lower
and upper
limit will result in a beta distribution fit, scaled to the interval [lower
, upper
].
library(SHELF) myfit <- fitdist(vals = v, probs = p, lower = 0, upper = 100)
The object myfit
includes parameters of the fitted distributions for each expert, sums of squared errors (elicited - fitted probabilities) for each distribution and expert, and the original elicited judgements:
names(myfit)
For example, the parameters of the fitted beta distributions are
myfit$Beta
Inspecting myfit$ssq
(sum of squared differences between elicited and fitted probabilities), we see that the normal distribution fits best for Expert 1, and the log Student-$t_3$ distribution fits best for Expert 2 (although the fit would probably be judged adequate for any the distributions, in this case). These best-fitting distributions are used as defaults in the plotfit
function, and the feedback
function when used with multiple experts
myfit$ssq
We plot the fitted distributions, including a linear pool.
plotfit(myfit, lp = TRUE)
Now we elicit a single a single 'consensus' distribution from the two experts. Suppose they agree $$P(X\le 25)=0.25, P(X\le30)=0.5, P(X\le40)=0.75.$$ We fit a single set of distributions to these judgements.
v <-c(25, 30, 40) p <-c(0.25, 0.5, 0.75) consensus <- fitdist(vals = v, probs = p, lower = 0, upper = 100)
The best fitting distribution is the log Student-$t_3$ distribution, but the (scaled) beta distribution would likely be adequate at this stage, and we may choose to use the beta instead, given it has the appropriate bounds.
We plot the fitted density, and mark the fitted 5th and 95th percentiles
plotfit(consensus, ql = 0.05, qu = 0.95, d = "beta")
We can obtain values of fitted percentiles and probabilities with the feedback
function. This will show fitted values for all the fitted distributions
feedback(consensus, quantiles = c(0.05, 0.95))
For the fitted beta distribution, we have $P(X<15.2)=0.05$, and we can also compare, for example, the fitted probability $P(X<25)=0.283$ with elicited probability $P(X<25)=0.25$.
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.