/**********************************************************************/
data tvals;
input t1-t8 @@;
cards;
18.0 15.9 24.8 19.7 20.0 12.4 18.5 21.2
;
data test;
set tvals;
array t{8} t1-t8;
array k{8};
array m{8};
array d{8};
do i=1 to 8;
k{i}=1/(t{i});
ksum+k{i};
end;
/* Harmonic mean calculation */
harmean=8/ksum;
/* Calculation of n harmonic means derived */
/* from samples of n-1 */
do i=1 to 8;
m{i}=7/(ksum-k{i});
msum+m{i};
end;
/* Mean of n harmonic means derived from */
/* samples of n-1 */
mean_nm1=msum/8;
/* Sum of squared deviations */
do i=1 to 8;
d{i}=(m{i}-mean_nm1)**2;
dsum+d{i};
end;
/* Jackknife standard deviation */
std=sqrt(dsum*7/8);
keep harmean mean_nm1 std;
title "Harmonic Mean and Standard Deviation Via Array";
proc print data=test;
run;
/**********************************************************************/
data tvals;
input t1-t8 @@;
cards;
18.0 15.9 24.8 19.7 20.0 12.4 18.5 21.2
;
proc iml;
use tvals;
read into t;
close tvals;
numt=ncol(t);
/* Inverse of half-lives */
k=1/t;
/* Harmonic mean calculation */
harmean=numt/(sum(k));
/* Calculation of n harmonic means derived */
/* from samples of n-1 */
h=j(numt,1,0);
do i=1 to numt;
h[i,1] =(numt-1)/(sum(k)-k[1,i]);
end;
/* Mean of n harmonic means derived from */
/* samples of n-1 */
mean_nm1=sum(h)/numt;
/* Jackknife standard deviation */
std=sqrt((ssq(h-mean_nm1))*(numt-1)/(numt));
/* Create data set of values */
create imldata
var{harmean mean_nm1 std};
append;
quit;
title "Harmonic Mean and Standard Deviation Via IML";
proc print data=imldata;
run;
/**********************************************************************/
data halflife;
input subject hl1 hl2 hl3 hl4 @@;
cards;
1 18.0 21.2 18.0 21.2
2 15.9 18.5 15.9 15.9
3 . . . .
4 24.8 12.4 24.8 24.8
5 19.7 20.0 19.7 12.4
6 20.0 19.7 20.0 20.0
7 12.4 24.8 12.4 19.7
8 . . . .
9 18.5 15.9 18.5 18.5
10 21.2 18.0 21.2 18.0
;
/* Macro variable for half-lives to be analyzed. */
%let vars=HL1 HL2 HL3 HL4;
proc iml;
use halflife;
read all var{&vars} into hlmat;
close halflife;
nobs=nrow(hlmat);
nvar=ncol(hlmat);
/* Code to find the number of nonmissing values. */
/* The number of nonmissing values are then used */
/* as the sample size for all of the relevant */
/* calculations. */
zero=j(nobs,nvar,0);
test=hlmat>zero;
nsamp=test[+,];
/* Inverse of half-lives */
k=1/hlmat;
/* Harmonic mean calculation */
harmean=1/k[:,];
/* Calculation of n harmonic means derived */
/* from samples of n-1 */
h=j(nobs,nvar,0);
do j=1 to nvar;
do i=1 to nobs;
h[i,j]=((nsamp[1,j])-1)/
((sum(k[,j]))-k[i,j]);
end;
end;
/* Mean of n harmonic means derived from */
/* samples of n-1 */
mean_nm1=h[:,];
/* Jackknife standard deviation */
dev=h-repeat(mean_nm1,nobs,1);
ssquares=dev[##,];
var=ssquares#(nsamp-j(1,nvar,1));
var=var/nsamp;
std=sqrt(var);
/* Name variables for data sets */
names={harmean mean_nm1 std n_sample};
id={&vars};
id=id`;
idname={ID};
final=harmean`|| mean_nm1`|| std`
|| nsamp`;
/* Create data set */
create final from final [colname=names];
append from final;
create id from id [colname=idname];
append from id;
quit;
data all;
merge id final;
run;
title1 "Harmonic Means and Standard Deviations";
proc print data=all noobs;
run;
/**********************************************************************/