restart: The problem presented here also exists with Maple 2018 and Maple 2019 interface(version) SWJwU3RhbmRhcmR+V29ya3NoZWV0fkludGVyZmFjZSx+TWFwbGV+MjAxNS4yLH5NYWN+T1N+WCx+RGVjZW1iZXJ+MjF+MjAxNX5CdWlsZH5JRH4xMDk3ODk1RzYi with(Statistics): X := RandomVariable(Normal(0, 1)): Generate a sample of X of size 10^6 Here is the R code for those who would like to verify the results and the performances N <- 10^6 R <- 20 Q <- c(0:5) M <- matrix(nrow=R, ncol=length(Q)) for (r in c(1:20)) { S <- rnorm(N, mean=0, sd=1) for (q in Q) { M[r, q+1] <-length(S[S>q]) } } Expected.Number.of.Outliers <- floor(N - pnorm(Q, mean=0, sd=1) * N) Observed.Number.of.Outliers <- round(colMeans(M)) Expected.Number.of.Outliers Observed.Number.of.Outliers Absolute.Differences <- M - kronecker(t(Expected.Number.of.Outliers), rep(1, R)) boxplot(differences) Remark the loop below takes about 30 seconds to run (to be compared to the 20 minutes it would take am I used the build-in procedure Select, and to the 3 seconds R demands). N := 10^6: R := 20: Q := [\$0..5]: nQ := numelems(Q): M := Matrix(R, nQ, 0): for r from 1 to R do S := Sample(X, N): for q in Q do M[r, q+1] := add(Heaviside~(S -~ q)): end do: end do: Expected_Number_Of_Outliers := Vector[row](nQ, i -> Probability(X > Q[i], numeric)*N); LSZJJ1ZlY3Rvckc2JCUqcHJvdGVjdGVkR0koX3N5c2xpYkc2IjYjSSRyb3dHRig2Iy9JJCVpZEdGKCI1bWxoIj15U3VZJT0= Observed_Number_Of_Outliers := Mean(M); LSZJJ1ZlY3Rvckc2JCUqcHJvdGVjdGVkR0koX3N5c2xpYkc2IjYjSSRyb3dHRig2Iy9JJCVpZEdGKCI1STMyJj15U3VZJT0= While the values of the Observed_Number_of_Outliers seem to be reasonably correct for q =0, 1, 2 and 3 there are highly suspicious for q = 4 and likely for q = 5 (but N is too small to be categorical) So let's examine more closely rhe results for q = 4 `q=4`:= convert(M[..,5], list); min(`q=4`) NzZdNDA0ODgwMDAwMDAwMDAwMF00MDRGODAwMDAwMDAwMDAwXTQwNDUwMDAwMDAwMDAwMDBdNDA0NzgwMDAwMDAwMDAwMF00MDQ2MDAwMDAwMDAwMDAwXTQwNTEwMDAwMDAwMDAwMDBdNDA0NzAwMDAwMDAwMDAwMF00MDQ5ODAwMDAwMDAwMDAwXTQwNEQwMDAwMDAwMDAwMDBdNDA0QjgwMDAwMDAwMDAwMF00MDQ5MDAwMDAwMDAwMDAwRipGKF00MDQ2ODAwMDAwMDAwMDAwRitGJ100MDREODAwMDAwMDAwMDAwRiNGJl00MDQzMDAwMDAwMDAwMDAw XTQwNDMwMDAwMDAwMDAwMDA= This means we obtained 20 estimations out of 20 of the probability that X exceeds q=4. There are only about on chance out of 1 million to get such a result (roughly (1/2)^20 ). Here is the result R returned for the same quantile q=4 M[,5]  24 37 23 38 30 43 23 30 25 39 25 29 29 36 39 41 35 37 36 34 (9 values out of 20 are less than the expected value of 31.67) A very simple test based on the Binomial distribution will prove, without any doubt, that the distribution of the 20 numbers of outliers for q = 4 implies the sampling algorithm is of is of very poor quality. Is it possible to obtain correct results with Maple ? (simulations done below for the case q=4 only) N := 10^6: R := 20: Menv := Vector[column](R, 0): for r from 1 to R do S := Sample(X, N, method=envelope): Menv[r] := add(Heaviside~(S -~ 4)): end do: `q=4`:= convert(Menv, list); Mean(Menv); Variance(Menv); NzZdNDAzQTAwMDAwMDAwMDAwMEYjXTQwNDE4MDAwMDAwMDAwMDBdNDA0MTAwMDAwMDAwMDAwMF00MDNFMDAwMDAwMDAwMDAwRiZdNDAzQzAwMDAwMDAwMDAwMF00MDM4MDAwMDAwMDAwMDAwXTQwM0IwMDAwMDAwMDAwMDBdNDAzRjAwMDAwMDAwMDAwMEYnXTQwMzUwMDAwMDAwMDAwMDBGJV00MDQwODAwMDAwMDAwMDAwXTQwM0QwMDAwMDAwMDAwMDBdNDA0MzAwMDAwMDAwMDAwMEYuRixdNDA0MzgwMDAwMDAwMDAwMF00MDM5MDAwMDAwMDAwMDAw XTQwM0U3MzMzMzMzMzMzMzM= XTQwMzhFNDYxMDJCMURBNDY= Conclusion, the "default" sampling strategy suffers strong flaw. Remark: it's well known that one of the best method to sample normal random variables is Box-Mueller's: why it has not been programmed by default ? TTdSMApJQVJUQUJMRV9TQVZFLzE4NDQ2NzQ0MDc4MTgxNjE2NTY2WColKWFueXRoaW5nRzYiRiVbZ2whJCUhISEiJyInXTQxMUU4NDgwMDAwMDAwMDBdNDEwMzVERkEwODBEMzczOV00MEQ2Mzc4ODcxRDZDNjlFXTQwOTUxNzk3OTU5QTg4NDZdNDAzRkFCRDY4MTM4RjUxOF0zRkQyNTg3RkQ1QjA4QzgwRiU=TTdSMApJQVJUQUJMRV9TQVZFLzE4NDQ2NzQ0MDc4MTg1MDcwODMwWColKWFueXRoaW5nRzYiRiVbZ2wnJCUhISEiJyInNDExRTg5MEQ5OTk5OTk5QTQxMDM2MUYyNjY2NjY2NjY0MEQ2M0E1NjY2NjY2NjY2NDA5NTE2MDAwMDAwMDAwMDQwNDlDQ0NDQ0NDQ0NDQ0QzRkZFNjY2NjY2NjY2NjY2RiU=