An Introduction to Generalized Linear Models by Annette J. Dobson Chapter 1: Introduction

1.6.5 Example: Tropical cyclones

Table 1.2 on page 13.

data table1_2;
  input years $1-7	season	number;
datalines;
1956/7	1	6
1957/8	2	5
1958/9	3	4
1959/60	4	6
1960/1	5	6
1961/2	6	3
1962/3	7	12
1963/4	8	7
1964/5	9	4
1965/6	10	2
1966/7	11	6
1967/8	12	7
1968/9	13	4
;
run;

Figure 1.1 on page 14.

%let start = 3;
%let end = 8;
proc iml;
  use table1_2;
  read all var {number} into y ;  
  n = nrow(y) ;  
  c = j(&end-&start+1, 2, 0);  
  do i = 1 to &end - &start + 1;
     c[i, 1]= sum(y)*log(i+&start - 1) - n*(i+&start - 1);
	 c[i, 2] = i+&start - 1;
  end;
  create fig1_1 from c[colname={'lstar' 'theta'}];
  append from c;
quit;
symbol i = spline v=none;
axis1 order = (35 to 55 by 5) minor = none label=none;
proc gplot data = fig1_1;
  plot lstar*theta /vaxis = axis1 hminor=0;
run;
quit;

Table 1.3 on page 15. We use proc iml to implement the bisection calculation.

%let criterion = .00001;
%let left = 5;
%let right = 6;
%let itnum = 20;
proc iml;
  use table1_2;
  read all var {number} into y ; /*read the number of count into a matrix*/
  n = nrow(y) ; /*number of rows*/
  c = j(&itnum, 3, 0); /*create a matrix containing informatin for Table 1.3*/
  left = &left;
  right = &right;

  c[1,1] = 1;
  c[1, 2] = left;
  c[1, 3]= sum(y)*log(left) - n*left;

  c[2,1] = 2;
  c[2, 2] = right;
  c[2, 3]= sum(y)*log(right) - n*right;
   lb = c[1,3];
   ub = c[2,3];
  diff = ub - lb;
   	do i = 1 to &itnum while(abs(diff)>&criterion);
    	    mid = (left + right)/2;
			lstar = sum(y)*log(mid) - n*mid;
		c[i+2, 1] = i+2;
		c[i+2, 2] = mid;
		c[i+2, 3] = lstar;
    		if ub > mid & mid > lb then do;
        		lb = lstar; 
				left = mid;
			end;
	    	else if lb > mid & mid > ub then do ;
				right = mid;
				ub = lstar;
			end;
			else if lb < ub then do;
				left = mid;
				lb = lstar;
			end;
			else do;
				right = mid;
				ub = lstar;
			end;
			diff = ub - lb;
  end;
	 create c  var {k theta_k lstar};
	 append from c;
quit;
proc print data = c noobs;
  where k ~=0;
  format lstar 8.5;
run;

 K    THETA_K       LSTAR
 1    5.00000    50.87953
 2    6.00000    51.00668
 3    5.50000    51.24186
 4    5.75000    51.19239
 5    5.62500    51.23491
 6    5.56250    51.24293
 7    5.53125    51.24355
 8    5.54688    51.24352
 9    5.53906    51.24361
10    5.53516    51.24359
11    5.53711    51.24360