% MATLAB SCRIPT FILE s322hw1.m  */
%    SOLUTION TO HW 1, ABOUT SEQEUNCES AND THEIR AVERAGES
%    DUE AUGUST 29, 1996
'Running MATLAB script file s322hw1.m' 

%    HERE IS THE SEQUENCE OF TRIALS THAT I DID TO SOLVE THIS
%    PROBLEM.  EACH ONE WAS CHOSEN AFTER VIEWING THE RESULTS OF THE LAST

format short ;
                %  THE DEFAULT OUTPUT FORMAT
format compact ;
                %  THIS MAKES MORE EFFICIENT USE OF THE SCREEN BY
                %  SUPRESSING LINE FEEDS

alpha = 1 ;
beta = .1 ;

jstart = 2 ;    %  STARTING VALUE
jstep = 1 ;     %  STEP TO TAKE IN CREATION OF SEQUENCE
nj = 20 ;       %  NUMBER IN SEQUENCE (20 KEEPS THE RESULT ON THE SCREEN)                                         

jend = jstart + jstep * (nj - 1) ;
                %  HERE IS ONE WAY TO SET THIS UP.  COULD ALSO CHOOSE
                %  JEND, WHICH THEN DETERMINES NJ.

vj = jstart:jstep:jend ;
                %  LETTING VARIABLE NAMES START WITH v CAN BE USEFUL
                %  TO REMIND THAT IT IS A VECTOR

vj = vj' ;      %  MAKE vj A COLUMN VECTOR

%vj' 
                %  IT IS VERY USEFUL TO CHECK THINGS DURING THE 
                %  DEVELOPMENT OF A SCRIPT, BY USING LINES LIKE 
                %  THIS (E.G. BY LOOKING AT THIS, I REALIZED I
                %  NEEDED THE (nj - 1), INSTEAD OF nj ABOVE)
                %  USUALLY THESE CAN DELETED AFTER THEY SERVE
                %  THEIR PURPOSE, BUT THIS ONE IS LEFT IN TO
                %  ILLUSTRATE THIS POINT.

vaj = vj.^(-alpha) ;
vbj = exp(-vj * beta) ;
                %  THIS IS THE CALCULATION OF THE SEQUENCES

'Part 1b: first try' 
[vj vaj vbj (vaj - vbj) (vaj ./ vbj)] 
                %  HERE THE SEQUENCES ARE PUT INTO A MATRIX FOR 
                %  VIEWING ON THE SCREEN

'Any key to continue' 
pause ;
                %  THESE LINES STOP THE PROGRAM, SO THAT ONE CAN
                %  LOOK AT THE RESULTS THUS FAR.  THE FIRST LINE
                %  WRITES A REMINDER AS TO HOW TO GO ON.

%    THESE RESULTS SHOW THAT THE SEQUENCES DON'T GO FAR ENOUGH YET.
%    THE NEXT STEP WOULD BE TO CHANGE THE PARAMETERS A BIT.  
%    ORDINARILY THIS IS DONE BY SIMPLE EDITING OF THE ABOVE.
%    BUT SO THAT YOU CAN SEE THE STEPS THAT I TOOK, THE ABOVE LINES 
%    WILL INSTEAD BE COPIED BELOW (WITH MOST THE COMMENTS REMOVED).

'Part 1b: 2nd try with bigger starting point' 
jstart = 23 ;    %  STARTING VALUE
jstep = 1 ;     %  STEP TO TAKE IN CREATION OF SEQUENCE
nj = 20 ;       %  NUMBER IN SEQUENCE (20 KEEPS THE RESULT ON THE SCREEN)                                         
jend = jstart + jstep * (nj - 1) ;
vj = jstart:jstep:jend ;
vj = vj' ;      %  MAKE vj A COLUMN VECTOR

vaj = vj.^(-alpha) ;
vbj = exp(-vj * beta) ;
[vj vaj vbj (vaj - vbj) (vaj ./ vbj)] 

'Any key to continue' 
pause ;

%    THE ABOVE SHOWS THAT THE FIRST TIME THAT aj > bj IS j = 36.
%    NOW PREPARE AN OUTPUT FILE, FOR PRINTING, AS DONE IN THE
%    SCRIPT s322eg1.m

'Part 1b: found answer, so print result to output file' 
fid = fopen('s322hw1.out','wt') ;
  titlstr1 = 'Output from the MATLAB Script s322hw1.m,   ' ;
  titlstr1 = [titlstr1,date] ;
  titlstr2 = '    Solution to Home Work Assignment 1, due August 29, 1996' ;
cntbytes = fprintf(fid,'%1s\n',titlstr1) ;
cntbytes = fprintf(fid,'%1s\n',titlstr2) ;

  headstr1 = 'Problem 1, part b' ;
  headstr2 = '           j         a_j          b_j     (a_j - b_j)' ;
