function combine_imprecise_priors_LR
% see: Morrison GS, Enzinger E (2016) Position Paper: What should a forensic practitioner’s likelihood ratio be? Science & Justice Virtual Special Issue on on measuring and reporting the precision of forensic likelihood ratios
% tested on Matlab version 8.6.0.267246 (R2015b)

x = (-7:.1:10)';

y_prior = f_prior(x);

y_LR = f_LR(x);

[y_posterior, xx] = convolve(x, y_prior, y_LR);

II_threshold_0 = xx >= 0;
p_threshold_0 = sum(y_posterior(II_threshold_0)) / sum(y_posterior)


figure(1);
plot(x, y_prior, '-r');
hold on
plot(x, y_LR, '-g');
plot(xx, y_posterior, '-b');
yy = [0, 0.6];
xlim([x(1) x(end)]);
set(gca, 'XTick', x(1):x(end), 'YLim', yy);
plot([0,0],yy, '-k');
hold off

end


function y = f_prior(x)
    boundary = [-6, -5, -3, -2];
    length_x = length(x);
    y = NaN(length_x,1);
    for Ix = 1:length_x
        z = x(Ix);
        if z <boundary(1)
            y(Ix) = 0;
        elseif z >= boundary(1) && z <boundary(2)
            y(Ix) = (z-boundary(1)) / (boundary(2)-boundary(1));
        elseif z >= boundary(2) && z <boundary(3)
            y(Ix) = 1;
        elseif z >= boundary(3) && z < boundary(4)
            y(Ix) = (boundary(4)-z) / (boundary(4)-boundary(3));
        elseif z >= boundary(4)
            y(Ix) = 0;
        end
    end
    y = y/3;
end


function y = f_LR(x)
    y = normpdf(x,6,.7653);
end


function [z, x_out] = convolve(x, y1, y2)
    length_x = length(x);
    zz = zeros(length_x, length_x);
    xz = zz;
    x2 = int16(x*10);
    x1 = x2;
    for Ix = 1:length_x
        xz(:,Ix) = x1 + x2;
        zz(:,Ix) = y1 .* y2;
        x1 = circshift(x1, 1);
        y1 = circshift(y1, 1);
    end
    unique_xz = unique(xz);
    length_unique_xz = length(unique_xz);
    z = zeros(length_unique_xz, 1);
    for Ix = 1:length_unique_xz
        II = xz == unique_xz(Ix);
        z(Ix) = sum(zz(II));
    end
    z = 10*z/sum(z);
    x_out = single(unique_xz)/10;
    II_delete = x_out < x(1) | x_out > x(end);
    z(II_delete) = [];
    x_out(II_delete) = [];
end