I want to do a QQ plot, where I compare a sample to multiple standard distributions, scaled, so that the points represented by the correct distribution are on the black line.
I expect to do basically what https://www.youtube.com/watch?v=okjYjClSjOg describes.
What I aim to do
- Sort the sample values.
- Consider each value to be the sample quantile of a distribution.
- From this distribution's cumulative probability function, take as many quantiles as one has sample values. Takes those expected quantiles in such a way, that the integral of the probability density function between two adjacent expected quantiles is always the same.
- Sort the expected quantiles and align them with the sample quantiles, so that one gets points.
- If the sample distribution is the same as the distribution, the x-coordinate and y-coordinate of a point should be approximately the same: the points are close to the line passing through the points (0,0) and (1,1).
I do this not only for one distribution, but for multiple distributions and put all the resulting points into the plot.
I implemented this, but the result is not what I thought it would be (as one can see below). Where did I go wrong in terms of my understanding?
- I sort my values and for each value, I expect the cumulative
probabilityto increase (thus the cumulative sum
cumsum) the same amount (namely
n()is the number of values that I have).
- For each of those probabilities I calculate the expected quantile for different
distributions. Below, I name them
At this point I am done with the core implementation, the rest is clean-up:
- I bring the data into a long format.
- I remove points with too large quantiles or sample values.
- For each distribution, I standardize the coordinates, so that the are easier to see.
- I plot the quantiles at the x-axis and the original values at the y-axis.
- I plot the points.
- I expect the resulting dots for the right distribution to be on the line passing (0,0) and (1,1).
I was sure I got it right and expressed this in code:
data.frame(eval_val = rlnorm(10000, 8, 6)) %>% arrange(eval_val) %>% mutate(prob_mass = 1/n(), probability = cumsum(prob_mass), percentile_norm = qnorm(probability), percentile_beta= qbeta(probability, shape1 = 1, shape2 = 2), percentile_lognorm = qlnorm(probability), ) %>% gather(distribution, percentile, starts_with("percentile_")) %>% filter(is.finite(percentile), eval_val < 1000) %>% group_by(distribution) %>% mutate(percentile = (percentile - mean(percentile)) / sd(percentile), eval_val = (eval_val - mean(eval_val)) / sd(eval_val)) %>% ggplot(aes(x=percentile, y=eval_val, color=distribution)) + geom_point() + geom_abline(slope = 1, intercept = 0)
but as one can see, the log-normal distribution is not on the line and neither is any of the others distributions.
qlnormwith the default parameters (0 and 1) but
rlnormwith 8 and 6. (If you simplify the code to one comparator distribution at a time and it still doesn't work, you may have a statistical issue. If you simplify it and it does work, you've got a coding issue. ) $\endgroup$