cntbytes = fprintf(fid,'\n\n %1s\n',headstr1) ;
cntbytes = fprintf(fid,'%1s\n',headstr2) ;

  j = 35 ; aj = j.^(-alpha) ; bj = exp(-j * beta) ;
  outm = [j; aj; bj; (aj-bj)] ;
  j = 36 ; aj = j.^(-alpha) ; bj = exp(-j * beta) ;
  outm = [outm; [j; aj; bj; (aj-bj)]] ;
cntbytes = fprintf(fid,'%12.0d %12.4f %12.4f % 12.4f\n',outm) ;



%    NOW MOVE ON TO PART C.  SOLVE AS ABOVE, EXCEPT NOW USE BIGGER 
%    STEP BECAUSE THIS WILL SURELY BE FURTHER OUT

alpha = 10 ;
beta = .0001 ;

'Part 1c: first try' 
jstart = 100 ;    %  STARTING VALUE
jstep = 100 ;     %  STEP TO TAKE IN CREATION OF SEQUENCE
nj = 20 ;       %  NUMBER IN SEQUENCE (20 KEEPS THE RESULT ON THE SCREEN)                                         
jend = jstart + jstep * (nj - 1) ;
vj = jstart:jstep:jend ;
vj = vj' ;      %  MAKE vj A COLUMN VECTOR

vaj = vj.^(-alpha) ;
vbj = exp(-vj * beta) ;
[vj vaj vbj (vaj - vbj) (vaj ./ vbj)] 

'Any key to continue' 
pause ;


%    THOSE DIDN'T GO FAST ENOUGH, TRY MUCH BIGGER VALUES

'Part 1c: 2nd try with bigger starting point' 
jstart = 10^5 ;    %  STARTING VALUE
jstep = 10^5 ;     %  STEP TO TAKE IN CREATION OF SEQUENCE
nj = 20 ;       %  NUMBER IN SEQUENCE (20 KEEPS THE RESULT ON THE SCREEN)                                         
jend = jstart + jstep * (nj - 1) ;
vj = jstart:jstep:jend ;
vj = vj' ;      %  MAKE vj A COLUMN VECTOR

vaj = vj.^(-alpha) ;
vbj = exp(-vj * beta) ;
[vj vaj vbj (vaj - vbj) (vaj ./ vbj)] 

'Any key to continue' 
pause ;

%    THAT OUTPUT FORMAT WAS HOPELESS, LET'S TRY ANOTHER

'Part 1c: 2nd try, but better output format' 
format short e ;
[vj vaj vbj (vaj - vbj) (vaj ./ vbj)] 
'Any key to continue' 
pause ;


%    NOW SEE ANSWER IS BETWEEN 1.4E006 AND 1.5E006, SO "ZOOM IN"

'Part 1c: 3nd try, zooming in on answer' 
jstart = 1.4*10^6 ;    %  STARTING VALUE
jstep = 10^4 ;     %  STEP TO TAKE IN CREATION OF SEQUENCE
nj = 10 ;       %  NUMBER IN SEQUENCE (20 KEEPS THE RESULT ON THE SCREEN)                                         
jend = jstart + jstep * (nj - 1) ;
vj = jstart:jstep:jend ;
vj = vj' ;      %  MAKE vj A COLUMN VECTOR

vaj = vj.^(-alpha) ;
vbj = exp(-vj * beta) ;
[vj vaj vbj (vaj - vbj) (vaj ./ vbj)] 

'Any key to continue' 
pause ;


%    NOW A COUPLE OF VERY SIMILAR ZOOM STEPS, WE GET TO:

'Part 1c: 4th try, after several such zoom steps' 
jstart = 1.4163*10^6 ;    %  STARTING VALUE
jstep = 10^1 ;     %  STEP TO TAKE IN CREATION OF SEQUENCE
nj = 10 ;       %  NUMBER IN SEQUENCE (20 KEEPS THE RESULT ON THE SCREEN)                                         
jend = jstart + jstep * (nj - 1) ;
vj = jstart:jstep:jend ;
vj = vj' ;      %  MAKE vj A COLUMN VECTOR

vaj = vj.^(-alpha) ;
vbj = exp(-vj * beta) ;
[vj vaj vbj (vaj - vbj) (vaj ./ vbj)] 

'Any key to continue' 
pause ;


%    PROBLEM HERE IS CAN'T READ THE VALUE OF J, 
%    AGAIN ADJUST THE FORMAT

format long e ;
'Part 1c: 4th try, with better format, so can see j''s' 
                %  NOTE THAT DOUBLE '' SAYS "PUT A QUOTE IN THE STRING"
