[Feature request] Conformal Predictive Distributions about conformalprediction.jl HOT 5 OPEN

using MLJ, EvoTrees, ConformalPrediction, Plots, Random, MLJLinearModels, Tables;
########################################
n= 100_000; p=70; σ=100.10;
X = [ones(n) randn(MersenneTwister(49), n, p-1)]
θ = randn(MersenneTwister(49), p)
CEF   = X*θ 
Noise = σ*randn(MersenneTwister(49), n)
y = CEF + Noise
train, calibration, test = partition(eachindex(y), 0.4, 0.4)
########################################
LinearRegressor = @load LinearRegressor pkg=MLJLinearModels
model = LinearRegressor(fit_intercept = false) 
mach = machine(model, Tables.table(X), y)
fit!(mach, rows=train)
pr_y = predict(mach, rows=test)
########################################
conf_mach = conformal_machine(mach)
calibrate!(conf_mach, selectrows(X, calibration), y[calibration])
pr = predict(conf_mach, X[test,:]; coverage=0.95)
pr_lower = [pr[j][1][2][] for j in 1:length(test)]
pr_upper = [pr[j][2][2][] for j in 1:length(test)]
mean(pr_lower .<= y[test] .<= pr_upper)   # 0.94975
###########################################
# recover the predicted distribution
xt = [X[test[1],:] ;;]'
c_grid = .01:.001:0.99 
LB = []; UB = [];
for j in eachindex(c_grid)
    pr = predict(conf_mach, xt; coverage=c_grid[j] )
    push!(LB, pr[1][1][2][])
    push!(UB, pr[1][2][2][])
end
plot(legend=:topleft)
plot!(LB, (1.0 .- c_grid)/2.0, lab="LB, quantile")
plot!(UB, (c_grid[end]/2.0) .+ (c_grid)/2.0, lab="UB, quantile")
plot!([pr_y[1]], seriestype = :vline, lab="y prediction, median", color="red") 

#AZ: recover the predicted distribution? 
xt = [X[test[1],:] ;;]'
c_grid = .01:.01:0.99 
LB = []; UB = [];
for j in eachindex(c_grid)
    pr = predict(conf_mach, xt; coverage=c_grid[j] )
    push!(LB, pr[1][1][2][])
    push!(UB, pr[1][2][2][])
end
plot(legend=:topleft)
plot!(LB, 1.0 .- c_grid, lab="LB at %-ile")
plot!(UB, 1.0 .- c_grid, lab="UB at %-ile")
plot!([pr_y[1]], seriestype = :vline, lab="y prediction point estimate", color="red")