function [output] = f_loss_QR(X,mod_params,loss_params)
    
    lower = loss_params.lower;
    upper = loss_params.upper;
    alpha = loss_params.alpha;
    beta  = loss_params.beta;
    w     = loss_params.w;

    yforeq            = X*mod_params;
    samples           = QR1_adjusted(0.05:0.05:0.95,yforeq);

    % Set integration bounds safely
    lb = max(min(samples), lower -5*std(samples));
    ub = min(max(samples), upper + 5*std(samples));

    % KDE once
    xgrid          = linspace(lb, ub, 400);
    [fgrid] = ksdensity(samples, xgrid, 'Function', 'pdf');
    %[xgrid, fgrid] = ksdensity(samples, 'Function', 'pdf', 'NumPoints', 1000);
    f_interp = @(x) interp1(xgrid, fgrid, x, 'linear', 0);
    %f_interp = @(x) arrayfun(@(xi) interp1(xgrid, fgrid, xi, 'linear', 0), x);

    % Define integrands
    DRf  = @(x) ((lower - x).^alpha) .* f_interp(x);
    EIRf = @(x) ((x - upper).^beta) .* f_interp(x);

    %DR  = quadgk(DRf, lb, lower);
    %EIR = quadgk(EIRf, upper, ub);

    %DRf   = @(x) (((lower-x)).^alpha).*(ksdensity(samples,x,'Function','pdf'));
    %EIRf  = @(x) (((x-upper)).^beta).*(ksdensity(samples,x,'Function','pdf'));

    %DR  = quadgk(DRf,-Inf,lower);
    %EIR = quadgk(EIRf,upper,Inf);
    DR  = quadgk(DRf,lb,lower, 'MaxIntervalCount', 1e5, 'AbsTol',1e-6, 'RelTol',1e-3);
    EIR = quadgk(EIRf,upper,ub, 'MaxIntervalCount', 1e5, 'AbsTol',1e-6, 'RelTol',1e-3);
    BR  = w.*EIR - (1-w).*DR;

    output = [DR,EIR,BR];
end       





 