Quantile vs HDI?

[DominiqueMakowski]

In correlation_test for instance, you use quantiles q = quantile(ts,[0.005,0.025,0.5,0.975,0.95]) to extract what I think is credible intervals. Is there a reason to use quantiles instead of HDI?

[tszanalytics]

Hey Dominique! I just found all these comments! No, no good reason, just the first attempt. I have add bounds to keep it in [-1,1], but haven't thought about HDI yet. Is there some code available already to do that, or should I think about how to do that?

I see that mamba.jl has a function, but it looks like it has some limitations. I believe Kruschke has some R code that would be very easy to port. Let me know if you had something in mind, otherwise I will have a go at it from scratch.

[tszanalytics]

Actually, for most of the tests implemented so far, since the distribution is unimodal and symmetric, the quantiles give the HPDI. For the correlation test, this will be true as long as the density is not near the -1 or 1 boundaries, in which case it is likely to be the area in one tail, so a one-sided quantile value will work. I have figured out a good way to impose the boundary conditions on the posterior kernel density plot, so I will have a go at the HPDI computation for that soon. I have a newer version of the bayes_correlation function (though I haven't renamed it to that yet - I will do so before uploading it), but I will try to incorporate the boundary constraints and compute HPDIs before I upload a newer version.

[tszanalytics]

To do asap:

1. Add ASH code for better kernel density evaluations.

2. Change function names to add bayes_ prefix.

3. Determine if separate HPDI function is needed yet (I think everything, with the possible exception of bayes_corr is unimodal and symmetric). Add a fn. if needed.

4. Naming: bayes_corr (one or two r's?), bayes_correlation, bayes_corr_test, bayes_cortest? Pick a name - I will ask a few people to get their reaction!

[fm444fm]

a simple HDI in Julia, adapted from a python equivalent:

function HDI_from_MCMC(posterior_samples; credible_mass=0.95)
    # Computes highest density interval from a sample of representative values,
    # estimated as the shortest credible interval
    # Takes Arguments posterior_samples (samples from posterior) and credible mass (normally .95)
    # Originally from https://stackoverflow.com/questions/22284502/highest-posterior-density-region-and-central-credible-region
    # Adapted to Julialang
    sorted_points = sort(posterior_samples)
    ciIdxInc = Int(ceil(credible_mass * length(sorted_points)))
    nCIs = length(sorted_points) - ciIdxInc
    ciWidth = repeat([0.0],nCIs)
    for i in range(1, stop=nCIs)
        ciWidth[i] = sorted_points[i + ciIdxInc] - sorted_points[i]
    end
    HDImin = sorted_points[findfirst(isequal(minimum(ciWidth)),ciWidth)]
    HDImax = sorted_points[findfirst(isequal(minimum(ciWidth)),ciWidth)+ciIdxInc]
    return(HDImin, HDImax)
end