modelProbComp.m

function prob = modelProbComp(Sims, Rexps, Rdotexps, Rsigma, Rdotsigma, tmatrix, maxTime)
    % Return a single scalar "probability-like" value in (realmin, 1] by
    % averaging per-trial likelihoods in log-space to avoid underflow.

    % basic shape checks (defensive)
    if size(Sims,2) < 3
        prob = realmin; return;
    end

    nTrials = size(Rexps, 2);
    logL = [];  % collect per-trial log-likelihoods

    for col = 1:nTrials
        idx = find(tmatrix(:, col) <= maxTime);
        if isempty(idx) || max(idx) > size(Sims, 1)
            continue;
        end

        R_exp    = Rexps(idx, col);
        Rdot_exp = Rdotexps(idx, col);
        R_sim    = Sims(idx, 2);
        Rdot_sim = Sims(idx, 3);

        sigmaR    = max(Rsigma(idx),    1e-12);
        sigmaRdot = max(Rdotsigma(idx), 1e-12);

        % per-trial Gaussian log-likelihood
        ll_R   = -0.5 * sum( log(2*pi*sigmaR.^2)    + ((R_exp  - R_sim   ).^2) ./ (sigmaR   .^2) );
        ll_Rd  = -0.5 * sum( log(2*pi*sigmaRdot.^2) + ((Rdot_exp - Rdot_sim).^2) ./ (sigmaRdot.^2) );
        ll_tot = ll_R + ll_Rd;

        if isfinite(ll_tot)
            logL(end+1,1) = ll_tot; 
        end
    end

    if isempty(logL)
        prob = realmin;
        return;
    end

    % log-mean-exp across trials (i.e., average likelihood, not sum)
    m   = max(logL);
    lme = m + log(mean(exp(logL - m)));   % log(mean(exp(logL)))
    prob = max(exp(lme), realmin);

    % simple non-physical guard (optional—keep if you like it)
    if min(Sims(:, 2)) >= Sims(end, 2) || isnan(prob)
        prob = realmin;
    end
end