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