function logL = custom_logLikePaperV2(params,y)
nchain = size(params,1);
nexp = size(y,1);
logL=zeros(nchain,1);
for i1 = 1:nchain
    for i2 = 1:nexp
        logIntMC = log(MC_int(params(i1,:),y(i2,:)));
        logL(i1) = logL(i1)+logIntMC;
    end
end
end
function logl_MC = MC_int(params,Y)
%% data for the montecarlo integration 
ran = [0 1; 0 1];   % integration range
abstol = 1e-6;      % abs error tolerance 
reltol = 1e-3;      % relative error tolerance
timelimit = inf;    % sets time limit (in sec) on integration
N = 5e5;            % ma no of function evaluations
n = 1e3;            % no of function calls for iteration
    %% solve loop
    t0 = clock;
    d = size(ran,1);    % dimension
    V=prod(range(ran'));
    E = NaN(1); 
    R = NaN(1);
    w = 1;
    points_used_locally = 0;
    s = zeros(1,2);
    t1 = inf;
    iter = 0;
    n_start = 1;
    Inv_var = inv(diag([1e-8,1e-8,1e-8]));
    while ~(E<abstol) && ~(R<reltol) && etime(clock,t0)<timelimit && sum(points_used_locally)<=N
        iter=iter+1;
        m = params(1); v = params(2)^2;
        mu = log((m^2)/sqrt(v+m^2));
        sigma = sqrt(log(v/(m^2)+1));
        E = lognrnd(mu,sigma,[n,1]);
        L = 1; h = 0.1; b = 0.1; I = 1/12*b*h^3; F = 30e3;
        disp = @(x) F.*x/2./(48*E*I).*(3*L^2-4*(x).^2);
        f = zeros(n,1);
        Y_forward = [disp(0.25) disp(0.5) disp(0.75)];
        for i1 = 1:n
            f(i1) = exp(-0.5*(Y-Y_forward(i1,:))*Inv_var*(Y-Y_forward(i1,:))');
        end
        s = s + [sum(f) sum(f.*f)];
        points_used_locally = n+points_used_locally;
        points_used_globally = n+points_used_locally;
        % evaluation of the integral
        % aggregate absolute error
        E=V/points_used_globally*sqrt( abs(sum(s(:,2))/points_used_globally - (sum(s(:,1))/points_used_globally).^2)+eps);
        %   full integral
        logl_MC = sum(1/points_used_locally .* s(:,1));
        %   aggregate relative error
        R = E/abs(logl_MC);
    end
end 