/* ed50.pgm computes the ed50 =-intercept/slope         */
/* and gives confidence limits from Fiellers theorem.   */
proc iml;
reset nolog nocenter;
/* print, 'ED50 and C.I.'; */

/* Put in alpha for confidence coeff. and df=degrees of freedom.   */

alpha=.05;
t=probit(1-alpha/2);       /* use this for normal critical values  */

/* Put in estimates of intercept and slope from binary regression  */

intercept=-6.3134 ;slope=4.3185;
ed50=-intercept/slope;r=ed50;

/* put in 2 by 2 covariance matrix of est. of  intercept and slope */

cov_mat={ 0.95423  -0.65166,
       -0.65166 0.45425};

print, alpha intercept slope ed50;
print, cov_mat;
der=j(2,1);
der[1]=-1/slope;
der[2]=intercept/(slope**2);
varrat=der`*cov_mat*der;     /* delta theorem var est. */
std=sqrt(varrat);
print,"ED50 and estimated std";
print,  ed50 std;

print,"Confidence Intervals from Fieller's Theorem";
c11=cov_mat[1,1];
c12=-cov_mat[1,2];
c22=cov_mat[2,2];
g=(t##2)*c22/(slope##2);
print,"G: must be less than 1.0 ";
print, g t;
ins=c11-2*r*c12+r*r*c22-g*(c11-c12*c12/c22);
plm=(t/slope)*sqrt(ins);
rl=(r-g*c12/c22-plm)/(1-g);
ru=(r-g*c12/c22+plm)/(1-g);
print,"ED50, Lower and Upper 1-alpha/2 Limits";
print,  ed50 rl ru;