1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
|
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
|
