First program :

proc iml;

options nonotes;

a=24.953;

b=4;

cal1=1+(1/a);

cal2=1+(2/a);

toi=(sqrt(gamma(cal2)-((gamma(cal1))##2)))/gamma(cal1);

print toi;

* The function module;

start f_weib2(x) global(g_x,thet); /* x[1]=b and x[2]=a and g_x is global variable, which can be x1 or x2*/

*The log-LH fun. using p obs of the parameters is the objective function to be maximized to obtain MLEs;

N = ncol(g_x);

temp = g_x;

sum1 = sum(log(temp));

sum2=sum((temp / x[1])##x[2]);

f = N*log(x[2]) - N*x[2]*log(x[1]) + (x[2]-1)*sum1 - sum2;

return(f);

finish f_weib2;

start g_weib2(x) global(g_x,thet); /* x[1]=b and x[2]=a */

N = ncol(g_x);

g = j(1,2,0.);

temp = g_x;

sum1 = sum(log(temp));

sum2=sum((temp / x[1])##x[2]);

sum3=sum(((temp / x[1])##x[2])#(log(temp / x[1])));

g[1] = ((-N * x[2]) / x[1]) + (sum2 * (x[2] / x[1]));/*the result of drivative MLE to b*/

g[2] = (N/ x[2]) - (N * log(x[1])) + sum1 - sum3;/*the result of drivative MLE to a*/

return(g);

finish g_weib2;/* end module */

NumSamples=5000;

N=15;

r=j(NumSamples,1);

s=j(NumSamples,1);

e=j(NumSamples,1);

*const=j(NumSamples,1,1);

count=0;

na=5000;

nb=5000;

Est1 = j(NumSamples,2); /* allocate space to store parameter estimates */

Est2 = j(NumSamples,2); /* allocate space to store parameter estimates */

x1 = j(1,15); /* allocate (1 x p) matrix*/

x2= j(1,15); /* each row is sample of size N */

*Set up constraints;

con = { 1.e-6 1.e-6 ,

. . }; /*lower bounds of parms >0 (1.e-6) and upper bounds of parms < infinity (.)*/

*Call an optimization routine, If you specify first-order derivatives analytically

with the "grd" module argument, you can drastically reduce the computation time

for numerical second-order derivatives.;

optn = {1 0}; /* 1: find max of a function and 0: print no information*/

do t=1 to na;

call randseed(12345);

call randgen(x1, 'weibull', a, b); /* fill it*/

g_x = x1;

x = a||b; /* initial guess for solution (a, b) */

*f=f_weib2(x);

*g=g_weib2(x);

call nlptr(rc, result, "f_weib2", x, optn, con,,,,"g_weib2"); /* i_th estiamtes */

Est1[t,] = result; /* save t_th estimates */

wd1=Est1[t,2];

do k= 1 to NumSamples;

call randgen(x2, 'weibull',1, 1); /* k_th sample */

g_x = x2;

x = 1||1; /* initial guess for solution (a, b) */

*f=f_weib2(x);

*g=g_weib2(x);

call nlptr(rc, result, "f_weib2", x, optn, con,,,,"g_weib2"); /* i_th estiamtes */

Est2[k,] = result; /* save i_th estimates */

wd2=Est2[k,2];

ratio=wd1/wd2;

gam1=1+(1/ratio);

gam2=1+(2/ratio);

cove=(sqrt(gamma(gam2)-(gamma(gam1))##2))/gamma(gam1);

r[k,1]=cove;

end; /*end of i loop*/

e=r;

r[rank(r),]=e;

low=r[125,];

up=r[4875,];

wid=up-low;

s[t,]=wid;

*Check if toi in the constructed FGCIs;

if low<=toi then if toi<=up then count=count+1;

end;/*end of t loop*/

percent=count/5000;

print percent;

q=1-percent;

print q;

meanwid=s[:,];

print meanwid;

quit;

run;

Second program:

proc iml;

options nonotes;

a=24.953;

b=4;

cal1=1+(1/a);

cal2=1+(2/a);

toi=(sqrt(gamma(cal2)-((gamma(cal1))##2)))/gamma(cal1);

print toi;

* The function module;

start f_weib2(x) global(g_x,thet); /* x[1]=b and x[2]=a and g_x is global variable, which can be x1 or x2*/

*The log-LH fun. using p obs of the parameters is the objective function to be maximized to obtain MLEs;

N = ncol(g_x);

temp = g_x;

sum1 = sum(log(temp));

sum2=sum((temp / x[1])##x[2]);

f = N*log(x[2]) - N*x[2]*log(x[1]) + (x[2]-1)*sum1 - sum2;

return(f);

finish f_weib2;

start g_weib2(x) global(g_x,thet); /* x[1]=b and x[2]=a */

N = ncol(g_x);

g = j(1,2,0.);

temp = g_x;

sum1 = sum(log(temp));

sum2=sum((temp / x[1])##x[2]);

sum3=sum(((temp / x[1])##x[2])#(log(temp / x[1])));

g[1] = ((-N * x[2]) / x[1]) + (sum2 * (x[2] / x[1]));/*the result of drivative MLE to b*/

g[2] = (N/ x[2]) - (N * log(x[1])) + sum1 - sum3;/*the result of drivative MLE to a*/

return(g);

finish g_weib2;/* end module */

/*module for sampling with replacement*/

start SampleReplace(A, n, k);

m = j(n, k); /* allocate result matrix */

call randgen(m, "Uniform"); /* fill with random U(0,1) */

m= ceil(nrow(A)*ncol(A)*m); /* integers 1,2,...,ncol(A)*/

return(shape(A[m], n)); /* reshape and return */

finish;

x2= j(1,15); /* each row is sample of size N */

Est2 = j(10000,2); /* allocate space to store parameter estimates */

wd3=j(10000,1);

*Set up constraints;

con = { 1.e-6 1.e-6 ,

. . }; /*lower bounds of parms >0 (1.e-6) and upper bounds of parms < infinity (.)*/

*Call an optimization routine, If you specify first-order derivatives analytically

with the "grd" module argument, you can drastically reduce the computation time

for numerical second-order derivatives.;

optn = {1 0}; /* 1: find max of a function and 0: print no information*/

do k= 1 to 10000;

call randseed(12345);

call randgen(x2, 'weibull',1, 1); /* k_th sample */

g_x = x2;

x = {1 1}; /* initial guess for solution (a, b) */

*f=f_weib2(x);

*g=g_weib2(x);

call nlpnra(rc, result, "f_weib2", x, optn, con,,,,"g_weib2"); /* i_th estiamtes */

Est2[k,] = result; /* save i_th estimates */

end;

wd3=Est2[,2];

NumSamples=5000;

r=j(NumSamples,1);

s=j(NumSamples,1);

e=j(NumSamples,1);

count=0;

na=5000;

nb=5000;

Est1 = j(NumSamples,2); /* allocate space to store parameter estimates */

x1 = j(1,15); /* allocate (1 x p) matrix*/

do t=1 to na;

call randseed(12345);

call randgen(x1, 'weibull', a, b); /* fill it*/

g_x = x1;

x = a||b; /* initial guess for solution (a, b) */

*f=f_weib2(x);

*g=g_weib2(x);

call nlpnra(rc, result, "f_weib2", x, optn, con,,,,"g_weib2"); /* i_th estiamtes */

Est1[t,] = result; /* save t_th estimates */

wd1=Est1[t,2];

x3 = wd3;

call randseed(54321);

wd2 = SampleReplace(x3, 5000, 1);

ratio=wd1/wd2;

gam1=1+(1/ratio);

gam2=1+(2/ratio);

cove=(sqrt(gamma(gam2)-(gamma(gam1))##2))/gamma(gam1);

e=cove;

cove[rank(cove),]=e;

low=cove[125,];

up=cove[4875,];

wid=up-low;

s[t,]=wid;

*Check if toi in the constructed FGCIs;

if low<=toi then if toi<=up then count=count+1;

end;/*end of t loop*/

percent=count/5000;

print percent;

q=1-percent;

print q;

meanwid=s[:,];

print meanwid;

quit;

run;