/**********************************************************************/


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;                                                    


/**********************************************************************/