#maple source file 
#  - to solve problem 6.10 of Efron & Tibshirani
#  - compare timings using Maple function time(...)
#
# Composed By: A.I. McLeod (aim@uwo.ca)
# Date: August 1997
#
#To use do the following:
#
#1. Start maple
#
#2. Input this file with read command,
#
#      read `C:\\maple\\prob610.mpl`;
#
#3. Execute:
#
#      prob610();
#
#4. To get cpu timing:
#
#      time(prob610());
#
#ILLUSTRATIVE TIMINGS WITH MAPLE V RELEASE 4
#Thinkpad, Pentium-120: 117.09 sec
#PC Pentium Pro:
#M40, Risc 6000 (lexis): 66.62
#IBM Risc 6000
#
interface(verboseproc=2);
readlib(showtime);
with(combinat, multinomial);
prob610 := proc()
  local p, x, z, j, k, i, ii, tm, probab, X1, X2, se;
  p   :=  evalf(8^(-8));
  x   :=  array([1.2,3.5,4.7,7.3,8.6,12.4,13.8,18.1]);
  z   :=  array([0$8]);
  j   :=  generate_partition([0$8]);
  X1  := 0; X2 := 0;
  while ([ ] <> j) do 
    k := 1;
    for i to 8 do for ii to j[i] do z[k] := x[i]; k := k+1 od; od;
    tm := trimmean(convert(z, list), 0.25);
    probab := p*multinomial(8,j[1],j[2],j[3],j[4],j[5],j[6],j[7],j[8]);
    X1 := X1 + probab*tm;
    X2 := X2 + probab*tm^2;
    j := generate_partition(j);
  od;
  se := sqrt(X2-X1^2);
  [se, X1];
end;
generate_partition := proc(j)
#generate all non-negative integer solutions to the equation
# j_1+j_2+...+j_n = n
#Algorithm:
#Step 1. Start with j == (n,0,...,0) = (j_1, ..., j_n)
#Step 2. Find the largest index, i, such that j_i + j_(i+1) + ... + j_n < n.
#Step 3. Increment j_i and then set everything to its left to 0 except for
#        j_1 which is set to the largest possible value.
#Step 4. Repeat Steps 2 & 3 until (0,...,n) is reached. This occurs when
#        there is no value of j_i from Step 2 which can be found.
#
  local i, n, k, J, j1, jsum;
  J := nops(j);
  k := array([0$J]);
  n := 0; 
  for i to J do n := n+j[i]; od;
  if evalb(n = 0) then k[1] := J; RETURN(convert(k,list));
  elif evalb(j[J]=J) then RETURN([ ]);
  else
    jsum := 0;
    j1 := J;
    while evalb(jsum+j[j1] < n) do jsum := jsum+j[j1]; j1 := j1-1 od;
    k[j1+1] := j[j1+1]+1;
    for i from (j1+2) to J do k[i] := j[i] od;
    k[1] := J-jsum-1;
  fi;
  convert(k,list);
end;
trimmean := proc(x, p) #evaluates p-trimined mean for list x
  local n, ip, sx, tm, i;
  n   := nops(x);
  ip  := round(n*p);
  sx  := sort(x);
  tm  := 0;
  for i from ip+1 to n-ip do tm := tm + sx[i] od;
  tm/(n-2*ip);
end;
on();