Analysis of inconsistent source sampling in monte carlo weight-window variance reduction methods

Chia sẻ: Thi Thi | Ngày: | Loại File: PDF | Số trang:9

Thêm vào BST

Báo xấu

14
lượt xem 0
download

Download Vui lòng tải xuống để xem tài liệu đầy đủ

This paper develops an original framework that mathematically expresses the coupling of the weight window and source biasing techniques, allowing the authors to explore the impact of inconsistent source sampling on the variance of MC results. A numerical experiment supports this new framework and suggests that certain classes of problems may be relatively insensitive to inconsistent source sampling schemes with moderate levels of splitting and rouletting.

Chủ đề:

Bình luận(0) Đăng nhập để gửi bình luận!

Lưu

Nội dung Text: Analysis of inconsistent source sampling in monte carlo weight-window variance reduction methods

Nuclear Engineering and Technology 49 (2017) 1172e1180 Contents lists available at ScienceDirect Nuclear Engineering and Technology journal homepage: www.elsevier.com/locate/net Original Article Analysis of inconsistent source sampling in monte carlo weight-window variance reduction methods David P. Griesheimer*, Virinder S. Sandhu Naval Nuclear Laboratory, P.O. Box 79, West Mifﬂin, PA 15122-0079, USA a r t i c l e i n f o a b s t r a c t Article history: The application of Monte Carlo (MC) to large-scale ﬁxed-source problems has recently become possible Received 31 May 2017 with new hybrid methods that automate generation of parameters for variance reduction techniques. Accepted 27 July 2017 Two common variance reduction techniques, weight windows and source biasing, have been automated Available online 14 August 2017 and popularized by the consistent adjoint-driven importance sampling (CADIS) method. This method uses the adjoint solution from an inexpensive deterministic calculation to deﬁne a consistent set of Keywords: weight windows and source particles for a subsequent MC calculation. One of the motivations for source Monte Carlo consistency is to avoid the splitting or rouletting of particles at birth, which requires computational Variance Reduction Importance Sampling resources. However, it is not always possible or desirable to implement such consistency, which results in CADIS inconsistent source biasing. This paper develops an original framework that mathematically expresses Weight Windows the coupling of the weight window and source biasing techniques, allowing the authors to explore the impact of inconsistent source sampling on the variance of MC results. A numerical experiment supports this new framework and suggests that certain classes of problems may be relatively insensitive to inconsistent source sampling schemes with moderate levels of splitting and rouletting. © 2017 Korean Nuclear Society, Published by Elsevier Korea LLC. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). 1. Introduction source particles produced with weights that lay outside of the weight window for the corresponding birth state of the particle. Many common Monte Carlo (MC) variance reduction techniques Source particles produced with an inconsistent birth weight are rely on weight windows to control the statistical weight of particles immediately subjected to weight adjustment (splitting or roulette), during transport in order to minimize the variance of ﬂux or re- which is widely believed to decrease the overall effectiveness of the action rate tallies in a speciﬁed region of phase space (e.g., position, weight-window variance reduction scheme. energy, direction). For these techniques, it is well known that the In 1998, Wagner and Haghighat [2] introduced the consistent optimal particle weight at each phase location in a problem is given adjoint-driven importance sampling (CADIS) method for creating by the objective-driven adjoint ﬂux for that location [1]. Particles adjoint-based sets of weight-window parameters based off of a with weight outside of a predeﬁned window about the optimal deterministic estimate for the adjoint ﬂux. In addition, Wagner and weight are subjected to rouletting or splitting (which adjust the Haghighat showed that the deterministic estimate of the adjoint weight in a fair manner) in order to maintain the weight within the ﬂux can also be used to deﬁne a biased source deﬁnition that is weight window. This weight adjustment applies at particle events consistent with the weight-window parameters. Here, consistent where the weight changes (e.g., births, collisions) as well as when means that source particles from the biased source are born with a particles move between regions of phase space with different weight that lies at the center point of the weight window corre- weight-window parameters. sponding to the initial (birth) phase state of the particle. The In early implementations of weight-window variance reduction development of a method for simultaneously creating a consistent methods, inconsistencies between the radiation source distribution source along with the weight-window parameters was a signiﬁcant and the weight window parameters for the simulation resulted in advancement and is a major advantage of the CADIS method. The same consistency is also found in the forward-weighted CADIS (FW-CADIS) method used to calculate global MC solutions [3]. However, even with the advancement of the CADIS method, there * Corresponding author. are some situations where it can be difﬁcult to ensure a completely E-mail address: david.griesheimer@unnpp.gov (D.P. Griesheimer). http://dx.doi.org/10.1016/j.net.2017.07.017 1738-5733/© 2017 Korean Nuclear Society, Published by Elsevier Korea LLC. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/ licenses/by-nc-nd/4.0/). D.P. Griesheimer, V.S. Sandhu / Nuclear Engineering and Technology 49 (2017) 1172e1180 1173 consistent source distribution, and, therefore, the CADIS method ever been performed to our knowledge. Although it appears self- cannot be applied as intended. evident that frequently adjusting the initial weight of source par- For example, the biased source produced by CADIS is based on ticles is counterproductive, it seems reasonable that a small initial an estimate of the adjoint ﬂux distribution produced from a weight adjustment for source particles may be acceptable for many deterministic solution method e typically the discrete ordinates applications. However, such a conclusion requires a thorough (SN) method. As a result, the adjoint ﬂux and the resulting biased characterization of the effect of inconsistent source sampling based source distribution are discretized over space, energy, and direc- on the degree of inconsistency. tion. In order to reduce the amount of time required to generate In this paper, we develop a mathematical framework for quan- weight-window parameters, a relatively coarse discretization may tifying the impact of inconsistent source sampling on the variance be used to estimate the adjoint ﬂux [3]. Although the discretized of tallied quantities in a MC simulation. The derived relationships biased source produced from CADIS is guaranteed to be consistent are supported with results from numerical experiments and pro- with the corresponding weight-window parameters, the CADIS vide a foundation for additional analyses tailored to a variety of source does not preserve higher-order information about the speciﬁc applications. source distribution, which causes discretization error within the sampled MC source. In practice, many MC codes that use CADIS 2. Expected variance by sample scheme simply assume that source particles are uniformly distributed within each discretized “bin” of the CADIS source. However, this In this section, we derive the expected variance in estimated assumption may still lead to a bias in the results from the MC response for several different source sampling schemes. Prior to transport simulation, especially for cases where there is detailed proceeding, it is useful to deﬁne notation and signiﬁcant statistical structure in the true source distribution, such as the energy spec- relationships that will be used throughout the remainder of the trum of a decay source. Any modiﬁcation of the source distribution paper. to reduce this bias may lead to inconsistencies between the adjusted source and the original weight-window parameters. 2.1. Notation and basic relationships In other situations, the radiation source may be “presampled” from a preceding calculation and stored as a census ﬁle containing In MC transport methods, each history can be viewed as the detailed state information about the source particles. This scenario combination of two separate realizations: the initial (birth) state of is common when generating secondary radiations during MC the source particle, denoted x, and the response of the history as transport (e.g., (n,g) or (g,n) reactions), exchanging information measured against some predetermined objective, denoted r. In this between MC eigenvalue and ﬁxed-source calculations, or in SN/MC context, we have assumed that the initial particle state, x, is a vector or MC/MC splice calculations where particles that reach a pre- that includes properties such as the birth energy, position, and deﬁned “trapping surface” are stored for use in a subsequent MC direction of the particle, and that the response, r, is a scalar value. simulation [4]. In these cases, it is straightforward to collapse the Note that these are arbitrary assumptions and may be changed particle census into a discretized source representation for use in without loss of generality. CADIS. However, replacing the census source by the discretized To an external observer, ignorant of the inner workings of the representation would eliminate valuable information stored in the MC transport algorithm, it appears that each history produces a census, such as the correlations between the position, energy, and realization (x,r) from the joint probability density function (PDF) direction of each particle, and is not a practical solution for many p(x,r). Based on the properties of joint probability distributions, it splice calculations. Therefore, retaining the particle census in- follows that the expected value and variance of any function f(x,r) troduces inconsistent source sampling into the subsequent MC applied to a realization of the joint PDF is given by calculation. In addition, for some types of analyses, it is desirable to generate Z∞ Z a single set of weight-window parameters that can be used with a E½f ðx; rÞ ¼ f ðx; rÞpðx; rÞdxdr (1) range of similar model conﬁgurations, often representing source, ∞ G geometry, or composition perturbations with respect to a single reference scenario. In these cases, the CADIS method is well suited and for determining the weight-window parameters and a consistent h i source for the reference conﬁguration, but it can become expensive Var½f ðx; rÞ ¼ E f 2 ðx; rÞ E½f ðx; rÞ2 ; (2) if the weight-window parameters and/or the consistent biased source must be regenerated for every model perturbation. In where G is the domain for the birth state of the particle. practice, a single set of weight-window parameters is often used for Note that the joint PDF p(x,r) can be written as the product of all of the model perturbations, regardless of whether each source conditional and marginal probability distributions, p(x,r) ¼ p(rjx) p(x). distribution is actually consistent with the weight-window In this case, the expected value of the function f(x,r) can be expressed parameters. as Finally, we note that, although the CADIS method has proven to Z be extremely successful, there are still weight-window variance E½f ðx; rÞ ¼ Ex ½Er ½f ðx; rÞjx ¼ Er ½f ðx; rÞjxpðxÞdx; (3) reduction techniques in use [5,6], and under development, that do not produce a consistent biased source distribution, for a variety of G reasons. where For any situation where an inconsistent source distribution may be used with weight-window variance reduction, it is important to Z∞ have a clear understanding of the effects of weight adjustment via Er ½f ðx; rÞjx ¼ f ðx; rÞpðrjxÞdr; (4) splitting or rouletting immediately after particle birth. Although ∞ the conventional wisdom maintains that any weight adjustment at birth will reduce the effectiveness of the weight-window variance and subscripts have been included on the expectation operators to reduction, no systematic, formal investigation of this conjecture has clarify which variable the expectation is taken with respect to. The 1174 D.P. Griesheimer, V.S. Sandhu / Nuclear Engineering and Technology 49 (2017) 1172e1180 " # " # relationship in Eq. (3) is commonly referred to as the fundamental 1 XN 1 XN Var½r property of conditional expected values, the law of total expecta- Var½m ^ r ¼ Var ri ¼ 2 Var ri ¼ : (10) N i¼1 N N tion, or the law of iterated expectation. i¼1 Similarly, the law of total variance can be used to express the Substituting Eq. (9) into Eq. (10) yields the ﬁnal expression variance of f(x,r) in terms of the conditional probability, Ex ½Var½rjx Var½Er ½rjx Var½m ^r ¼ þ : (11) Var½ f ðx; rÞ ¼ Ex ½Var½f ðx; rÞjx þ Var½Er ½f ðx; rÞjx; (5) N N where 2.3. Importance sampling scheme 2 Var½ f ðx; rÞjx ¼ Er ðf ðx; rÞ Er ½f ðx; rÞjxÞ x : (6) The previous section established the basic estimators for the expected response and associated variance using only unbiased Note that Eq. (5) demonstrates that the total variance of f(x,r) realizations from the joint distribution p(x,r). However, it is well includes two components: the ﬁrst term gives the variance due known that the variance of the expected response can be reduced only to the randomness of r for a ﬁxed value of x, referred to as the via importance sampling [7]. In adjoint-driven importance sam- transport variance, and the second term gives the variance due to pling, the underlying source distribution (in this case, p(x)) is the effect of the randomness of x on the conditional distribution altered so that the initial particle state is drawn from a distribution p(rjx), which is referred to as the source variance. The use of the that is proportional to the corresponding response of the source. In conditional probability is often used for analysis of MC radiation order to preserve the original source distribution, the observed transport algorithms because it explicitly shows the dependence of response of the particles is weighted in proportion to the ratio of the response on the initial state of a particle. the probabilities of the initial state in the original and modiﬁed Now that the notation and relevant statistical relationships have probability distributions. been established, we will consider the variance of the estimated For the MC radiation transport process, let us deﬁne a biased response for several different source sampling schemes. source distribution p0 (x) such that 2.2. Unbiased sampling scheme p0 ðx; rÞ ¼ pðrjxÞp0 ðxÞ: (12) Consider a fair MC transport process where each particle history The corresponding weighting factor for each state point is given by i begins with an initial particle state x ~ i sampled from the distri- bution p(x) and produces a corresponding response ~ri according to wðxÞ ¼ pðxÞ=p0 ðxÞ: (13) the probability distribution p(rjx). As described in the previous With importance sampling, the unbiased statistic for estimating section, the initial state and corresponding response can be viewed the expected response is as either sequential (conditional) realizations, or as a single reali- zation of the ordered pair (x ~ i ,~ri ) from the joint PDF p(x,r). 1 XN 0 For a simulation with N independent particle histories, an esti- b imp m r ¼ ~i ~r i ; w x (14) mate for the expected response can be computed with the unbiased N i¼1 sample statistic ~0i denotes a source state realization sampled from the where x 1 XN modiﬁed distribution p0 (x0 ). m br ¼ ~r : (7) N i¼1 i Similarly, the sample variance statistic N 2 The variance of the response can be estimated using the sample 1 X b 2 s r;imp ¼ ~0i ~r i m w x b imp r ; (15) variance statistic N1 i¼1 1 X N is an unbiased estimator for Var[w(x0 )r], where 2 b sr ¼ ð~r m b r Þ2 ; (8) N 1 i¼1 i Var½wðx0 Þr ¼ Ex0 ½Var½wðx0 Þrjx0 þ Var½Er ½wðx0 Þrjx0 : (16) which is an unbiased estimator for the variance of the response r, Recognizing that the weight factor w(x0 ) is a constant with respect to the variance of r j x0 , the ﬁrst term on the right-hand side Var½r ¼ Ex ½Var½rjx þ Var½Er ½rjx (9) of Eq. (16) can be rewritten as h i by the law of total variance [Eq. (5)]. Ex0 ½Var½wðx0 Þrjx0 ¼ Ex0 w2 ðx0 ÞVar½rjx0 : (17) In this context, the ﬁrst term in Eq. (9) has the physical inter- pretation as the variance due to the randomness of the transport Applying the deﬁnition for the expectation operator Ex0 to Eq. process (transport variance), whereas the second term is the vari- (17) yields ance due to the randomness of the initial birth state of each source Z particle (source variance). Ex0 ½Var½wðx0 Þrjx0 ¼ Var½rjx0 w2 ðx0 Þp0 ðx0 Þdx0 : (18) Finally, we note that the variance of the sum of uncorrelated G variables is equal to the sum of the variance of the individual random variables (a relationship referred to in some texts as the Using Eq. (13) to relate p(x0 ) ¼ w(x0 )p0 (x0 ) and switching the Bienayme formula), which allows the variance of the estimator for dummy variable of integration from x0 to x allows Eq. (18) to be the mean response m b r to be written as rewritten as D.P. Griesheimer, V.S. Sandhu / Nuclear Engineering and Technology 49 (2017) 1172e1180 1175 Z Z Ex0 ½Var½wðx0 Þrjx0 ¼ wðxÞVar½rjxpðxÞdx; (19) Var½wðxÞr Var½rjx ¼ pðxÞdx; (24) E½r Er ½rjx G G which can be written in expectation notation as which shows that the relative variance of the total response is equal to the average of the relative variance for all of the conditional responses taken with respect to the true source distribution p(x). Ex0 ½Var½wðx0 Þrjx0 ¼ Ex ½wðxÞVar½rjx: (20) Note that the use of importance sampling for the source state Substituting Eq. (18) into Eq. (16) gives an expression for the does not eliminate the variability in observed response owing to response variance under importance sampling. the transport process (transport variance). In order to achieve a true zero-variance process, it is necessary to also adjust the transport process so that each initial state has the same observed response, Var½wðx0 Þr ¼ Ex ½wðxÞVar½rjx þ Var½wðx0 ÞEr ½rjx0 : (21) not just the same expected response as shown above. Note that the ﬁrst term on the right-hand side of Eq. (21) is In consistent-source variance reduction methods, the source expressed in terms of the expectation and variance with respect to weighting function is chosen based on Eq. (23), in order to mini- the original source distribution, p(x), rather than the modiﬁed mize (or eliminate) the source variance term. It follows that this distribution p0 (x). importance weighting function does not necessarily minimize the formula for the variance of the total variance of the response, as demonstrated by the example Finally, applying the Bienayme problem described in the results section. For simpliﬁed problems, it sum of uncorrelated variables to Eq. (14) and using the result for is possible to show that a weighting function obtained by mini- Var[w(x′) r] gives the ﬁnal variance for the statistic m br mizing both the transport and source variance terms in Eq. (21) produces a lower total variance. Generalized conditions for the h i E ½wðxÞVar½rjx Var½wðx0 Þ ½rjx0 optimum weighting function have not been derived. However, it x Var m^ imp r ¼ þ r : (22) should be noted that results from preliminary testing covering a N N range of possible model conditions suggest that the zero-source- Comparing Eqs. (22) and (11) illustrates the effect of importance variance weighting function given in Eq. (23) produces responses sampling via the presence of the weight parameter w(x) in each within a few percent of the optimal variance for realistic situations. term. It is particularly interesting to note that the weighting parameter affects both the transport variance and source variance 2.4. Source splitting scheme terms. The effect on the transport variance term is somewhat sur- prising, as it suggests that importance sampling schemes that Next, let us consider a generalization of the expected response emphasize initial birth states with above-average response estimator in which source particles are split at birth and multiple (transport) variance may diminish the effectiveness of importance independent realizations are generated for each initial state point. sampling for reducing total variance. For an MC simulation with N0 independently sampled state The objective of the importance sampling scheme described in points and M independent realizations for each initial state point, this section is to reduce the total variance of the response by the unbiased importance sampling statistic for estimating the ex- altering the distribution of source states in an unbiased way. When pected response is importance sampling is applied to a process that involves a prob- ability distribution of only one variable, say p(x), it is known that N0 N0 1 X 1 XM 0 ~ 1 X ~zsplit ; the optimal weighting function for calculating the expected value b split m r ¼ wðx Þ r j;i ≡ (25) of any response function g(x) is given by w(x) ¼ E[g(x)]/g(x). In fact, N 0 i¼1 M j¼1 i N 0 i¼1 M;i this optimal weight gives Var[w(x)g(x)] ¼ 0 by ensuring that each sample gives a response that is exactly equal to the expected where ~rj;i is a realization from the conditional distribution p(rjx0 i), response E[g(x)]. split However, the identiﬁcation of an optimal weighting function for and ~zM;i is a realization of the random variable obtained by taking importance sampling is not as straightforward when dealing with a the sample mean of the responses due to transporting M replicates joint probability distribution. By analogy with the single-variable ~0i . Note that the variable zsplit of initial source state x M is actually a case, it seems reasonable to assume that the weight should be function of a realization taken from the joint probability distribu- deﬁned as tion p(x,r1, …, rM). Such a realization is not physically realistic, because it implies that a single source site will induce M response Z realizations. However, the joint realization is an accurate repre- Er ½rjxpðxÞdx sentation of MC particle transport with source particle splitting. G Furthermore, we note that the response variables r1, …, rM are wðxÞ ¼ ; (23) conditionally independent, which means that each variable is Er ½rjx conditionally dependent on the random variable x, but indepen- which is the ratio of the mean total response to the mean response dent from all of the other response variables. The property of conditioned on the initial state x. In fact, it can be immediately conditional independence allows the joint PDF for the birth state shown that the weighting function deﬁned in Eq. (23) optimizes and M responses to be written as the variance associated with the selection of an initial source state Y M pðx; r1 ; …; rM Þ ¼ pðr1 ; …; rM jxÞpðxÞ ¼ pðxÞ p rj x : (26) by eliminating the source variance term (Var[Er[w(x′)rjx′] ¼ 0). j¼1 Again, this is because the weighting function ensures that each source state contributes exactly the same expected response. Further recognizing that all of the response realizations are Applying the weight function in Eq. (23) to the equation for total taken from a common conditional PDF, p(rj j x) ¼ p(r j x) for all rj, variance of the response, Eq. (21), yields allows the joint PDF in Eq. (26) to be written as 1176 D.P. Griesheimer, V.S. Sandhu / Nuclear Engineering and Technology 49 (2017) 1172e1180 h i E ½wðxÞVar½rjx þ Var½wðx0 ÞEr ½rjx0 : split x pðx; r1 ; …; rM Þ ¼ pðxÞðpðrjxÞÞM : (27) Var zM ¼ (37) M Returning to the expected response estimator, m b split r , deﬁned in Inspection of Eq. (37) indicates that source splitting only affects Eq. (25), it follows that the corresponding sample variance statistic the transport variance term, and that the total response variance per source sample will decrease as the splitting factor M increases. 0 12 Equation (37) gives the variance for the mean response of M 2 1 X N0 @1 X M 0 split A replicates for a single source state. The variance for the mean b s r;split ¼ w xi r j;i m ~ br ; (28) N 0 1 i¼1 M j¼1 response over N0 independent trials, mb split , can be determined by r applying the Bienayme formula to Eq. (25), then substituting the is an unbiased estimator for Var [zsplit M ], where expression in Eq. (37), to yield h i h h ii h h ii h i E ½Var½wðxÞrjx Var½wðx0 ÞE ½rjx0 split split Var mb split ¼ x þ r : ¼ Ex0 Var zM x0 þ Var Er zM x0 : split Var zM (29) r MN 0 N0 (38) Fortunately, Eq. (29) can be recast into a more meaningful form. A comparison of Eq. (38) with the variance for importance We begin by expanding the ﬁrst term (transport variance) on the sampling without splitting [Eq. (22)] shows that the only difference right-hand side of Eq. (29) between the expressions is the presence of the 1/M factor in the ﬁrst term of Eq. (38). Note that if the number of source samples are 2 2 33 h h ii X held constant between importance sampling with and without split 0 1 M Ex0 Var zM x ¼ Ex0 4Var4 wðx0 Þrj x0 5 5: (30) source splitting (e.g., N0 ¼ N), source splitting will cause the overall M j¼1 response variance to go down, because of a decrease in the trans- port variance term. This makes sense, as the source replicates are Recognizing that the number of replicates, M, and the weight reducing the statistical uncertainty in the response associated with factor w(x) are independent of the replicate number and response transporting radiation from each sampled initial particle state. variance allows Eq. (30) to be rewritten as However, if N0 and M are constrained such that the total amount 2 2 33 of transport work is held constant (MN0 ¼ N), a question arises h h ii regarding optimal allocation of resources between M and N0 . In split 0 w 2 ðx 0 Þ XM 0 Ex0 Var zM x ¼ Ex0 4 Var4 r 5 5: j x (31) M2 cases where the source variance is larger than the transport vari- j¼1 ance, the optimal allocation is to maximize the number of source Because the realizations of ri,j (i.e., rjjx) are conditionally inde- particles, N0 ¼ N, and use only one replicate for each source. formula to write However, as the magnitude of the source variance becomes small pendent, it is possible to apply the Bienayme relative to the transport variance, the penalty (i.e., increase in total h h i i E 0 w2 ðx0 ÞVar½rjx0 variance) associated with source splitting decreases. In the limiting 0 x Ex0 Var zsplit M x ¼ : (32) case, when the source variance is equal to zero (i.e., all weighted M source sites produce the same expected response), the total Finally, we apply the procedure outlined in Eqs. (18e20) to write response variance is not affected by source splitting. Eq. (32) in terms of the expected value with respect to the unbiased This observation leads to several interesting conclusions. First, source distribution, p(x) this result conﬁrms the long-held conventional wisdom that source splitting will increase the variance of tallied quantities in simula- h h i i E ½wðxÞVar½rjx 0 x tions for which it is used, assuming that total work is held constant. Ex0 Var zsplit M x ¼ : (33) M However, the result also shows that the increase in variance due to splitting may be small in cases where the source variance is small Returning to the second term (source variance) in Eq. (29), we relative to the transport variance. This second result is especially can again expand this term and factor out constants to give intriguing when recalling that importance sampling is used to 2 2 33 minimize the source variance term. Thus, in this regime, source h h ii 0Þ XM split 0 wðx splitting may be used to correct for an inconsistent source distri- Var Er zM x ¼ Var4 Er 4 rj x x0 5 5: 0 (34) M bution without a signiﬁcant increase in total response variance. j¼1 2.5. Source rouletting scheme Applying the deﬁnition of the conditional expected value for the response, Eq. (4), gives Let us now consider an expected response estimator in which 2 3 source particles are subjected to Russian roulette at birth. For an MC M Z ∞ h h ii 0Þ X simulation with N0 independently sampled state points, the unbi- wðx Var Er zM x0 ¼ Var4 rj p rj x0 drj 5: split (35) ased importance sampling statistic for estimating the expected M j¼1 0 response with Russian roulette is Because the random variables rj are conditionally independent and N0 ~ 0 drawn from the same probability distribution p(r j x), it follows that b roulette 1 X ~ i ~r i bi w x 1 XN0 ~zroulette ; m r ¼ 0 ≡ (39) N i¼1 M N i¼1 M;i 0 2 3 h h ii Z∞ 0 4wðx0 Þx0 Var Er zsplit M x ¼ Var rpðrjx0 Þdr 5 ~0i are realizations from the joint distribution p0 (x0 ,r), where ~ri and x (36) 0 ~ and b is a realization from a Bernoulli distribution with probability i ¼ Var½wðx0 ÞEr ½rjx0 : roulette of success given by M (for M 1), and ~zM;i is a realization of the Substituting Eqs. (33) and (36) into Eq. (29), gives the ﬁnal random variable obtained by setting the response equal to zero split simpliﬁed expression for the variance of zM with probability (1 e M). Note that the random variable b D.P. Griesheimer, V.S. Sandhu / Nuclear Engineering and Technology 49 (2017) 1172e1180 1177 determines whether each particle will survive the roulette process, common response estimator that can be used for both source and is completely independent from the random variables x and r. rouletting and splitting, including noninteger splitting ratios. In this case, the sample variance statistic For an MC simulation with N0 independently sampled state !2 points, the generalized sampling statistic for estimating the ex- 1 X N0 ~ wðx b ~ i 0 Þ~r i 2 pected response is b s r;roulette ¼ 0 i m roulette b r ; (40) N 1 i¼1 M 0 0 1 XN 1 X QMS ~0 wx 1 XN ~zgeneral ; is an unbiased estimator for Var[zroulette M ], where b general m r ¼ 0 b j ~ 0i ~r j;i ≡ 0 M;i (46) N i¼1 M N i¼1 " 2 # 2 j¼1 h i bwðx0 Þr bwðx0 Þr roulette Var zM ¼E E ; (41) M M where ~rj;i is a realization from the conditional distribution p(rjx0 i), ~ 0 is a realization from a Bernoulli distribution with probability of bj which can be factored to yield general success given by M/QMS, ~zM;i is a realization of the random var- h i 1 h i h i iable obtained by taking the sample mean of the responses due to ¼ 2 E b2 E ðwðx0 ÞrÞ2 E½b2 E½wðx0 Þr2 ; roulette Var zM M transporting QMS replicates of the initial source site x ~ 0 , and (42) the symbol QS denotes a ceiling operation such that QMS is the smallest integer that is larger than M. As described in Section 2.4, because the realizations of (x0 ,r) are independent from the re- general the variable zM is actually a function of a realization taken from alizations of b. Recognizing that E[b2] ¼ M and E[b]2 ¼ M2 for a the conditionally independent joint probability distribution Bernoulli distributed random variable enables Eq. (42) to be p(x,r ,…,r ). In addition, the Bernoulli realizations b~ 0 are assumed simpliﬁed, giving 1 M i to be unconditionally independent from the random variables x h i 1 1M and r. ¼ Var½wðx0 Þr þ roulette Var zM E½r2 : (43) M M Note that in Eq. (46) the parameter M determines how the weight adjustment is performed, with M < 1 causing source Again, the variance for the mean response over N0 independent roulette rouletting, M > 1 causing source splitting, and M ¼ 1 reverting to trials, m br , can be determined by applying the Bienayme formula standard importance sampling. It is also important to recognize to Eq. (39), then substituting the expression in Eq. (43) to yield that the estimator deﬁned in Eq. (46) automatically accounts for h i 1 1M noninteger splitting ratios by adjusting the effective splitting ratio Var mb roulette r ¼ Var½wðx0 Þr þ E½r2 : (44) upward to the integer value QMS and then applying rouletting to MN 0 MN0 eliminate the excess particles produced during splitting. For comparison with the corresponding variance for importance Based on the form of Eq. (46), it follows that the corresponding sampling with [Eq. (38)] and without splitting [Eq. (38)], the pre- sample variance statistic vious result for Var[w(x0 )r] given in Eq. (21) can be applied to Eq. 0 12 0 (44), yielding an alternate form 2 1 X N @ 1 X QMS 0 0 ~w x general A b s r;general ¼ 0 b j ~ i ~r j;i m br (47) h i E ½wðxÞVar½rjx N 1 i¼1 M j¼1 Var mb roulette ¼ x h i r MN0 general (45) is an unbiased estimator for Var zM , where 0 0 Var½wðx ÞEr ½rjx 1 M þ þ E½r2 : h i h h ii h h ii MN 0 MN 0 general general 0 general 0 Var zM ¼ Ex0 Var zM x þ Var Er zM x : From Eq. (44), it is clear that applying Russian roulette will al- ways cause the response variance to increase relative to importance (48) sampling for a ﬁxed number of source samples (N0 ¼ N). This in- Again, Eq. (48) can be simpliﬁed into a more intuitive form, crease is attributable both to a decrease in the denominator of each following the same general approach as applied in Section 2.4. term because M < 1, as well as the presence of an extra additive As usual, we begin by expanding the ﬁrst term on the right-hand term proportional to the square of the expected response. side of Eq. (48) and moving constant factors out of the variance If N0 and M are constrained such that the total amount of operator, which gives transport work is held constant (M N0 ¼ N), Eq. (44) shows that 2 2 33 rouletting will always produce an increase in response variance h h ii general 0 w2 ðx0 Þ X QMS 0 relative to a traditional importance sampling scheme [Eq. (22)]. Ex0 Var zM x ¼ Ex0 4 Var4 b0 5 5: 2 j x rj (49) This variance penalty appears as an additive term, which is pro- M j¼1 portional to (1 e M)$E[r]2. Unlike the splitting scheme considered previously, which was dependent on the biased source distribution Recognizing that the realizations from b0 are independent and used for importance sampling, the variance penalty for rouletting the rj variables are conditionally independent allows application of only depends on the roulette survival probability and the expected the Bienayme formula to Eq. (49), resulting in response. Interestingly, Eq. (44) indicates that there is no increase 2 0 h h ii w ðx ÞQMS in variance when rouletting particles in a system where the ex- general 0 0 0 Ex0 Var zM x ¼ Ex0 Var½b rjx : (50) pected response is zero. M2 2.6. Generalized weight adjustment scheme By the deﬁnition of variance and the independence of