Microsoft Word - Article 4 - 77-747-2-GA.doc ACTA IMEKO December 2013, Volume 2, Number 2, 13 – 19 www.imeko.org ACTA IMEKO | www.imeko.org December 2013 | Volume 2 | Number 2 | 13 Optimization of checking interval and adjustment limit for random-walking process including effects of adjustment uncertainty and time lag Katsuhiro Shirono, Hideyuki Tanaka, Masanori Shiro, Kensei Ehara National Metrology Institute of Japan (NMIJ), National Institute of Advanced Industrial Science and Technology (AIST), Tsukuba 3058565, Japan Section: RESEARCH PAPER Keywords: Checking Interval; Adjustment Limit; Calibration Interval; Random Walk Citation: Katsuhiro Shirono, Hideyuki Tanaka, Masanori Shiro, Kensei Ehara, Optimization of checking interval and adjustment limit for random-walking process including effects of adjustment uncertainty and time lag, Acta IMEKO, vol. 2, no. 2, article 4, December 2013, identifier: IMEKO-ACTA-02 (2013)-02-04 Editor: Paolo Carbone, University of Perugia Received February 15th, 2013; In final form September 12th, 2013; Published December 2013 Copyright: © 2013 IMEKO. This is an open-access article distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited Corresponding author: Katsuhiro Shirono, e-mail: k.shirono@aist.go.jp 1. INTRODUCTION Control charts are useful for monitoring a process in various types of property control, such as a property of products, bias of a measuring device, and so on. In designing a control chart, optimization of the control parameters in terms of economic performance has been discussed so far as reviewed by Keats et al. [1]. When examining optimization, it is necessary to separate static processes, in which the property is time invariant, and dynamic processes, in which the property is time variant. Taguchi et al. [2] proposed an online quality control procedure for dynamic processes by which the checking interval and adjustment limit are economically optimized. (The checking interval corresponds to the calibration interval for the control of a measuring device. The calibration does not contain the adjustment of the measuring device in this case.) However, the procedure is approximately derived, so that the optimization is not precisely given even when the property is in a random- walking process, which is the simplest stochastic process. Adams and Woodall [3] confirmed the property of Taguchi’s procedure and proposed a more accurate optimization procedure for random-walking processes. Srivastava and Wu subsequently reported a simpler solution than that of Adams and Woodall [4], followed by an improved version [5]. Box and Kramer [6] reported a procedure applicable to a more complex type of stochastic process. However, there has been little discussion of the adjustment uncertainty and time lag in the previous reports. With regard to the adjustment uncertainty, when a process is adjusted by replacing some parts, it is uncertain whether the idealistic property is realized. Environmental or human factors can be sources of an adjustment uncertainty. The time lag refers to the difference between the time at which an item is produced or employed and the time at which the item is checked. This was dealt with in Taguchi’s procedure in Taguchi et al. [2]. On the other hand, no time lag was discussed in the above studies, so the effect is unclear even in a random-walking process. Of course, there is some degree of time lag in a realistic process. Clarifying the effects of these two parameters on optimization is, therefore, a prerequisite in order to apply the optimization procedure to a realistic process. ABSTRACT Taguchi’s online quality control aims at the optimization of the checking interval and adjustment limit of a process from the economic point of view. (The checking interval corresponds to the calibration interval for the control of a measuring device.) The present study provides a mathematical formalism and an efficient computation method to take adjustment uncertainty and time lag into consideration in the optimization of a random-walking property. In cases where either the checking cost or the adjustment cost is free, the approximate analytical solution is given. Qualitatively, a larger uncertainty of adjustment yields a shorter optimized checking interval and a larger adjustment limit. On the other hand, a larger time lag yields a smaller checking interval and adjustment limit. Information of approximate equations and a data set useful for efficient computation is provided. ACTA IMEKO | www.imeko.org December 2013 | Volume 2 | Number 2 | 14 In the present study, the effects of the adjustment uncertainty and time lag on optimization are expressed mathematically and an algorithm to optimize the checking interval and adjustment limit is given. Moreover, a case in which the cost of checking is free, which seems interesting in terms of reality, as well as a case in which the cost of adjustment is free, are investigated here. In actuality, analytical solutions can be given for these cases. In the present study, only simple random-walking processes are discussed. Processes with a shift due to assignable causes, e.g., Lorenzen and Vance [7], Banerjee and Rahim [8], Nayebpour and Woodall [9], Surtihadi and Raghavachar [10], and Wu et al. [11], are not within the scope of our interest here. Measurement uncertainty in checking is also ignored. When measurement errors are constant, the effect of a measurement uncertainty can be assessed as in Taguchi’s procedure described in Taguchi et al. [12]. Although CUSUM or EWMA, which are presented in Zhang [13], Reynolds and Stoumbos [14, 15], and other reports, might be useful to attenuate random measurement errors, the application of these techniques raises complexity in the optimization procedure. The discussion of processes with measurement uncertainty will be a future task. This paper is organized as follows: In Section 2, the symbols employed in this paper are listed. Section 3 gives the mathematical formalism of the optimization problem including the effects of the adjustment uncertainty and time lag. Sections 4 and 5 show the numerical procedure for the optimization and its application, respectively. Section 6 provides analytical solutions for cases in which the checking cost or the adjustment cost is free. The contents of the paper are briefly summarized in Section 7. Appendix 1 gives an explanation of the derivation of the equation in the main text. Information of approximate equations and a data set useful for efficient computation is provided in Appendix 2. Appendix 3 explains the derivation of the equation as well as Appendix 1. A data set useful for efficient computation is also provided as an Online Appendix. 2. SYMBOLS t: Time elapsed since the last adjustment. x = xt: Difference from the targeted value at time t. T: Checking interval. D: Adjustment limit. i: Number of checkings since the last adjustment. e: Error in an adjustment (x0 = e). : Population standard deviation of e. t: Time lag, meaning the time lapse between t = kT and the adjustment time. B: Normalized Brownian motion ((dB)2 = dt). : Diffusion coefficient. E[.]: Operator for expectation with respect to both B and e. k: Number satisfying |xkT| > D. Cq: Loss (= cost of poor quality) during t = 0 ~ kT. Cl: Loss during time lag (t = kT ~ kT + t). Cc: Cost of a single checking. Ca: Cost of a single adjustment. CT: Total cost during t = 0 ~ kT+t. L: Expected cost per unit time. T*: Optimum checking interval. D*: Optimum adjustment limit. L*: Optimum expected cost per unit time. : Dimensionless adjustment limit equal to D2/2T. : Dimensionless standard deviation of adjustment error equal to /D. f(, ): Dimensionless expectation of time between adjustments equal to 2/D2·E[kT]. g(, ): Dimensionless expectation of integration of xt2 with t during t = 0 ~ kT equal to 2/D·E[∫ kT0 xt2dt] . h(, ): Dimensionless expectation of square of measured value in the last checking equal to E[xkT2]/D2. n: Iteration number in computation procedure. N: Maximum iteration number in computation procedure. R0: Random variable derived from the normal distribution with a mean of 0 and a variance of 2. R: Random variable derived from the normal distribution with a mean of 0 and a variance of 1. D0: Tolerance limit. L0: Cost per time when the property is out of the tolerance limit. D: Tentative adjustment limit for the parameter determination. T: Tentative checking interval for the parameter determination. ARL(D, T): Average run length with D and T. p: Approximate conditional probability such that |xiT| ≤  on the condition that |x(i-1)T| ≤  (1 ≤ i ≤ k－1). 3. MATHEMATICAL FORMALISM Since there is an uncertainty in adjustment, the difference from the targeted value, x, at t = 0 is not 0 but e, which is derived from the normal distribution with a mean of 0 and a variance of 2. It is assumed that x is along with a random walk for t > 0. In other words, the variation of x, dx, is given by the following stochastic differential equation: dBdx  . (1) B is a normalized Brownian motion such that (dB)2 = dt. Note that it is not necessary for t to be actual time. In a case where the property varies with the number of produced items, t can be the number of items. The difference, x, is checked at t = iT. When |xiT| ≤ D, the process is not adjusted. After the k-th checking satisfying |xkT| > D, the process is stopped and the difference is adjusted. Although the magnitude of the difference at t = kT, xkT, is found to be larger than D, the process cannot be adjusted at t = kT but at t = kT + t because of the existence of the time lag, t. This is a process realized by a Brownian motion, B. Examples of a process without and with the adjustment error and the time lag are shown in Figure 1. A loss during an interval is expressed by the integration of the loss function proportional to the square of the difference, x2, as applied by Srivastava and Wu [4, 5]. Hence, the loss during t = 0 ~ kT is given as follows:  kT dtxqC 0 2 q . (2) The loss during the time lag, t = kT ~ kT+t, is also given as follows:    tkT kT dtxqC 2l . (3) The total cost during t = 0 ~ kT+t, CT, is given as CT = Cq + Cl＋kCc + Ca. What is minimized in the optimization is the expectation of cost per the expectation of time, L = E[CT]/ E[kT+t]. It should be noted that this expectation operator, E[.], does not mean the expectation only with respect to B but also with respect to e. ACTA IMEKO | www.imeko.org December 2013 | Volume 2 | Number 2 | 15 The expectation of the loss during the time lag, i.e., t = kT ~ kT + t, is given as     2222l 2 t q txqEdtxqECE kT TkT kT          . (4) See Appendix 1 for a detailed derivation of this equation. Thus, L can be expressed as follows:           tkTE C T C kTEt q txqEdtxqE tkTE CE L kT kT            ac2220 2 T 2  . (5) Due to the time lag, the above equation is more complicated than the corresponding equation (22) in Box and Kramer [6]. L is the function of the combination of (D, T), and the optimization means the determination of the combination of (D*, T*) by which L is at its minimum. 4. COMPUTATION D and T are the parameters of L, and it is necessary to compute equation (5) to obtain the optimized combination of (D*, T*). The calculation of the three expectations E[kT], E[∫ kT0 xt2dt], and E[xkT2] is a computationally heavy task and it is not realistic to directly conduct the numerical integration involving these expectations every time D and/or T changes in the optimization procedure. To avoid this, the key is to transform these three expectations as follows:      , 2 2 f D kTE  , (6)    , 2 4 0 2 g D dtxE kT        , (7)    ,22 hDxE kT  , (8) where  = D2/2T and  = /D are dimensionless numbers considered as the dimensionless adjustment limit and the dimensionless standard deviation of the adjustment error, respectively. Once f, g, and h are computed for the various combinations of (, ) and references are prepared with tables or approximate equations, the burden of computation of the expectations is considerably reduced by avoiding direct numerical integration. In the present study, a computation procedure entailing the following five steps was conducted to obtain the values of f, g, and h with N = 1000: 1) Set  and . T = 1/. n = 0. 2) n = n + 1. Derive R0 from the normal distribution with a mean of 0 and a variance of 2 where x0 = R0. 3) Repeat to derive R from the normal distribution with a mean of 0 and a variance of 1 where xjT = x(j-1)T + √T·R until |xjT| > 1 4) k = j. Fn = kT, Gn = ∑ 1 0   k i {(T/3)·(xiT2+ x(i+1)T2+ xiT·x(i+1)T) +2T2/6} and Hn = xkT2. Go to 2) until n = N. 5) f(x,) = (1/N)·∑ 10   k i Fn, g(x,) = (1/N)·∑ 1 0   k i Gn, and h(x,) = (1/N)·∑ 10   k i Hn. See Srivastava and Wu [4] (the third equation on p. 2152) for the calculation of Gn in Step 4). Although in the case of  = 0, approximate equations for f and g (or those with minor mathematical differences) have been given in previous studies such as Adams and Woodall [3], Srivastava and Wu [4, 5], Box and Kramer [6], and Box and Jenkins [16], f and g with  = 0 are also recalculated in this study. The approximate equations for f, g and h are provided in Appendix 2. The approximate equations are applicable for 0.1 <  < 10 and 0 <  < 1. Moreover, since the utilization of tables and second-order interpolation seems to be helpful in obtaining more precise values for f, g, and h than approximate equations with a small number of parameters, tables were made for the range of 0.001 <  < 1000 and 0 <  < 1 with smoothing of the computed results. The information on the reference tables for log(f), log(g), and log(h) are also provided in Appendix 2. Once the (a) (b) Figure 1. Examples of a random-walking process (a) without and (b) with the adjustment error, e, and time lag, t. The open circles denote the starting points and the filled circles denote the checking points. ACTA IMEKO | www.imeko.org December 2013 | Volume 2 | Number 2 | 16 values of f, g, and h are obtained, L can be derived from equations (5)–(8) as follows:         tfD C T C f D tqthqDg D q L                       , , 2 1 ,, 2 2 a c 2 2 222 2 4 . (9) Based on this equation, we can obtain (D*, T*) satisfying both ∂L/∂D|D = D* = 0 and ∂L/∂T|T= T* = 0 by an appropriate numerical method, such as Newton's method. 5. EXEMPLIFICATION 5.1. Calculation conditions The example in subsection 23.3 in Taguchi et al. [12] and modified cases are employed in the exemplification in the present study. Although this example is basically for the property control of products, the mathematical formalism is the same as that in the property control of a measurement device. It should be noted that the symbols in this study are different from those in Taguchi et al. [12]. Since the parameters of q and  are not employed in the original example, it seems preferable to touch on the derivation method of them here. According to Taguchi et al. [12], the parameter q can be determined as follows: 2 0 0 D L q  , (10) where D0 is the tolerance limit for the property, and L0 is the cost per time when the property is out of the tolerance limit. Although,  = 0 is the prerequisite in the parameter determination of  in Taguchi et al. [12], the determination method can be modified for the case of  > 0 as follows:     TDARL D , 22 2  , (11) where D and T is the tentative adjustment limit and checking interval for the parameter determination respectively, and ARL(D, T) is the average run length with D and T. However, this parameter determination is appropriate only when T << D2/2. To determine , it might be generally applicable to employ the relationship between the variation in time, (ti－tj), and the variation in difference, (xi－xj), that is E[(xi－xj)2] =  |ti－tj|. The calculation conditions are summarized in Table 1. Case 1 is the same example as that in subsection 23.3 in Taguchi et al. [12], while Cases 2 and 3 are examples emphasizing the adjustment uncertainty and time lag, respectively. The minimizations of equation (9) for all the cases are conducted with f, g and h obtained from the reference tables and second- order interpolation. The information on the reference tables for log(f), log(g), and log(h) are provided in appendix 2. 5.2. Calculation results The calculation results are also summarized in Table 1. The optimized combination for Case 1 is (D*, T*) = (2.98 m, 288 units). This differs from the optimized combination of (3.80 m, 201 units) by Taguchi’s procedure in Taguchi et al. [12]. This result shows the necessity of taking the random- walking behaviour of the property into consideration for accurate optimization. On the other hand, Adams and Woodall’s procedure using equation (4.2) in Adams and Woodall [3] gives the optimized combination of (2.91 m, 275 units). This implies that our procedure is basically the same as that of Adams and Woodall, because the slight gap depends on the marginal difference in the mathematical expression and the accuracy of the approximate equation. Comparing the results of Case 2 with Case 1, the effect of the adjustment uncertainty can be confirmed. Qualitatively, the larger the adjustment uncertainty is, the larger the optimized adjustment limit is and the smaller the optimized checking interval is. It is interesting that the larger adjustment uncertainty derives a larger adjustment limit, because the larger adjustment limit basically means looser control. The optimized cost per time, L*, becomes larger with a larger adjustment uncertainty, so it can be said that the adjustment uncertainty has a negative Table 1. Calculation conditions and calculation results. Case 1 Case 2 Case 3 Calculation conditions  /m·unit-1/2 0.144 q /$·m-2 0.003556 Cc /$ 1.5 Ca /$ 12  /m 0 1.0 0 t /unit 1 1 50 Calculation results D* /m 2.98 3.14 2.85 T* /unit 288 278 281 L* /$·unit-1 0.0342 0.0356 0.0361 ACTA IMEKO | www.imeko.org December 2013 | Volume 2 | Number 2 | 17 impact in terms of the long-term cost. Comparing the results of Case 3 with Case 1, the effect of the time lag can be confirmed. Qualitatively, the larger the time lag is, the smaller both the optimized adjustment limit and checking interval are. The trend of the optimized adjustment limit differs from that of Case 2. In terms of the long-term cost, the time lag has a negative impact as does the adjustment uncertainty. Figure 2 shows the contours of the cost per time, L, for Cases 1 and 2. The qualitative features in both cases are basically the same. It can be seen that the gradient of the variation of L is small in a certain direction. In other words, when making D larger and T smaller or D smaller and T larger from (D*, T*), the variation of cost is relatively insignificant. When certain conditions are imposed on the combination of (D, T), the utilization of a contour might be helpful to determine the optimized combination. 6. ANALYTICAL SOLUTIONS FOR SPECIFIC CASES 6.1. Case of Cc = 0 When the checking cost is essentially zero (or simply pushed aside), the difference, x, will be constantly monitored. This realistically means the total inspection of products or the everyday checking of a measuring device. Although this case has been discussed in previous studies such as Adams and Woodall [3], Srivastava and Wu [4, 5], and Box and Kramer [6], the effects of the adjustment uncertainty and the time lag have not yet been clarified. On the condition of  << D, calculation of two of the expectations, E[kT] and E[∫ kT0 xt2dt], can be replaced by solving the corresponding differential equations as shown in Øksendal [17]. The solutions of the corresponding differential equations are easily obtained: E[kT] = (D2－2)/2 and E[∫ kT0 xt2dt] = (D4 －34)/62. Since E[xkT2] = D2 is obvious, L in equation (5) can be given as follows:   tD CtqtqDDq L    222 a 2242244 2 1 3 6 1   . (12) On the additional condition of t << E[kT] and t << D2/2, D* satisfying ∂L/∂D|D = D* = 0 can be approximately given as follows: t q C D a          22 212 * 6   . (13) This is the exact solution in the case of  = 0 and t = 0. Moreover, it is worth noting that in the case of  = 0 and t = 0, the distribution of x is an isosceles triangle with a range from x = － D to x = D. As suggested in the numerical exemplification in Section 5, equation (13) also implies that the adjustment uncertainty makes the optimized adjustment limit large. 6.2. Case of Ca = 0 According to Miyagawa [18], when |xiT| ≤ , the adjustment makes the variance of x larger. This gives rise to the phenomenon called hunting. Adjusting the process after every instance of checking is not the best means of reducing the variance and the cost per time, even when the adjustment cost is zero or pushed aside. To prevent hunting, the process should be adjusted only after the checking when |xiT| > . Thus, when the adjustment cost is zero or pushed aside, the optimized adjustment limit, D*, is not 0 but . Concerning the optimized checking interval, T*, although Srivastava and Wu [4] reported an analytical equation for this case when  = 0 and t = 0, the effects of the adjustment uncertainty and time lag have not yet been clarified. When 2 << 2T, the conditional probability such that |xiT| ≤  (i = 1, 2, …, k－1) on the condition that |x(i-1)T| ≤  is considered almost constant irrespective of i. If the conditional probability is set as p (p << 1), L in equation (5) can be approximated as follows: (a) (b) Figure 2. Contours of the cost per time, L, with the calculation conditions of (a) Case 1 and (b) Case 2 in Table 1. ACTA IMEKO | www.imeko.org December 2013 | Volume 2 | Number 2 | 18     t p T p C t q p tTq p Tq L          1 12112 c22 222   . (14) See Appendix 3 for a detailed derivation of this equation. When t << T, it is easily found that the value of p is not significant to L because p << 1. On the condition that 2 << 2T and t << T, T* satisfying ∂L/∂T|T= T* = 0 can be derived approximately as follows: t q C T  2 c* 2  . (15) It might be unexpected that  does not have an effect on the determination of T* according to equation (15). This is the exact solution in the case of  = 0 and t = 0. 7. SUMMARY This paper offers a theoretical formalism and an efficient numerical computation algorithm for taking the adjustment uncertainty and time lag into consideration in an online quality control procedure with a random-walking process model. Although the adjustment uncertainty and time lag have been disregarded in most of the previous studies, these two parameters can have an impact on a realistic process. Qualitatively, it can be said that a larger adjustment uncertainty derives a larger optimized adjustment limit and a smaller optimized checking interval, while a larger time lag derives a smaller optimized adjustment limit and checking interval. REFERENCES [1] J.B. Keats, E. Castillo, E.V. Collani, E.M. Saniga, Economic Modeling for Statistical Process Control, J. Qual. Tech. 29 (1997) pp.144–147. [2] G. Taguchi, E.A. Elsayed, T. Hsiang, Quality Engineering in Production Systems, McGraw-Hill, New York, 1989, ISBN 0- 07-100358-4. [3] B.M. Adams, W.H. Woodall, An analysis of Taguchi’s on-line process-control procedure under a random-walk model, Technometrics 31 (1989) pp.401–413. [4] M.S. Srivastava, Y. Wu, A second order approximation to Taguchi’s on-line control procedure, Commun. Statist. － Theory Meth. 20 (1991) pp.2149–2168. [5] M.S. Srivastava, Y. Wu, An improved version of Taguchi’s on- line control procedure, J. Stat. Plan. Infer. 43 (1995) pp.133-145. [6] G. Box, T. Kramer, Statistical process monitoring and feedback adjustment—A discussion, Technometrics 34 (1992) pp.251-267. [7] T.J. Lorenzen, L.C. Vance, The economic design of control charts: A unified approach, Technometrics 28 (1986) pp. 3–10. [8] P.K. Banerjee, M.A. Rahim, Economic design of x -control charts under Weibull shock models, Technometrics 30 (1988) pp. 407-414. [9] M.R. Nayebpour, W.H. Woodall, An analysis of Taguchi's on- line quality monitoring procedures for attributes, Technometrics, 35 (1993) pp. 53–60. [10] J. Surtihadi, M. Raghavachar, Exact economic design of X charts for general time in-control distributions, Int. J. Prod. Res. 32 (1994) pp. 2287-2302. [11] Z. Wu, M. Shamsuzzaman, E. S. Pan, Optimization design of control charts based on Taguchi’s loss function and random process shifts, Int. J. Prod. Res. 42 (2004) pp. 379–390. [12] G. Taguchi, S. Chowdhury, Y. Wu, Taguchi’s Quality Engineering Handbook, John Wiley & Sons, Inc, New Jersey, 2004, ISBN 0-47-141334-8. [13] N.F. Zhang, A statistical control chart for stationary process data, Technometrics 40 (1998) pp. 24–38. [14] M.R. Reynolds Jr., Z.G. Stoumbos, Should exponentially weighted moving average and cumulative sum charts be used with Shewhart limits?, Technometrics 47 (2005) pp. 409–424. [15] M.R. Reynolds Jr., Z.G. Stoumbos, Comparisons of some exponentially weighted moving average control charts for monitoring the process mean and variance, Technometrics 48 (2006) pp. 550–567. [16] G.E.P. Box, G.M. Jenkins, “Further contributions to adaptive quality control: Simultaneous estimation of dynamics: Non-zero costs”, in: Bulletin of the International Statistical Institute, Proceedings of the 34th Session Ottawa 1963, H. Segal (editor), University of Toronto Press, Toronto, 1964, pp. 943–974. [17] B.K. Øksendal, Stochastic differential equations: An introduction with applications (in Japanese), Springer Japan, Tokyo, 1999, ISBN 4-43-170804-9. [18] M. Miyagawa, Technology for getting quality–What the Taguchi method has brought us (in Japanese). Union of Japanese Scientists and Engineers (JUSE), Tokyo, 2006, ISBN 4-81-710339-6. APPENDIX 1: DERIVATION OF EQUATION (4) The expectation of Cl, E[Cl], can be expressed as follows:     2 l 2 , kT T kT kT T kT t kTkT E C E q x dt qE q x y dt              (A1) where yt is the stochastic process satisfying y0 = 0 and dyt = dB. Equation (A1) can be decomposed into two terms as follows:   2 2l 2 2 0 , kT T kT T kT t kTkT kT T kT t E C qE x dt qE y dt qE x T qE y dt                       (A2) because E[∫ TkTkT  xkT·ytdt] = E[xkT]·E[∫ T 0 ytdt] and both E[xkT] and E[∫ T0 ytdt] are 0. As found from the equation for Gn in Step 4) of Section 4 in the main manuscript,   2 2 2 2 2 0 0 2 2 2 2 2 3 6 . 3 6 2 kT T tkT T T T E y dt T T E y y y y T T T E y                              (A3) Note that y0 = 0 and E[yT2] =  T. Substituting equation (A3) into equation (A2) yields     222l 2 T q TxqECE kT   . (A4) This is equation (4). APPENDIX 2: APPROXIMATE EQUATIONS AND REFERENCE TABLES The approximate equations for f, g and h are provided in Table A1. It should be noted that the approximate equations are applicable only for the range of 0.1 <  < 10 and 0 <  < 1. The reference tables for log(f), log(g), and log(h) for the range of 0.001 <  < 1000 and 0 <  < 1 are provided in the Online appendix at http://acta.imeko.org/index.php/acta-imeko/ article/view/IMEKO-ACTA-02%20%282013%29-02-04/198. In the calculation in the main manuscript, the reference tables and second-order interpolation are employed. ACTA IMEKO | www.imeko.org December 2013 | Volume 2 | Number 2 | 19 APPENDIX 3: DERIVATION OF EQUATION (14) Computation of the three expectations E[kT], E[xkT2], and E[∫ kT0 xt2dt] is the key to the derivation of equation (14). As described in the main manuscript, when 2 << 2T, the probability that |xiT| ≤  (i = 1, 2, …, k－1) on the condition that |x(i-1)T| ≤  is considered to be almost constant irrespective of i. That probability is set as p (p << 1). The probability such that |xiT| ≤  (i = 1, 2, …, k－1) and |xkT| >  is therefore given as pk － 1·(1 － p). Hence, E[kT] can be calculated as follows:     p T ppkTkTE k k      1 1][ 1 1 . (A5) Regarding E[xkT2], the distribution of xkT is almost approximated to the normal distribution with a mean of 0 and a variance of 2T, but it should be kept in mind that the magnitude of xkT is larger than . Let the variance of the value that is derived from the normal distribution with a mean of 0 and a variance of 2T and the magnitude of that is larger than  be p2. It can be said that 2T ≈ p2p + E[xkT2](1－p). Since the first term on the right-hand side is marginal compared with the second term, the following expression can be derived: p T xE kT   1 ][ 2 2  . (A6) As found from the equation for Gn in Step (4) of Section 4 in the main manuscript, E[∫ kT0 xt2dt] is given as the sum of E[T/3·(xiT2+ x(i+1)T2+ xiT·x(i+1)T)+2T2/6] for i = 1, …, k－1. Because xiT2 < 2 (i = 1, …, k－1) and 2 << 2T, E[T/3·(xiT 2+ x(i+1)T2 + xiT·x(i+1)T)] is negligibly smaller than 2T2/6 except when i = k － 1. E[xkT2] cannot be ignored but should be assessed by equation (A6). Thus, E[∫ kT0 xt2dt] can be approximated as follows:       2 2 2 2 0 2 2 2 2 2 3 6 . 3 1 6 2 1 kT kT T T E x dt E x E k E kTT T T p p                       (A7) Substituting equations (A5) to (A7) into equation (5) results in equation (14). Table A1. The approximate equations for f(,), g(,), and h(,). These approximate equations are applicable only for the range 0.1 <  < 10 and 0 <  < 1.                  4 0 4 0 lnexp, i j ji ijff  fij j = 0 j = 1 j = 2 j = 3 j = 4 i = 0 1.028 －2.386×10－2 －1.785×10－1 －8.825×10－3 3.794×10－2 i = 1 －4.683×10－1 －2.215×10－3 －2.583×10－1 1.462×10－1 －2.548×10－2 i = 2 9.580×10－2 5.409×10－3 －8.584×10－2 7.445×10－2 －2.008×10－2 i = 3 －2.487×10－3 －6.615×10－4 1.406×10－2 －1.429×10－2 5.733×10－3 i = 4 －2.071×10－3 －5.527×10－4 6.702×10－3 －7.450×10－3 2.685×10－3                  4 0 4 0 lnexp, i j ji ijgg  gij j = 0 j = 1 j = 2 j = 3 j = 4 i = 0 6.610×10－1 －1.503×10－2 2.771×10－1 5.225×10－2 －2.334×10－2 i = 1 －1.135 －2.939×10－4 －1.889×10－1 1.703×10－1 －3.026×10－2 i = 2 2.050×10－1 4.571×10－3 －1.114×10－1 1.167×10－2 2.088×10－2 i = 3 1.051×10－3 1.463×10－3 1.341×10－2 －2.623×10－2 1.073×10－2 i = 4 －4.765×10－3 1.966×10－4 7.269×10－3 －3.900×10－3 1.812×10－4                  4 0 4 0 lnexp, i j ji ijhh  hij j = 0 j = 1 j = 2 j = 3 j = 4 i = 0 1.028 －1.308×10－2 1.257×10－1 9.406×10－2 －2.560×10－2 i = 1 －4.678×10－1 －1.008×10－2 －8.802×10－2 1.961×10－1 －5.909×10－2 i = 2 9.608×10－2 －8.025×10－3 －2.901×10－2 3.029×10－2 －2.293×10－3 i = 3 －2.540×10－3 6.767×10－4 4.578×10－3 －7.100×10－3 2.066×10－3 i = 4 －2.113×10－3 1.066×10－3 －1.137×10－4 9.414×10－4 －1.133×10－3