Probability of a value given a dataset I am trying to figure out a way to determine how extreme a particular value is from single variable, when the distribution of the variable is unknown. 
I have tried doing this with a kind of boot-strap approach using the observed data (which presumably has a proper name but I can't figure out what it is), which I coded in R and is below. But there feels something wrong with this approch and I am hoping someone here can point it out. I also think there is probably a good bayesian approach to this problem but I am not sure what it is.
What I would like to know is whether this approach is on the right track, and what kind of terms I should be searching to learn more about this. I think my question is quite similar to this one but this seems to correcting to a normal distribution, which I think I want to avoid. 
Thanks,
Sam
d = rbeta(100, 1, 10) # pretend you don't know this. 
#plot(density(d))
x = 0.01

score = c()
for(i in 1:1000){
  value = sample(d, 1)
  score[i] = value >= x
}
mean(score)

 A: I'm no expert but I'll try to give a shot at your question. If I understand it
correctly, you have a sample of datapoints d but you don't know the
parameters of the distribution that generated it. You want to know how extreme
the value x is in that sample. So you have:
set.seed(1234)
d <- rbeta(100, 1, 10) # You don't know the parameters
x <- 0.01

To know how extreme x is you could simply check where x ranks in d and turn
the rank into a percentage:
sum(d <= x) / length(d) # x is in the lowest 9% of d

You could then bootstrap d and get some confidence interval for the ranking
of x:
score <- rep(NA, 10000)
for(i in 1:length(score)){
  d1 <- sample(d, size= length(d), replace= TRUE)
  score[i] <- sum(d1 <= x) / length(d1)
}
quantile(score, p= c(0.05, 0.95)) 
#   5%  95% 
# 0.05 0.14   # x is between the 5% to 14% lowest of d


Bayesian approach
A Bayesian approach could go along these lines. x is anywhere between 0% and
100% of d (i.e. x can be anywhere between being the lowest or largest value
of d.) Since this range is between 0 and 1, we use a Beta prior to
describe our guess for where x could be. For example, choose:
alpha <- 2
beta <- 3
prior <- dbeta(seq(0, 1, by= 0.01), alpha, beta) 
plot(seq(0, 1, by= 0.01), prior, type= 'l', xlim= c(0, 1))

Then we consider the values of d $<=$ x as "number of successes" and those
d $>$ x as "number of failures". I.e. our data comes from a binomial
distribution. In this case updating the prior becomes very simple since the
posterior is just $Beta(alpha + successes, beta + failures)$ (reference
needed, please Google it up). So in our case we have:
succ <- sum(d <= x)
fail <- sum(d > x)

posterior <- dbeta(seq(0, 1, by= 0.01), alpha + succ, beta + fail)

We can use the posterior distribution and the updated parameters to plot the
results:
plot(seq(0, 1, by= 0.01), posterior, type= 'l', xlim= c(0, 1), col= 'red', ylab= 'Density', xlab= 'x')
points(seq(0, 1, by= 0.01), prior, type= 'l', col= 'blue')
hist(d, freq= FALSE, add= TRUE, breaks= 20)
legend(x= 'topright', legend= c('prior', 'posterior', 'data'), col= c('blue', 'red', 'black'), lty= 'solid')

You could sample from the posterior distribution to get an estimate of how extreme x is: 
mean(rbeta(10000, alpha + succ, beta + fail)) # ~ 10%


Hope this is correct and makes sense - I'd like to hear other views on the problem...
A: For a set of values that are NON-INTEGER, one can only compute interval probabilities. If your data is cross-sectional (i.e. no autocorrelative structure) then simply create a histogram and compute the cdf to get your interval probabilities. 
If your data is either chronological or spatial , one can use the distribution of the model residuals as the basis for a monte-carlo to bootstrap the cdf .
If the data set is all integers then count the number of each kind of observation and estimate the percentage of the total to get probabilities.
