function ROC = NormalROC(logbetas,us,vars,un,varn) % ROC = NormalROC(logbetas,us,vars,un,varn) % % Compute ROC curve given the parameters. Done % analytically for speed. % % The trick here is to solve for the values of x % such that log(l(x)) > logbetas. When the variances % are unequal, there can be two distinct regions. % I basically handle this by brute force, since there % are not too many distinct cases. % % This routine is not carefully tested and may % contain errors. % % xx/xx/xx dhb Wrote it. % 10/4/00 dhb Added caveat about bug possibilities, based % in part on disagreement between this routine % and some Monte Carlo simulations. % Set up room for answer [mbetas,null] = size(logbetas); if ( mbetas == 1 && null ~= 1) logbetas = logbetas'; mbetas = null; end ROC = zeros(mbetas,2); % By convention, will assume that us > un, switch if % false if ( un > us ) ut = us; vart = vars; us = un; vars = varn; un = ut; varn = vart; end % Set up the quadratic equation to find the zeros of % the likelyhood ratio. This expression is arrived at % by writing out the likelihood ratio for Normal random % variables, taking logarithms, and getting a quadratic % expression for the values of x that generate the desired % log beta. a = (vars-varn)*ones(mbetas,1); b = 2*(varn*us - vars*un)*ones(mbetas,1); c = (vars*(un^2)-varn*(us^2))*ones(mbetas,1) + ... varn*vars*log(varn/vars)*ones(mbetas,1) - ... 2*varn*vars*logbetas; if (vars == varn) x = (-1*c) ./ b; ROC(:,1) = ones(mbetas,1) - NormalCumulative(x,us,vars); ROC(:,2) = ones(mbetas,1) - NormalCumulative(x,un,varn); else % Solve for the zeros of the quatdratic discrim = sqrt( b.^2 - 4.*a.*c ); x1 = (-b + discrim) ./ (2*a); x2 = (-b - discrim) ./ (2*a); % Handle the case where ratio heads to +infinity if ( (vars-varn) > 0 ) % Imaginary roots mean 1,1 index = find( imag(discrim) ~= 0 ); [m,null] = size(index); if ( m ~= 0 ) ROC(index,1) = ones(m,1); ROC(index,2) = ones(m,1); end index = find( imag(discrim) == 0 ); [m,null] = size(index); if ( m ~= 0 ) ROC(index,1) = ones(m,1) - ... NormalCumulative(x1(index),us,vars) - ... NormalCumulative(x2(index),us,vars); ROC(index,2) = ones(m,1) - ... NormalCumulative(x1(index),un,varn) - ... NormalCumulative(x2(index),un,varn); end % Handle the case where the ratio heads to -infinity else % Imaginary roots mean 0,0 index = find( imag(discrim) ~= 0 ); [m,null] = size(index); if ( m ~= 0 ) ROC(index,1) = zeros(m,1); ROC(index,2) = zeros(m,1); end index = find( imag(discrim) == 0 ); [m,null] = size(index); if ( m ~= 0 ) ROC(index,1) = NormalCumulative(x2(index),us,vars) - ... NormalCumulative(x1(index),us,vars); ROC(index,2) = NormalCumulative(x2(index),un,varn) - ... NormalCumulative(x1(index),un,varn); end end end