'' Interactive version that displays only Sqrt[B*W(L,B)]. ''*------------------------------------------------------------------------------* ''|------------------------------------------------------------------------------| ''|www.sim | ''| | ''| | ''| L A BBBB A TTTTT CCC H H 222 | ''| L A A B B A A T C C H H 2 2 | ''| L A A B B A A T C H H 2 | ''| L AAAAA BBBB AAAAA T C HHHHH 2 | ''| L A A B B A A T C H H * 2 | ''| L A A B B A A T C C H H * * 2 | ''| LLLLL A A BBBB A A T CCC H H * 22222 | ''| | ''| | ''| | ''| (Revision of LABATCH Version 1.0) | ''| | ''| SIMSCRIPT II.5 Implementation | ''| | ''| | ''| | ''| G. S. Fishman | ''| | ''| October, 1997 | ''| | ''| | ''| Department of Operations Research | ''| 210 Smith Building CB # 3180 | ''| University of North Carolina at Chapel Hill | ''| Chapel Hill, NC 27599-3180 | ''| (919) 962-8401 | ''| email: gfish@fish.or.unc.edu | ''| fax: (919) 962-0391 | ''| | ''| | ''|------------------------------------------------------------------------------| ''|------------------------------------------------------------------------------| ''| | ''| LABATCH.2 is described in: | ''| | ''| Fishman, G.S. (1997). LABATCH.2: An enhanced implementation of the batch | ''| means method, Technical Report No. 97/04,Operations Research Department, | ''| University of North Carolina at Chapel Hill. | ''| | ''| The original LABATCH analysis package is based on: | ''| | ''| Yarberry, L.S. (1993). Incorporating a dynamic batch size selection mechanism| ''| in a fixed-sample-size batch means procedure. Ph.D. Dissertation, Department | ''| of Operations Research, University of North Carolina at Chapel Hill. | ''| | ''| | ''| Use of Version 1.0 is illustrated in: | ''| | ''| Fishman, G.S. (1996). Monte Carlo: Concepts, Algorithms, and Applications, | ''| Springer-Verlag, New York. | ''| | ''| and | ''| | ''| Fishman, G.S. and L.S. Yarberry (1997). An implementation of the batch means | ''| method, INFORMS Journal on Computing, 9, 296-310. | ''| | ''|------------------------------------------------------------------------------| ''*------------------------------------------------------------------------------* '' '' ''------------------------------------------------------------------ ''------------------------------------------------------------------ '' Preamble ''------------------------------------------------------------------ ''------------------------------------------------------------------ preamble '' normally mode is undefined '' '' Function Definitions '' define PHI as double function define PHI_QUANTILE as double function define STUDENT_T_QUANTILE as double function '' ''------------------------------------------------------------------ '' PSI - the sampling function ''------------------------------------------------------------------ '' define PSI as double function '' end ''------------------------------------------------------------------ ''------------------------------------------------------------------ '' Main ''------------------------------------------------------------------ ''------------------------------------------------------------------ main '' '' '' This program serves 3 purposes: '' '' 1. If IN_UNIT = 0 and OUT_UNIT > 0, '' it runs the M/M/1 simulation, callS BATCH_MEANS '' for each DATA vector, does analysis, and writes '' results to unit OUT_UNIT. '' '' 2. If IN_UNIT > 0 and OUT_UNIT = 0, '' it runs the M/M/1 simulation and writes each DATA vector '' to unit IN_UNIT. '' '' 3. If IN_UNIT > 0 and OUT_UNIT > 0, '' it calls BATCH_MEANS once which then reads DATA '' vectors one at a time from unit IN_UNIT, does '' analysis, and writes results to unit OUT_UNIT. '' '' It uses the PSI function to simulate waiting times in the M/M/1 '' queueing model. ''------------------------------------------------------------------ '' Read values for input parameters. Definitions are the same as '' in routine BATCH_MEANS. '' SEED := initial pseudorandom number seed ''------------------------------------------------------------------ '' define BETA as double variable define DATA as 1-dim double array define DELTA as double variable define I as integer variable define IN_UNIT as integer variable define L_UPPER as integer variable define OUT_UNIT as integer variable define SCREEN as integer variable define T as integer variable define RULE as integer variable define SEED as integer variable '' open unit 11 for input read IN_UNIT,OUT_UNIT,T,DELTA,RULE,BETA,L_UPPER,SCREEN,SEED using unit 11 close unit 11 seed.v(1) = SEED reserve DATA(*) as 2 '' if IN_UNIT > 0 if OUT_UNIT > 0 open unit IN_UNIT for input open unit OUT_UNIT for output call BATCH_MEANS(IN_UNIT,OUT_UNIT,T,2,DATA(*),DELTA,RULE, BETA,L_UPPER,SCREEN) else open unit IN_UNIT for output for I = 1 to T do DATA(2) = 0. DATA(1) = PSI if DATA(1) > 0 DATA(2) = 1 endif write DATA(1),DATA(2) as 2 D(15,6) using unit IN_UNIT loop endif else open unit OUT_UNIT for output for I = 1 to T do DATA(2) = 0. DATA(1) = PSI if DATA(1) > 0 DATA(2) = 1 endif call BATCH_MEANS(IN_UNIT,OUT_UNIT,T,2,DATA(*),DELTA,RULE, BETA,L_UPPER,SCREEN) loop endif '' if OUT_UNIT > 0 close unit OUT_UNIT endif if IN_UNIT > 0 close unit IN_UNIT endif end ''------------------------------------------------------------------ ''------------------------------------------------------------------ '' PSI function for the M/M/1 Queueing Example ''------------------------------------------------------------------ ''------------------------------------------------------------------ function PSI '' ''------------------------------------------------------------------ '' This function determines waiting times in queue for an M/M/1 '' queueing system. The first waiting time is produced from the '' steady-state distribution: With probability 1 - TAU the waiting '' time is 0 and with probability TAU the waiting time is '' exponential with rate nu - zeta. Subsequent waiting times are '' generated via Lindley's recursion: W(I+1) = max(0,W(I)+S(I)-T(I)), '' where S(I) is exponentially distributed with rate nu and T(I) is '' exponentially distributed with rate zeta. We assume a unit '' service rate; that is, nu = 1. This means that the arrival rate, '' zeta, and the traffic intensity, TAU, are identical. ''------------------------------------------------------------------ '' define RC as double variable define SS as saved integer variable define INTERARRIVAL_TIME as double variable define SERVICE_TIME as double variable define TAU as saved double variable define WAITING_TIME as saved double variable '' if SS = 0 SS = 1 open unit 25 for input read tau using unit 25 close unit 25 ''------------------------------------------------------------------ '' Generate the first waiting time. ''------------------------------------------------------------------ if random.f(1) < TAU WAITING_TIME = exponential.f(1.0/(1.0-TAU),1) else WAITING_TIME = 0.0 endif else ''------------------------------------------------------------------ '' Generate a subsequent waiting time. ''------------------------------------------------------------------ INTERARRIVAL_TIME = exponential.f(1.0/TAU,1) SERVICE_TIME = exponential.f(1.0,1) WAITING_TIME = max.f(0.0, WAITING_TIME + SERVICE_TIME - INTERARRIVAL_TIME) endif rc = WAITING_TIME return with RC end ''------------------------------------------------------------------ ''------------------------------------------------------------------ '' Batch Means '' ''******************************************************************* ''* * ''* To Reduce the potential for error introduced by unauthorized * ''*modifications,it is recommended that this Source Code be obtained* ''*by ftp directly from ftp@or.unc.edu. * ''* * ''******************************************************************* ''------------------------------------------------------------------ ''------------------------------------------------------------------ routine BATCH_MEANS (IN_UNIT,OUT_UNIT,T,S_NUM,PSI_VECTOR,DELTA,RULE, BETA,L_UPPER,SCREEN) '' ''------------------------------------------------------------------ '' '' Input Parameters (All are local to routine BATCH_MEANS.) '' '' IN_UNIT := 0 if data are to be repeatedly '' transferred from a main program '' := 30 if data are to be read one vector '' at a time from input file SIMU30 '' (n.b., 30 is merely an example) '' OUT_UNIT := unit to be used for output '' '' T := total number of observations '' '' S_NUM := number of sample series '' '' '' PSI_VECTOR := vector of S_NUM sample values '' '' DELTA := desired significance level for '' confidence intervals (.01 or .05 suggested) '' '' RULE := 1 if ABATCH rule is used '' := 0 if LBATCH rule is used '' '' BETA := significance level for testing for independence '' ( .10 suggested) '' N.B. If BETA = 1, each review uses the FNB rule. '' = -1, each review uses the SQRT rule. '' '' L_UPPER := upper bound on initial number of batches '' (3 <= L_UPPER <= 100) '' '' SCREEN := 1 if interim review estimates are to be '' displayed on the screen '' := 0 otherwise '' '' '' Preamble of calling program must contain declarations: '' '' define PHI as a double function '' define PHI_QUANTILE as a double function '' define STUDENT_T_QUANTILE as a double function '' '' '' Array Dimensions '' '' S_MAX := maximum number of sample series '' '' V_MAX := maximum number of reviews + 1 '' (default is 60 in program) '' '' '' Important Parameters and Variables '' '' B_1 := initial batch size '' B_2_SQRT := minimal permissible batch size using '' SQRT rule '' '' L_1 := initial number of batches '' L_2_SQRT := minimal permissible number of batches '' using SQRT rule '' '' '' HEAD := 1 to produce headings in the computed '' tableau file (default is 1 in program) '' := 0 to omit headings '' '' T_PRIME - Observations 1,...,T_PRIME+1 are used '' to comput B*W(L,B) on final review '' '' ''------------------------------------------------------------------ '' The following section contains array dimensions, parameter '' declarations, working storage and function definitions for the '' LABATCH.2 package. Do not modify them. ''------------------------------------------------------------------ define B_1 as saved integer variable define B_2_SQRT as saved integer variable define BETA as saved double variable define DELTA as saved double variable define HEAD as saved integer variable define IN_UNIT as saved integer variable define L_1 as saved integer variable define L_2_SQRT as saved integer variable define L_UPPER as saved integer variable define OUT_UNIT as saved integer variable define RULE as saved integer variable define S_MAX as saved integer variable define S_NUM as saved integer variable define SCREEN as saved integer variable define T as saved integer variable define T_PRIME as saved integer variable define V_MAX as saved integer variable '' define BETA_VECTOR as saved 1-dim double array define METHOD_VECTOR as saved 1-dim integer array define PSI_VECTOR as saved 1-dim double array '' ''------------------------------------------------------------------ '' '' Working Storage '' '' In the following definitions, the subscripts I and J are used '' to denote the value of a quantity for series I on review J. '' ROW_VECTOR(I) := number of completed reviews '' B_MATRIX(I,J) := batch size '' L_MATRIX(I,J) := number of batches '' P_MATRIX(I,J) := p-value '' V_MATRIX(I,J) := sample variance of a batch mean '' (W(L,B) in the output tableaus) '' X_BAR_MATRIX(I,J) := sample mean '' '' Statistics for the independent case of series I are maintained '' in column # ROW_VECTOR(I)+1 of B_MATRIX, L_MATRIX, P_MATRIX, '' V_MATRIX and X_BAR_MATRIX. '' ''------------------------------------------------------------------ '' '' define ROW_VECTOR as 1-dim saved integer array define B_MATRIX as 2-dim saved integer array define IND_SUM as 1-dim saved double array define L_MATRIX as 2-dim saved integer array define P_MATRIX as 2-dim saved double array define V_MATRIX as 2-dim saved double array define X_BAR_MATRIX as 2-dim saved double array '' define B as integer variable define B_SQRT as integer variable define I as saved integer variable define II as integer variable define J as integer variable define J_SQRT as integer variable define JJ as saved integer variable define KK as alpha variable define L as integer variable define L_SQRT as integer variable define OLD_ROW as saved integer variable define R as integer variable define R_MAX as saved integer variable define ROW as integer variable define R_SQRT as integer variable define SERIES as integer variable define SWAP_INT as integer variable define SS as saved integer variable define TEST as integer variable define C as double variable define RHO as double variable define S as double variable define SWAP_DBL as double variable define TAU as double variable define X as double variable define Y as double variable define Y_SQRT as double variable define B_VECTOR as saved 1-dim integer array define B_SQRT_VECTOR as saved 1-dim integer array define J_VECTOR as saved 1-dim integer array define J_SQRT_VECTOR as saved 1-dim integer array define L_VECTOR as saved 1-dim integer array define L_SQRT_VECTOR as saved 1-dim integer array define R_SQRT_VECTOR as saved 1-dim integer array define R_VECTOR as saved 1-dim integer array define TEST_VECTOR as saved 1-dim integer array define S_VECTOR as saved 1-dim double array define W_IND_VECTOR as saved 1-dim double array define Y_IND_VECTOR as saved 1-dim double array define Y_SQRT_VECTOR as saved 1-dim double array define Y_VECTOR as saved 1-dim double array define Z_IND_VECTOR as saved 1-dim double array define S_MATRIX as saved 2-dim double array define THETA_MATRIX as saved 2-dim double array define W_MATRIX as saved 2-dim double array define XI_MATRIX as saved 2-dim double array '' if S_MAX = 0 S_MAX = S_NUM V_MAX = 60 reserve BETA_VECTOR as S_MAX reserve METHOD_VECTOR as S_MAX reserve ROW_VECTOR as S_MAX reserve B_MATRIX as S_MAX by V_MAX reserve L_MATRIX as S_MAX by V_MAX reserve P_MATRIX as S_MAX by V_MAX reserve V_MATRIX as S_MAX by V_MAX reserve X_BAR_MATRIX as S_MAX by V_MAX reserve IND_SUM as S_MAX reserve B_VECTOR as S_MAX reserve B_SQRT_VECTOR as S_MAX reserve J_VECTOR as S_MAX reserve J_SQRT_VECTOR as S_MAX reserve L_VECTOR as S_MAX reserve L_SQRT_VECTOR as S_MAX reserve R_SQRT_VECTOR as S_MAX reserve R_VECTOR as S_MAX reserve TEST_VECTOR as S_MAX reserve S_VECTOR as S_MAX reserve W_IND_VECTOR as S_MAX reserve Y_IND_VECTOR as S_MAX reserve Y_SQRT_VECTOR as S_MAX reserve Y_VECTOR as S_MAX reserve Z_IND_VECTOR as S_MAX reserve S_MATRIX as S_MAX by V_MAX reserve THETA_MATRIX as S_MAX by V_MAX reserve W_MATRIX as S_MAX by V_MAX reserve XI_MATRIX as S_MAX by V_MAX '' ''------------------------------------------------------------------ '' HEAD = 1 for SERIES = 1 to S_NUM do METHOD_VECTOR(SERIES) = RULE BETA_VECTOR(SERIES) = BETA loop '' '' Compute B_1,B_2_SQRT,L_1,L_2_SQRT based on T observations. '' call SIZE given T,L_UPPER yielding B_1,B_2_SQRT,L_1,L_2_SQRT,T_PRIME '' T_PRIME = 2**trunc.f(log.e.f(real.f(T/B_1/L_1))/log.e.f(2.))*B_1*L_1 endif '2' if IN_UNIT > 0 read PSI_VECTOR using unit IN_UNIT endif for SERIES = 1 to S_NUM do if SS < S_NUM ''------------------------------------------------------------------ '' Begin initialization. ''------------------------------------------------------------------ if SS = 0 if S_NUM > S_MAX print 1 line thus Maximum number of SERIES exceeded. stop endif I = 1 endif SS = SS + 1 B = B_1 B_SQRT = B_2_SQRT L = L_1 L_SQRT = L_2_SQRT J = 0 Y = 0.0 R = 1 J_SQRT = 0 Y_SQRT = 0.0 R_SQRT = 2 R_MAX = max.f( 2*trunc.f(log.e.f(real.f(T/B))/log.e.f(2.0))+1, 2*trunc.f(log.e.f(real.f(T/B_SQRT))/log.e.f(2.0))+2) if R_MAX+1 > V_MAX print 1 line thus Maximum vector size exceeded. stop endif TEST = 1 S = 0.0 ROW_VECTOR(SERIES) = 0 ''------------------------------------------------------------------ '' End initialization. ''------------------------------------------------------------------ else ''------------------------------------------------------------------ '' Begin restoration. ''------------------------------------------------------------------ B = B_VECTOR(SERIES) L = L_VECTOR(SERIES) S = S_VECTOR(SERIES) B_SQRT = B_SQRT_VECTOR(SERIES) J = J_VECTOR(SERIES) J_SQRT = J_SQRT_VECTOR(SERIES) L_SQRT = L_SQRT_VECTOR(SERIES) R = R_VECTOR(SERIES) R_SQRT = R_SQRT_VECTOR(SERIES) TEST = TEST_VECTOR(SERIES) Y = Y_VECTOR(SERIES) Y_SQRT = Y_SQRT_VECTOR(SERIES) ''------------------------------------------------------------------ '' End restoration. ''------------------------------------------------------------------ endif X = PSI_VECTOR(SERIES) if I < T_PRIME S = S + X endif IND_SUM(SERIES) = IND_SUM(SERIES) + X ''------------------------------------------------------------------ '' Update summary data for independent case. ''------------------------------------------------------------------ Y_IND_VECTOR(SERIES) = Y_IND_VECTOR(SERIES) + X * X if I = 1 Z_IND_VECTOR(SERIES) = 0.5 * Y_IND_VECTOR(SERIES) else Z_IND_VECTOR(SERIES) = Z_IND_VECTOR(SERIES) + X * W_IND_VECTOR(SERIES) endif W_IND_VECTOR(SERIES) = X if mod.f(I,B_SQRT) = 0 ''------------------------------------------------------------------ '' If a batch of size B_SQRT is complete: ''------------------------------------------------------------------ J_SQRT = J_SQRT + 1 call BATCH_UPDATES(J_SQRT,R_SQRT, S, W_MATRIX(SERIES,*), S_MATRIX(SERIES,*), THETA_MATRIX(SERIES,*), XI_MATRIX(SERIES,*), R_MAX) Y_SQRT = Y_SQRT + W_MATRIX(SERIES,R_SQRT) * W_MATRIX(SERIES,R_SQRT) endif if mod.f(I,B) = 0 ''------------------------------------------------------------------ '' If a batch of size B is complete: ''------------------------------------------------------------------ J = J + 1 call BATCH_UPDATES(J, R, S, W_MATRIX(SERIES,*), S_MATRIX(SERIES,*), THETA_MATRIX(SERIES,*), XI_MATRIX(SERIES,*), R_MAX) Y = Y + W_MATRIX(SERIES,R) * W_MATRIX(SERIES,R) if J = L ''------------------------------------------------------------------ '' If the current review is over, compute statistics. ''------------------------------------------------------------------ OLD_ROW = ROW_VECTOR(SERIES) ROW_VECTOR(SERIES) = ROW_VECTOR(SERIES) + 1 ROW = ROW_VECTOR(SERIES) B_MATRIX(SERIES,ROW) = B L_MATRIX(SERIES,ROW) = L X_BAR_MATRIX(SERIES,ROW) = S / real.f(B*L) TAU = (S * S) / real.f(L) RHO = Y - TAU V_MATRIX(SERIES,ROW) = RHO / (real.f(B) * real.f(B) * real.f(L-1)) if RHO > 0.0 C = (-TAU + THETA_MATRIX(SERIES,R) + XI_MATRIX(SERIES,R) + 0.5 * W_MATRIX(SERIES,R) * W_MATRIX(SERIES,R)) / RHO P_MATRIX(SERIES,ROW) = 1.0 - PHI(C/sqrt.f(real.f(L-2) / (real.f(L) * real.f(L) - 1.0))) else P_MATRIX(SERIES,ROW) = 0.0 endif ''------------------------------------------------------------------ '' End compute statistics. ''------------------------------------------------------------------ ''------------------------------------------------------------------ '' If J is odd: ''------------------------------------------------------------------ if mod.f(J,2) = 1 Y = Y - W_MATRIX(SERIES,R) * W_MATRIX(SERIES,R) endif Y = Y + 2.0 * XI_MATRIX(SERIES,R) J = trunc.f(J/2) R = R + 2 if TEST = 1 and P_MATRIX(SERIES,ROW) <= BETA_VECTOR(SERIES) if B > 1 B_SQRT = 2 * B_SQRT ''------------------------------------------------------------------ '' If J_SQRT is odd: ''------------------------------------------------------------------ if mod.f(J_SQRT,2) = 1 Y_SQRT = Y_SQRT - W_MATRIX(SERIES,R_SQRT) * W_MATRIX(SERIES,R_SQRT) endif Y_SQRT = Y_SQRT + 2.0 * XI_MATRIX(SERIES,R_SQRT) J_SQRT = trunc.f(J_SQRT / 2) R_SQRT = R_SQRT + 2 endif ''------------------------------------------------------------------ '' FNB rule ''------------------------------------------------------------------ B = 2 * B else if B = 1 ''------------------------------------------------------------------ '' FNB rule ''------------------------------------------------------------------ B = 2 * B else SWAP_INT = 2 * B ''------------------------------------------------------------------ '' SQRT rule ''------------------------------------------------------------------ B = B_SQRT B_SQRT = SWAP_INT SWAP_INT = 2 * L ''------------------------------------------------------------------ '' SQRT rule ''------------------------------------------------------------------ L = L_SQRT L_SQRT = SWAP_INT SWAP_INT = J J = J_SQRT J_SQRT = SWAP_INT SWAP_DBL = Y Y = Y_SQRT Y_SQRT = SWAP_DBL SWAP_INT = R R = R_SQRT R_SQRT = SWAP_INT endif TEST = METHOD_VECTOR(SERIES) endif ''------------------------------------------------------------------ '' Compute statistics for independent case. ''------------------------------------------------------------------ if I > 2 TAU =(IND_SUM(SERIES) * IND_SUM(SERIES)) / real.f(I) RHO = Y_IND_VECTOR(SERIES) - TAU V_MATRIX(SERIES,ROW_VECTOR(SERIES)+1) = RHO/real.f(I-1) if RHO > 0.0 C = (-TAU + Z_IND_VECTOR(SERIES) + 0.5 * W_IND_VECTOR(SERIES) * W_IND_VECTOR(SERIES))/RHO P_MATRIX(SERIES,ROW_VECTOR(SERIES)+1) = 1.0 - PHI(C/sqrt.f(real.f(I-2)/(real.f(I) * real.f(I) - 1.0))) else P_MATRIX(SERIES,ROW_VECTOR(SERIES)+1)=0.0 endif endif endif endif ''------------------------------------------------------------------ '' Begin storage. ''------------------------------------------------------------------ if I < T B_VECTOR(SERIES) = B B_SQRT_VECTOR(SERIES) = B_SQRT J_VECTOR(SERIES) = J J_SQRT_VECTOR(SERIES) = J_SQRT L_VECTOR(SERIES) = L L_SQRT_VECTOR(SERIES) = L_SQRT R_VECTOR(SERIES) = R R_SQRT_VECTOR(SERIES) = R_SQRT S_VECTOR(SERIES) = S TEST_VECTOR(SERIES) = TEST Y_VECTOR(SERIES) = Y Y_SQRT_VECTOR(SERIES) = Y_SQRT endif ''------------------------------------------------------------------ '' End storage. ''------------------------------------------------------------------ '10' loop if I = T call WRITE_TABLEAU_FILE (DELTA,HEAD,S_NUM,T,T_PRIME,BETA_VECTOR(*),IND_SUM(*), METHOD_VECTOR(*),ROW_VECTOR(*),B_MATRIX(*),L_MATRIX(*), P_MATRIX(*),V_MATRIX(*),X_BAR_MATRIX(*),OUT_UNIT) endif ''------------------------------------------------------------------ '' Display selected interim results on screen. ''------------------------------------------------------------------ if SCREEN = 1 if ROW > OLD_ROW if OLD_ROW < 1 write as /,/,/, "LABATCH.2 INTERIM REVIEW STATISTICAL ANALYSIS",/,/,/ using unit 6 JJ = min.f(7,S_NUM) write JJ as "Estimation of Sqrt[sigma^2_infinity] by ", "Sqrt[B*W(L,B)] for Series 1 Through ",I 2,/, "****************************************", "**************************************",/,/,/, " No. of",/, "Review Obs. " using unit 6 for II = 1 to JJ write II as I 1, S 11 using unit 6 write as /,/ using unit 6 endif write ROW,I as I 3,I 10 using unit 6 for SERIES = 1 to JJ write sqrt.f(B_MATRIX(SERIES,ROW)* V_MATRIX(SERIES,ROW)) as E(12,3) using unit 6 '14' write as /,"continue[y/n]? ",+ using unit 6 read KK using unit 5 if KK ne "y" and KK ne "n" go to 14 else if KK = "n" call WRITE_TABLEAU_FILE (DELTA,HEAD,S_NUM,I,T_PRIME,BETA_VECTOR(*), IND_SUM(*),METHOD_VECTOR(*),ROW_VECTOR(*), B_MATRIX(*),L_MATRIX(*),P_MATRIX(*),V_MATRIX(*), X_BAR_MATRIX(*),OUT_UNIT) stop endif OLD_ROW = ROW endif endif ''------------------------------------------------------------------ '' End screen display for this review. ''------------------------------------------------------------------ I = I + 1 if IN_UNIT > 0 if I <= T go to 2 endif endif return end ''------------------------------------------------------------------ ''------------------------------------------------------------------ '' routine BATCH_UPDATES ''------------------------------------------------------------------ ''------------------------------------------------------------------ routine BATCH_UPDATES(J_DUMMY, R_DUMMY, S, W_VECTOR, S_VECTOR, THETA_VECTOR, XI_VECTOR, R_MAX) '' define J_DUMMY as integer variable define J as integer variable define R_DUMMY as integer variable define R as integer variable define R_MAX as integer variable define S as double variable define SS as saved integer variable define W as double variable define W_VECTOR as 1-dim double array define S_VECTOR as 1-dim double array define THETA_VECTOR as 1-dim double array define XI_VECTOR as 1-dim double array '' '' if SS = 0 SS = 1 reserve W_VECTOR as R_MAX reserve S_VECTOR as R_MAX reserve THETA_VECTOR as R_MAX reserve XI_VECTOR as R_MAX endif J = J_DUMMY R = R_DUMMY ''------------------------------------------------------------------ '' While J is even: ''------------------------------------------------------------------ while mod.f(J,2) = 0 do W = S - S_VECTOR(R) XI_VECTOR(R) = XI_VECTOR(R) + W * W_VECTOR(R) W_VECTOR(R) = W S_VECTOR(R) = S J = J / 2 R = R + 2 loop if J = 1 W = S THETA_VECTOR(R) = 0.5 * W * W XI_VECTOR(R) = 0.0 else W = S - S_VECTOR(R) THETA_VECTOR(R) = THETA_VECTOR(R) + W * W_VECTOR(R) endif W_VECTOR(R) = W S_VECTOR(R) = S return end ''------------------------------------------------------------------ ''------------------------------------------------------------------ '' Routine WRITE_TABLEAU_FILE ''------------------------------------------------------------------ ''------------------------------------------------------------------ Routine WRITE_TABLEAU_FILE (DELTA,HEAD,S_NUM,T,T_PRIME,BETA_VECTOR,IND_SUM,METHOD_VECTOR, ROW_VECTOR,B_MATRIX,L_MATRIX,P_MATRIX,V_MATRIX,X_BAR_MATRIX,OUT_UNIT) '' ''------------------------------------------------------------------ '' '' Array Dimensions '' '' S_MAX := the maximum number of sample SERIES '' V_MAX := the maximum number of reviews + 1 '' (one for the independent case) '' '' Input Parameters '' '' DELTA := desired significance level for '' confidence intervals '' HEAD := 1 to produce headings in the computed '' tableau file (default in program) '' := 0 to omit headings '' S_NUM := number of sample series '' T := number of observations '' T_PRIME := number of observations used for variance '' analysis '' BETA_VECTOR := vector of significance levels for '' testing independence for each sample '' IND_SUM := sum of the values of the observations '' METHOD_VECTOR := Vector determining which rule to use '' set to 1 for ABATCH, 0 for LBATCH '' OUT_UNIT :=- unit used for output '' '' Working Storage '' '' In the following definitions, the subscripts I and J are used '' to denote the value of a quantity for series I on review J. '' '' ROW_VECTOR(I) := number of completed reviews '' B_MATRIX(I,J) := batch size '' L_MATRIX(I,J) := number of batches '' P_MATRIX(I,J) := p-value '' V_MATRIX(I,J) := sample variance of a batch average '' (W(L,B) in the output tableaus) '' X_BAR_MATRIX(I,J) := sample mean '' '' Statistics for the independent case of series I are maintained '' in column # ROW_VECTOR(I)+1 of B_MATRIX, L_MATRIX, P_MATRIX, '' V_MATRIX and X_BAR_MATRIX. '' ''------------------------------------------------------------------ '' define BETA_VECTOR as 1-dim double array define IND_SUM as 1-dim double array define METHOD_VECTOR as 1-dim integer array define ROW_VECTOR as 1-dim integer array define B_MATRIX as 2-dim integer array define L_MATRIX as 2-dim integer array define P_MATRIX as 2-dim double array define V_MATRIX as 2-dim double array define X_BAR_MATRIX as 2-dim double array '' define B as integer variable define DELTA as double variable define H as double variable define HEAD as integer variable define L as integer variable define LINES as integer variable define OUT_UNIT as integer variable define P as double variable define ROW as integer variable define S as integer variable define SERIES as integer variable define S_NUM as integer variable define T as integer variable define T_PRIME as integer variable define V as double variable define X_BAR as double variable '' '' if HEAD = 1 write as /, "*---------------------------------------", "---------------------------------------*",/, "|---------------------------------------", "---------------------------------------|",/, "| ", " |",/, "| ", " |",/, "| ", " |",/, "| L A BBBB A T", "TTTT CCC H H 222 |",/, "| L A A B B A A ", " T C C H H 2 2 |",/, "| L A A B B A A ", " T C H H 2 |",/ using unit OUT_UNIT write as "| L AAAAA BBBB AAAAA ", " T C HHHHH 2 |",/, "| L A A B B A A ", " T C H H * 2 |",/, "| L A A B B A A ", " T C C H H * * 2 |",/, "| LLLLL A A BBBB A A ", " T CCC H H * 22222 |",/, "| ", " |",/ using unit OUT_UNIT write as "| ", " |",/, "| ", " |",/, "| (Revision of LABA", "TCH Version 1.0) |",/, "| ", " |",/, "| SIMSCRIPT II.5 ", "Implementation |",/, "| ", " |",/, "| ", " |",/, "| G. S. Fi", "shman |",/, "| ", " |",/, "| October,", " 1997 |",/, "| ", " |",/, "| ", " |",/, "| ", " |",/ using unit OUT_UNIT write as "| Department of Op", "erations Research |",/, "| 210 Smith Buil", "ding CB # 3180 |",/, "| University of North C", "arolina at Chapel Hill |",/, "| Chapel Hill, ", "NC 27599-3180 |",/, "| (919) 9", "62-8401 |",/, "| email: gfish@f", "ish.or.unc.edu |",/, "| fax: (919) 9", "62-0391 |",/, "| ", " |",/, "| ", " |",/ using unit OUT_UNIT write as "|---------------------------------------", "---------------------------------------|",/, "|---------------------------------------", "---------------------------------------|",/, "| ", " |",/, "| ", " |",/, "| LABATCH.2 is described in: ", " |",/, "| ", " |",/, "| Fishman, G.S. (1997). LABATCH.2: An en", "hanced implementation of the batch |",/, "| means method, Technical Report No. 97/", "04, Operations Research Department, |",/, "| University of North Carolina at Chapel", " Hill. |",/, "| ", " |",/, "| The original LABATCH analysis package ", "is based on: |",/, "| ", " |",/, "| Yarberry, L.S. (1993). Incorporating a", " dynamic batch size selection mechanism|",/, "| in a fixed-sample-size batch means pro", "cedure. Ph.D. Dissertation, Department |",/, "| of Operations Research, University of ", "North Carolina at Chapel Hill. |",/ using unit OUT_UNIT write as "| ", " |",/, "| ", " |",/, "| Use of Version 1.0 is illustrated in: ", " |",/, "| ", " |",/, "| Fishman, G.S. (1996). Monte Carlo: Con", "cepts, Algorithms, and Applications, |",/, "| Springer-Verlag, New York. ", " |",/ using unit OUT_UNIT write as "| ", " |",/, "| and ", " |",/, "| ", " |",/, "| Fishman, G.S. and L.S. Yarberry (1997)", ". An implementation of the batch means |",/, "| method, INFORMS Journal on Computing, ", "9, 296-310. |",/, "| ", " |",/ using unit OUT_UNIT write as "|---------------------------------------", "---------------------------------------|",/, "*---------------------------------------", "---------------------------------------*",/ using unit OUT_UNIT endif '' write T,(1.0-DELTA)*100 as /,/, "Final Tableau ",/, " Mean ", "Estimation ",/, " *****", "********** ",/, " (t =",I 9, " ) ",/,/, " ", S 6,D(4,1), "% ",/, " _ Standard Error C", "onfidence Interval _ ",/, " Series X sqrt[B*W(L,B)/t] L", "ower Upper (Upper-Lower)/|X|",/,/ using unit OUT_UNIT '' for SERIES = 1 to S_NUM do ROW = ROW_VECTOR(SERIES) B = B_MATRIX(SERIES,ROW) L = L_MATRIX(SERIES,ROW) S = L*B if S = T_PRIME S = T endif X_BAR = IND_SUM(SERIES)/real.f(S) V = V_MATRIX(SERIES,ROW) H = STUDENT_T_QUANTILE(1.0-DELTA/2., L-1) * sqrt.f(B*V/S) '' write SERIES,X_BAR,sqrt.f(B*V/S),X_BAR-H,X_BAR+H, 2*H/max.f(1,abs.f(X_BAR)) as I 4, S 4, E(11,3), S 4, E(11,3), S 3, E(11,3), S 2, E(11,3), S 5, E(11,3) ,/,/ using unit OUT_UNIT loop '' write 100 * (L*B/T) as /, S 3, "_ ", " ",/, S 3, "X(t) is based on all t observations. ",/, S 3, "W(L,B) is based on first", D(7,2), "% of the t observations. " using unit OUT_UNIT '' if HEAD = 1 for LINES = 1 to 66 - 11 -2*S_NUM do write as / using unit OUT_UNIT loop endif '' for SERIES = 1 to S_NUM do if HEAD = 1 if METHOD_VECTOR(SERIES) = 1 write SERIES as "Interim Review Tableau ",/, S 23, " ABATCH Data Analysis for Series ",I 2 using unit OUT_UNIT else write SERIES as "Interim Review Tableau ",/, S 23, " LBATCH Data Analysis for Series ",I 2 using unit OUT_UNIT endif '' write (1. - DELTA) * 100.0 as /,/, " ", " ",D(4,1), "% ",/, " _ ", " Confidence Interval ",/, "Review L*B L B X ", " Lower Upper Sqrt[B*W(L,B)] p-value",/,/ using unit OUT_UNIT endif '' for ROW = 1 to ROW_VECTOR(SERIES) do B = B_MATRIX(SERIES,ROW) L = L_MATRIX(SERIES,ROW) P = P_MATRIX(SERIES,ROW) X_BAR = X_BAR_MATRIX(SERIES,ROW) S = L*B if S = T_PRIME S = T X_BAR = IND_SUM(SERIES)/real.f(S) endif V = V_MATRIX(SERIES,ROW) H = STUDENT_T_QUANTILE(1.0-DELTA/2., L-1) * sqrt.f(B*V/S) '' write ROW,L*B, L, B, X_BAR, X_BAR-H, X_BAR+H, sqrt.f(real.f(B)*V), P as I 3, S 1, I 8, S 1, I 8, S 1, I 8, S 1, E(11,3), S 1, E(11,3), S 1, E(11,3), S 1, E(11,3), S 2, D(6,4),/ using unit OUT_UNIT loop '' if HEAD = 1 write as /,S 3,"If data are independent:",/,/ using unit OUT_UNIT endif L = L*B B = 1 P = P_MATRIX(SERIES,ROW_VECTOR(SERIES)+1) V = V_MATRIX(SERIES,ROW_VECTOR(SERIES)+1) H = STUDENT_T_QUANTILE(1.0-DELTA/2., S-1) * sqrt.f(B*V/S) '' write T,T,B,X_BAR,X_BAR-H,X_BAR+H,sqrt.f(real.f(B)*V),P as S 4, I 8, S 1, I 8, S 1, I 8, S 1, E(11,3), S 1, E(11,3), S 1, E(11,3), S 1, E(11,3), S 2, D(6,4),/ using unit OUT_UNIT '' write max.f(0,BETA_VECTOR(SERIES)) as /,S 1,D(5,2), " significance level for independence testing." using unit OUT_UNIT '' write row-1,100 * (L*B/T) as /, S 3, "Review", I 3, " used the first",D(7,2), "% of the t observations for W(L,B). " using unit OUT_UNIT '' if HEAD = 1 for LINES = 1 to 66-13-ROW_VECTOR(SERIES) do write as / using unit OUT_UNIT loop endif loop return end ''------------------------------------------------------------------ ''------------------------------------------------------------------ '' PHI function ''------------------------------------------------------------------ ''------------------------------------------------------------------ function PHI(X) '' ''------------------------------------------------------------------ '' This function computes the value of the normal c.d.f. at the '' point X. '' '' Reference: Abramowitz, M. and I.A. Stegun (1964). Handbook of '' Mathematical Functions with Formulas, Graphs and Mathematical '' Tables. Washington, U.S. Government Printing Office ''------------------------------------------------------------------ '' define X as double variable define RC as double variable define Y as double variable '' if X > 0.0 if X > 45.0 RC = 1.0 else Y = X RC = 1.0-0.50/(1.0+(0.0498673470+(0.0211410061+ (0.0032776263+(0.0000380036+(0.0000488906+ 0.0000053830*Y)*Y)*Y)*Y)*Y)*Y)**16 endif else if X < -45.0 RC = 0.0 else Y = -X RC = 0.50/(1.0+(0.0498673470+(0.0211410061+(0.0032776263 +(0.0000380036+(0.0000488906+0.0000053830*Y)*Y)*Y)*Y) *Y)*Y)**16 endif endif return with RC end ''------------------------------------------------------------------ ''------------------------------------------------------------------ '' PHI_QUANTILE function ''------------------------------------------------------------------ ''------------------------------------------------------------------ function PHI_QUANTILE(BETA) '' ''------------------------------------------------------------------ '' This function computes the quantile of order BETA from the normal '' distribution. '' '' Reference: Abramowitz, M. and I.A. Stegun (1964). Handbook of '' Mathematical Functions with Formulas, Graphs and Mathematical '' Tables. Washington, U.S. Government Printing Office ''------------------------------------------------------------------ '' define BETA as double variable define DELTA as double variable define HIGH as double variable define LOW as double variable define MID as double variable define RC as double variable define W as double variable '' if BETA >= 0.5 DELTA = BETA else DELTA = 1.0-BETA endif W = Sqrt.f(-log.e.f((1.0-DELTA)*(1.0-DELTA))) LOW = W-(2.515517+(0.802853+0.010328*W)*W)/(1.0+(1.432788+(0.189269+ 0.001308*W)*W)*W) HIGH = LOW + .00045 LOW = LOW - .00045 while abs.f(MID - (LOW+HIGH)/2.0) > 0.000000000000001 do MID = (LOW + HIGH)/2.0 if PHI(MID) < DELTA LOW = MID else HIGH = MID endif loop if BETA >= 0.5 RC = MID else RC = -MID endif return with RC end ''------------------------------------------------------------------ ''------------------------------------------------------------------ '' Student t quantile function ''------------------------------------------------------------------ ''------------------------------------------------------------------ function STUDENT_T_QUANTILE(BETA,L) '' ''------------------------------------------------------------------ '' This function computes the quantile of order BETA from the '' Student t distribution with L degrees of freedom. '' '' Reference: Abramowitz, M. and I.A. Stegun (1964). Handbook of '' Mathematical Functions with Formulas, Graphs and Mathematical '' Tables. Washington, U.S. Government Printing Office ''------------------------------------------------------------------ '' define BETA as double variable define L as integer variable '' define DELTA as double variable define DFL as double variable define HIGH as double variable define LOW as double variable define MID as double variable define RC as double variable define TEST as double variable define W as double variable define X as double variable define Y as double variable '' if BETA = 0.5 RC = 0.0 else if L >= 7 X = PHI_QUANTILE(BETA) Y = X*X DFL = real.f(L) RC = X*(1.0+((-945.0+(-3600.0+(2880.0+23040.0*DFL)*DFL)* DFL)+((-1920.0+(4080.0+(15360.0+23040.0*DFL)*DFL)* DFL)+((1482.0+(4560.0+4800.0*DFL)*DFL)+((776.0+720.0* DFL)+79.0*Y)*Y)*Y)*Y)/(DFL**4*92160.0)) else select case L case 1 RC = tan.f((BETA-0.5)*3.1415926535898) case 2 RC = (2.0*BETA-1.0)/sqrt.f(2.0*BETA*(1.0-BETA)) default if BETA >= 0.5 DELTA = BETA else DELTA = 1.0-BETA endif LOW = PHI_QUANTILE(DELTA) HIGH = (2.0*DELTA-1.0)/ sqrt.f(2.0*DELTA*(1.0-DELTA)) while abs.f(MID - (LOW+HIGH)/2.0) > 0.000000000000001 do MID = (LOW + HIGH)/2.0 W = MID * MID select case L case 3 TEST = sqrt.f(3.0)*tan.f(3.1415926535898* (DELTA-0.5)-(MID*sqrt.f(3.0))/(3.0+W)) case 4 TEST = (2.0*DELTA-1.0)* sqrt.f(4.0+W)**3/(6.0+W) case 5 TEST = sqrt.f(5.0)*tan.f(3.1415926535898* (DELTA-0.5)-(MID*sqrt.f(5.0))/(3.0* (5.0+W)**2)*(25.0+3.0*W)) case 6 TEST = (2.0*DELTA-1.0)*sqrt.f(6.0+W)**5/ (W*W+15.0*W+67.5) endselect if TEST > mid LOW = MID else HIGH = MID endif loop if BETA >= 0.5 RC = MID else RC = -MID endif endselect endif endif return with RC end '' ''---------------------------------------------------------------- ''---------------------------------------------------------------- '' routine SIZE finds L_1 and B_1 first to maximize T_PRIME (<=100) '' and then to maximize the number of reviews. ''---------------------------------------------------------------- ''---------------------------------------------------------------- '' routine SIZE given T,L_UPPER yielding B_1,B_2_SQRT,L_1,L_2_SQRT, T_PRIME define ALPHA as integer variable define ALPHA_MAX as integer variable define B and L as saved 1-dim integer arrays define B_1 as integer variable define B_2_SQRT as integer variable define I as integer variable define I_MAX as saved integer variable define L_1 as integer variable define L_2_SQRT as integer variable define L_UPPER as integer variable define T as integer variable define T_PRIME as integer variable define TEMP_1 as integer variable define TEMP_2 as integer variable '' reserve B(*) and L(*) as 46 '' '' Assign values for B() and L(). '' B( 1) = 1 L( 1) = 3 B( 2) = 5 L( 2) = 7 B( 3) = 5 L( 3) = 21 B( 4) = 7 L( 4) = 15 B( 5) = 7 L( 5) = 25 B( 6) = 7 L( 6) = 35 B( 7) = 8 L( 7) = 11 B( 8) = 9 L( 8) = 13 B( 9) = 12 L( 9) = 17 B(10) = 12 L(10) = 51 B(11) = 12 L(11) = 85 B(12) = 15 L(12) = 21 B(13) = 17 L(13) = 36 B(14) = 17 L(14) = 60 B(15) = 17 L(15) = 84 B(16) = 19 L(16) = 27 B(17) = 21 L(17) = 25 B(18) = 21 L(18) = 35 B(19) = 22 L(19) = 31 B(20) = 22 L(20) = 93 B(21) = 26 L(21) = 37 B(22) = 27 L(22) = 57 B(23) = 29 L(23) = 41 B(24) = 31 L(24) = 66 B(25) = 32 L(25) = 45 B(26) = 33 L(26) = 47 B(27) = 36 L(27) = 51 B(28) = 36 L(28) = 68 B(29) = 36 L(29) = 85 B(30) = 39 L(30) = 55 B(31) = 41 L(31) = 87 B(32) = 43 L(32) = 61 B(33) = 46 L(33) = 65 B(34) = 49 L(34) = 69 B(35) = 50 L(35) = 71 B(36) = 51 L(36) = 60 B(37) = 51 L(37) = 84 B(38) = 53 L(38) = 75 B(39) = 56 L(39) = 79 B(40) = 60 L(40) = 85 B(41) = 63 L(41) = 89 B(42) = 66 L(42) = 93 B(43) = 67 L(43) = 95 B(44) = 70 L(44) = 99 B(45) = 84 L(45) = 85 B(46) = 85 L(46) = 96 '' I_MAX = 1 ALPHA_MAX = 0 TEMP_2 = 0 '' for I = 1 to 46 do if T >= B(I)*L(I) and L(I) <= L_UPPER ALPHA = trunc.f(log.e.f(T/(L(I)*B(I)))/log.e.f(2.)) TEMP_1 = 2**ALPHA*B(I)*L(I) if TEMP_1 >= TEMP_2 if TEMP_1 > TEMP_2 or L(I)*B(I) < L(I_MAX)*B(I_MAX) I_MAX = I ALPHA_MAX = ALPHA TEMP_2 = TEMP_1 endif endif endif 'AAA' loop '' B_1 = B(I_MAX) B_2_SQRT = 3 '' if B_1 > 1 B_2_SQRT = trunc.f(sqrt.f(2)*B_1 + .5) endif '' L_1 = L(I_MAX) L_2_SQRT = trunc.f(sqrt.f(2)*L_1 + .5) T_PRIME = TEMP_2 return end