jstart = 1.4163*10^6 ;    %  STARTING VALUE
jstep = 10^1 ;     %  STEP TO TAKE IN CREATION OF SEQUENCE
nj = 10 ;       %  NUMBER IN SEQUENCE (20 KEEPS THE RESULT ON THE SCREEN)                                         
jend = jstart + jstep * (nj - 1) ;
vj = jstart:jstep:jend ;
vj = vj' ;      %  MAKE vj A COLUMN VECTOR

vaj = vj.^(-alpha) ;
vbj = exp(-vj * beta) ;
[vj vaj vbj (vaj - vbj) (vaj ./ vbj)] 

'Any key to continue' 
pause ;


%    THAT RAN INTO TWO MUCH OUTPUT, SO CUT IT DOWN A LITTLE

'Part 1c: 4th try, with better format, and less output' 
[vj (vaj - vbj) (vaj ./ vbj)] 
'Any key to continue' 
pause ;


%    AFTER A LITTLE MORE OF THE ABOVE, WE SEE THAT THE FIRST 
%    TIME THAT aj > bj IS j = 36,  SO WRITE TO OUTPUT FILE

'Part 1c: found answer, so print result to output file' 

  headstr1 = 'Problem 1, part c' ;
  headstr2 = '         j         a_j          b_j     (a_j - b_j)' ;
cntbytes = fprintf(fid,'\n\n %1s\n',headstr1) ;
cntbytes = fprintf(fid,'%1s\n',headstr2) ;

  j = 1416360 ; aj = j.^(-alpha) ; bj = exp(-j * beta) ;
  outm = [j; aj; bj; (aj-bj)] ;
  j = 1416361 ; aj = j.^(-alpha) ; bj = exp(-j * beta) ;
  outm = [outm; [j; aj; bj; (aj-bj)]] ;
cntbytes = fprintf(fid,'%12.0d %12.4e %12.4e % 12.4e\n',outm) ;
          %  NOTE "EXPONENTIAL" FORMAT IS MORE USEFUL FOR THIS PART


%    NOW MOVE ON TO PART D.  SOLVE AS ABOVE, EXCEPT NOW USE BIGGER 
%    STEPS BECAUSE THIS WILL SURELY BE FURTHER OUT
%    BECAUSE THEY ARE SO SIMILAR TO THE STEPS ABOVE, ONLY A FEW OF
%    THE IMPORTANT ONES ARE SHOWN HERE

alpha = 1000 ;
beta = .0001 ;

'Part 1d: first try' 
jstart = 10^8 ;    %  STARTING VALUE
jstep = 10^8 ;     %  STEP TO TAKE IN CREATION OF SEQUENCE
nj = 20 ;       %  NUMBER IN SEQUENCE (20 KEEPS THE RESULT ON THE SCREEN)                                         
jend = jstart + jstep * (nj - 1) ;
vj = jstart:jstep:jend ;
vj = vj' ;      %  MAKE vj A COLUMN VECTOR

vaj = vj.^(-alpha) ;
vbj = exp(-vj * beta) ;
[vj (vaj - vbj) (vaj ./ vbj)] 

'Any key to continue' 
pause ;


%    THE PROBLEM WAS "UNDERFLOW", I.E. MATLAB IS TRYING TO WORK WITH
%    NUMBERS CLOSER TO 0 THAN IT CAN HANDLE (ABOUT 10^(-300)), SO IT
%    SIMPLE REPLACES THEM WITH THE "BEST APPROXIMATION IT CAN FIND",
%    WHICH IS ZERO ITSELF.  SINCE BOTH VALUES ARE 0, THIS RESULTS IN
%    THE RATIO BEING 0/0, WHICH MATLAB CALLS "NAN", FOR "NOT A NUMBER"
%    LET'S TACKLE THIS PROBLEM BY COMPARING LOGS.
           
'Part 1d: first try, log scale' 
jstart = 10^8 ;    %  STARTING VALUE
jstep = 10^8 ;     %  STEP TO TAKE IN CREATION OF SEQUENCE
nj = 20 ;       %  NUMBER IN SEQUENCE (20 KEEPS THE RESULT ON THE SCREEN)                                         
jend = jstart + jstep * (nj - 1) ;
vj = jstart:jstep:jend ;
vj = vj' ;      %  MAKE vj A COLUMN VECTOR

vlogaj = -alpha * log(vj) ;
vlogbj = -vj * beta ;
[vj (vlogaj - vlogbj)] 

'Part 1c: found answer, so print result to output file' 

  headstr1 = 'Problem 1, part d' ;
  headstr2 = '         j         a_j          b_j     (a_j - b_j)' ;
cntbytes = fprintf(fid,'\n\n %1s\n',headstr1) ;
cntbytes = fprintf(fid,'%1s\n',headstr2) ;

  j = 190660024 ;
  logaj = -alpha * log(j) ;
  logbj = -j * beta ;
outm = [j; logaj; logbj; (logaj-logbj)] ;
  j = 190660025 ;
  logaj = -alpha * log(j) ;
  logbj = -j * beta ;
outm = [outm; [j; logaj; logbj; (logaj-logbj)]] ;
cntbytes = fprintf(fid,'%12.0d %12.4e %12.4e % 12.4e\n',outm) ;



%    NOW WORK ON PROBLEM 2 
valpha = [.01 .1 .5 .9] ;
            %  PUT THEM AS A ROW VECTOR
vn = [10 ; 100 ; 1000] ;
            %  PUT THEM AS A COLUMN VECTOR

%    FIRST CONSTRUCT MATRIX OF AVERAGES
n = vn(1) ;
vj = 2:1:n ;
            %  VECTOR OF j
malpha = ones(n-1,1) * valpha ;
            %  THIS MATRIX MULTIPLICATION GIVES MATRIX WITH
            %  CONSTANT COLUMNS, THAT ARE THE -alpha
mj = vj' * ones(1,4) ;
            %  THIS MATRIX MULTIPLICATION GIVES MATRIX WITH
            %  CONSTANT ROWS, THAT ARE THE j
maj = mj .^ (-malpha) ;
            %  MATRIX OF a_j TO AVERAGE COLUMS
mavg = mean(maj) ;
            %  WILL BUILD THIS INTO MATRIX OF AVERAGES

%    THESE ARE REPEAT OF ABOVE LINES FOR n = 100
n = vn(2) ;
vj = 2:1:n ;
malpha = ones(n-1,1) * valpha ;
mj = vj' * ones(1,4) ;
maj = mj .^ (-malpha) ;
mavg = [mavg; mean(maj)] ;

%    THESE ARE REPEAT OF ABOVE LINES FOR n = 1000
n = vn(3) ;
vj = 2:1:n ;
malpha = ones(n-1,1) * valpha ;
mj = vj' * ones(1,4) ;
maj = mj .^ (-malpha) ;
mavg = [mavg; mean(maj)] ;

%    NOW CONSTRUCT MATRIX OF APPROXIMATIONS
mn = vn * ones(1,4) ;
malpha = ones(3,1) * valpha ;
mapprox = mn .^ (-malpha) ;
mapprox = mapprox ./ (1 - malpha) ;

%    NOW CONSTRUCT MATRICES OF BOUNDS
mterm = mn .^ (1 - malpha) - 1 ;
mterm = mterm ./ (1 - malpha) ;
mlobd = (mterm + mn .^ (-malpha) - 1) ./ (mn - 1) ;
mupbd = mterm ./ (mn - 1) ;

%    NOW WRITE TO OUTPUT FILE
'Part 2: Write result to output file' 

  headstr1 = 'Problem 2' ;
  headstr2 = '                      alpha = .01' ;
  headstr3 = '   n     Lower Bound      Value    Upper Bound    Estimate' ;
cntbytes = fprintf(fid,'\n\n %1s\n',headstr1) ;
cntbytes = fprintf(fid,'\n %1s\n',headstr2) ;
cntbytes = fprintf(fid,'%1s\n',headstr3) ;
  outm = [vn mlobd(:,1) mavg(:,1) mupbd(:,1) mapprox(:,1)] ;
cntbytes = fprintf(fid,'%6.0d %12.4f %12.4f %12.4f %12.4f\n',outm') ;

  headstr2 = '                     alpha = .1' ;
cntbytes = fprintf(fid,'\n %1s\n',headstr2) ;
cntbytes = fprintf(fid,'%1s\n',headstr3) ;
  outm = [vn mlobd(:,2) mavg(:,2) mupbd(:,2) mapprox(:,2)] ;
cntbytes = fprintf(fid,'%6.0d %12.4f %12.4f %12.4f %12.4f\n',outm') ;

  headstr2 = '                     alpha = .5' ;
cntbytes = fprintf(fid,'\n %1s\n',headstr2) ;
cntbytes = fprintf(fid,'%1s\n',headstr3) ;
  outm = [vn mlobd(:,3) mavg(:,3) mupbd(:,3) mapprox(:,3)] ;
cntbytes = fprintf(fid,'%6.0d %12.4f %12.4f %12.4f %12.4f\n',outm') ;

  headstr2 = '                     alpha = .9' ;
cntbytes = fprintf(fid,'\n %1s\n',headstr2) ;
cntbytes = fprintf(fid,'%1s\n',headstr3) ;
  outm = [vn mlobd(:,4) mavg(:,4) mupbd(:,4) mapprox(:,4)] ;
cntbytes = fprintf(fid,'%6.0d %12.4f %12.4f %12.4f %12.4f\n',outm') ;



fclose('all')
                
