{'fraction_non_alphanumeric': 0.028112922989914226, 'fraction_numerical': 0.009755867659534357, 'mean_word_length': 4.728536717062635, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Monte Carlo convergence of rival samplers

It is often necessary to make sampling-based statistical inference about many probability distributions in parallel. Given a finite computational resource, this article addresses how to optimally divide sampling effort between the samplers of the different distributions. Formally approaching this decision problem requires both the specification of an error criterion to assess how well each group of samples represent their underlying distribution, and a loss function to combine the errors into an overall performance score. For the first part, a new Monte Carlo divergence error criterion based on Jensen-Shannon divergence is proposed. Using results from information theory, approximations are derived for estimating this criterion for each target based on a single run, enabling adaptive sample size choices to be made during sampling.


Introduction
Let X 1 , X 2 , . . . , X n denote n independent random samples to be drawn from a probability distribution π with support X for the purpose of approximating π, when the form of this distribution is unknown. The corresponding random measure of this sample is the Monte Carlo estimator of π, This is a maximum likelihood estimator of π and so is consistent, lim n→∞Πn (B) = π(B) for all π-measurable sets B. Sometimes estimating the entire distribution π is of intrinsic inferential interest, with no limit on the functionals of π which might be of future interest. Alternatively, the random sample might be an intermediary update of a sequential Monte Carlo samplerwhere it is desirable that the samples represent the current target distribution well at each step.
Pointwise Monte Carlo errors are inadequate for capturing the overall rate of convergence of the realised empirical measureπ n to π. This consideration is particularly relevant if π is an infinite mixture of distributions of unbounded dimension: In this case it becomes necessary to specify a degenerate, fixed dimension function of interest before Monte Carlo error can be assessed (Sisson and Fan, 2007). This seems unnecessarily restrictive, particularly since the assessment of convergence will vary depending on which function is selected and that choice might be somewhat arbitrary.
In the work presented here, multiple target distributions will be considered in parallel. This scenario is frequently encountered in real-time data processing, where streams of data pertaining to different statistical processes are collected and analysed in fixed time-window updates. Decisions on how much sampling effort to allocate to each target will be made sequentially, based on the apparent relative complexity of the targets; the complexities of the targets will not be known a priori, but can be estimated from the samples which have been obtained so far. As a consequence, the size of the sample n drawn from any particular target distribution will be a realisation of a random variable N, determined during sampling by an implied random stopping rule governed by the history of samples drawn from that target and those obtained from the other targets.
To extend the utility of Monte Carlo error to entire probability measures, the following question is addressed: If a new sample of random size N were drawn from π, how different toπ n might the new empirical measure be? If repeatedly drawing samples in this way led to relatively similar empirical measures, this suggests that the target is relatively well represented by N samples; whereas if the resulting empirical measures were very different, then there would be a stronger desire to obtain a (stochastically) larger number of samples. To formally address this question, a new Monte Carlo divergence error is proposed.
Correctly balancing the sample sizes is a non-trivial problem. Apparently sensible, but ad hoc, allocation strategies can lead to extremely poor performance, much worse than simply assigning the same number of samples to each target. Here, a sample-based estimate of the current Monte Carlo divergence error of an empirical measure is derived; the errors of each sampler are combined across samplers via a loss function, leading to a fully-principled, sequential sample allocation strategy. Good performance is demonstrated in both univariate and variable dimensional problems.
Section 2 formally defines and justifies Monte Carlo divergence as an error criterion. Next, Section 3 examines the implications of two different loss functions for combining sampler errors into a single performance score. Section 4 introduces some alternative sample size selection strategies; some are derived by adapting related ideas in the existing literature, and some are ad hoc. The collection of strategies are compared on univariate and variable dimension target distributions in Section 5 before a brief discussion in Section 6.

Monitoring Monte Carlo convergence
In this section the rate of convergence of the empirical distribution to the target will be assessed by information theoretic criteria. In information theory, it is common practice to discretise distributions of any continuous random variables (see Paninski, 2003). Without this discretisation (or some alternative smoothing method) the overlap of any two sets of independent samples would be empty, and distribution-free comparisons of the empirical measures from samples therefore rendered meaningless. For example, the Kullback-Leibler divergence between two independent realisations of (1) will be always be infinite.
When a target distribution relates to that of a continuous random variable, a common discretisation of both the empirical measure and notionally the target will be performed. For the rest of the article, bothπ and π should be regarded as suitably discretised approximations to the true distributions when the underlying variables are continuous. For univariate problems a large but finite grid with fixed spacing will be used to partition X into bins; for mixture problems with unbounded dimension, the same strategy will be used for each component of each dimension, implying an infinite number of bins. Later in Section 3.4, consideration will be given to how the number of bins for each dimension should be chosen.

Monte Carlo divergence error
For a discrete target probability distribution π, letΠ n be the estimator (1) for a prospective sample of size n to be drawn from π, and letΠ be the same estimator when the sample size n is a random stopping time.
For n ≥ 1, the Monte Carlo divergence error of the estimatorΠ n will be defined as where H is Shannon's entropy function; recall that if p = (p 1 , . . . , p K ) is a probability mass function, H(p) = − K i=1 p i log(p i ). The Monte Carlo divergence error e KL,n has a direct interpretation as the expected Kullback-Leibler divergence of the empirical distribution of a sample of size n from the target π, and so provides a natural measure of the adequacy ofΠ n for estimating π.
The Monte Carlo divergence error of the estimatorΠ when n is a random stopping time is defined as the expectation of e KL,n with respect to the stopping rule, or equivalently where the expectation is now with respect to both π and the stopping rule. This more general definition of Monte Carlo divergence error should be interpreted as the expected Kullback-Leibler divergence of the empirical distribution of a sample of random size from the target π.
To provide a sampling based justification for this definition of Monte Carlo divergence error, for M > 1 consider the empirical distribution estimatesπ 1 , . . . ,π M which would be obtained from M independent repetitions of sampling from π, where the sample size of each run is a random stopping time from the same rule.
The Jensen-Shannon divergence (Lin and Wong, 1990;Lin, 1991) ofπ 1 , . . . ,π M , measures the variability in these distribution estimates by calculating their average Kullback-Leibler divergence from the closest dominating measure, which is their average. The Jensen-Shannon divergence is a popular quantification of the difference between distributions, and its square root has the properties of a metric on distributions (Endres and Schindelin, 2003). Just as Monte Carlo variance is the infinite limit of the sample variance of M sample means, the Monte Carlo divergence error defined in (3) is easily seen to be the limit of (4) as M → ∞: By the strong law of large numbers, lim M →∞ , of the runs. It follows that (4) is a finite sample, consistent estimate of (3).
Finally, it should be noted that there is a second interpretation of the Monte Carlo divergence error: e KL is also the negative bias of the maximum likelihood estimator of the entropy of π given a random sample. In the next section it will be shown that this alternative interpretation is very useful, since it leads to a mechanism for estimating e KL from a single sample.

Estimating Monte Carlo divergence error
Whilst H(Π n ) is the maximum likelihood estimator of H(π), it is known to be a negatively biased since H is a concave function (Miller, 1955). Various bias corrections have been proposed in the information theory literature, and these corrections themselves can serve as approximately unbiased estimates of e KL,n . Also note that any unbiased estimate of e KL,n is also an unbiased estimate of e KL , the error under the random stopping rule.
Given a sample of size n, the popular Miller-Madow method estimates the negative bias of the maximum likelihood estimate H(π) to be (K − 1)/(2n), where K is the number of nonempty bins inπ. This estimate proves to be too crude for the purpose here of estimating e KL,n , since any two distributions with the same number of represented bins would be estimated to have the same divergence, regardless of how uniform the corresponding bin probabilities might be.
An improvement on the Miller-Madow estimate was provided by Grassberger (1988Grassberger ( , 2003, where n i is the number of samples in the ith nonempty bin, such that K i=1 n i = n, and φ(n i ) = n i {log(n i ) − ψ(n i )}, where ψ is the digamma function. Here (5) will provide an approximately unbiased estimate of the expected Kullback-Leibler divergence of the empirical distribution from the target, e KL,n , based on a single run of the sampler.

Efficient calculation
Calculation of (5) during sampling can be updated at each iteration very quickly, using the following equations. Let i ′ be the bin in which the nth observation falls and ∆ j be the jth forward difference operator. Thenê As well as estimating the current Monte Carlo divergence error of a distribution estimate after n samples from π, it will be of interest to estimate the expected reduction in error that would be achieved from obtaining one more sample, Assuming the probability of this new sample falling into the ith bin is approximated by the empirical, maximum likelihood estimate n i /n, the expected reduction in error from a further sample can be estimated asδ Again considering this calculation iteratively, if the nth observation falls into bin i ′ then δ KL,n = (n − 1)nδ KL,n−1 − n i ′ ∆ 2 φ(n i ′ − 1) n(n + 1) .

Alternative formulations of bias estimation
It should be noted that further refinements (additive terms) to the bias estimate of (5) are provided by Grassberger (2003), such as However, these additional terms, which arise as part of a second order approximation of an integral, are unstable, oscillating between positive and negative values. In this context, without careful treatment such terms can incorrectly suggest that the expected error might very slightly increase by taking further samples, which in practice is not true. Furthermore, due to their oscillating sign, these terms do not affect the overall drift of the function, which will be the quantity of longer term interest when deciding whether more sampling effort should be afforded.

Rival samplers
Consider m probability distributions π 1 , . . . , π m . For each π j , suppose there is a corresponding sampler which might be drawing independent samples directly from π j , or indirectly via suitably thinned Markov chain Monte Carlo simulation. importance Random samples are to be drawn from each of these samplers, with the empirical distributions then serving as approximations of the corresponding target distributions. Given a fixed computational resource, which might simply correspond to a final total number of random samples n, the decision problem to be addressed is how best to divide those n samples between the m samplers. That is, how sample sizes n (1) , . . . , n (m) for the m targets should be chosen subject to the constraint m j=1 n (j) = n. The samplers can be viewed as rivals to one another for the same fixed computational resource.
Without this constraint, the problem would be ill-posed: for all j, n (j) should be chosen to be as large as possible, since Monte Carlo errors are monotonically decreasing with sample size. A constraint is required for sample size choice to be practically meaningful. In contrast, results which establish the minimum sample size for which errors should fall within a (typically arbitrary) desired level of precision are theoretically interesting, but are perhaps best viewed in the reverse direction in this context; given the inevitable usage of the maximum allowable computation time, understanding the level of error this limit implies.
The default choice is for equal samples sizes, n (j) = n/m, but such an approach disregards any differences in the complexities of the target distributions, which in general could be arbitrarily different. Here the aim is to adaptively determine how much sampling effort should be afforded to each sampler. The preceding section established a general method for assessing the error of each sampler. The choice of how to balance sample sizes between the samplers will be made according to a loss function of those errors.

Loss functions for Monte Carlo errors
Suppose that the decision to assign n (j) of the total n samples to the sampler of π j implies a Monte Carlo error level e Two natural alternative loss functions for combining the individual errors e (j) n (j) into an overall performance error are considered here. Other choices, or indeed linear combinations of these two losses, could be examined. One possibility is that utility could be derived from controlling the average error across the samplers, suggesting an (expected) loss function of the form Alternatively it might be important to control the maximum error of the m samplers, suggesting a different loss function These two loss functions will imply different sample size strategies. As a simple illustration, consider an example of estimating the means of two target distributions with known variances σ 2 1 , σ 2 2 , with error measured by the Monte Carlo variance of the sample means; under L ave , the optimal ratio of the sample sizes from the two targets, n (1) /n (2) , would be given by the ratio of the standard deviations σ 1 /σ 2 , whereas under L max the optimal sample size ratio would be given by the ratio of the variances σ 2 1 /σ 2 2 . Therefore care should be exercised in specifying the correct form of loss function required.

Sequential allocation of samples
Away from the stylised example of the previous section, in general it is more likely that little will be known a priori about the target distributions being sampled. Instead, the aim will be to dynamically decide, during sampling, which samplers should be afforded more sampling effort, conditional on the information learnt so far about the targets. A sequential decision approach is taken. Having taken n ′ < n samples, with n (j) of these allocated to the jth sampler, the decision problem is to choose from which sampler to draw the (n ′ + 1)th sample, such that the chosen loss function of the estimated Monte Carlo errors of the samplers {ê In this sequential setting, the operational difference between the loss functions L max and L ave becomes clearer. If the aim is to minimise L max , then the optimal decision for allocating one more sample is to allocate it to the sampler with the highest estimated error, since error is a decreasing function of sample size. Alternatively, if minimising L ave , then the new sample should be allocated to the sampler for which the estimated decrease in error is highest, since this will minimise the overall sum. Note that this is a myopic decision rule, looking only one step ahead. There are three reasons why this is preferred; first, considering optimal sequences of future allocations leads to a combinatorial explosion unsuitable for a method intended for optimising the use of a fixed computational resource; second, the final number of samples may even be unknown; third, the estimated error or expected change in error under the Monte Carlo divergence criterion (5), (8) can be updated very efficiently via (6) and 9: after one more sample, only one bin count n (j) i for one sampler j will have changed.
Any sequential sampling allocation scheme which depends on the outcomes of the random draws will imply a random stopping rule for the number of samples eventually allocated to each sampler. This adds an extra complication, since some stopping rules will introduce bias into Monte Carlo estimates (Mendo and Hernando, 2006). Here this bias arises if the first samples taken from a target distribution have a particularly low estimated Monte Carlo error, as this will cause the other rival samplers to share all of the remaining samples; without corrective action, this phenomenon causes Monte Carlo estimators to be biased towards estimates of this character.
When the Monte Carlo error is the divergence measure (3), low error estimates correspond to low entropy empirical measures, which can spuriously arise if the first random samples happen to fall into the same bin. Therefore, to eradicate this bias, a minimum number of samples ℓ (j) is recommended for each target distribution, to prevent degenerate sample sizes. To examine stability, these minima can be chosen in increasing steps until the resulting samples sizes converge. For the examples in Section 5, due to the relatively fine grid used for binning samples it was enough to set ℓ (j) = 500 to obtain convergence.

Algorithm: Rival sampling
The full algorithm for sequential sampling from rival target distributions to minimise estimated loss is now presented. For m samplers of target distributions π 1 , . . . , π m , let ℓ (j) ≥ 1 be the minimum number of samples that should be drawn from π j . Let L ∈ {L max , L ave } be the chosen loss function for combining Monte Carlo errors across the samplers.
The algorithm proceeds as follows: 1. Initialise-For j = 1, . . . , m: (a) Draw ℓ (j) samples from π j and calculateπ j , the binned empirical estimate of π j ; let K j be the number of non-empty bins inπ j , and n (j) 1 , . . . , n K j be the corresponding bin counts; set n (j) = ℓ (j) .

Choosing a bin width for discretisation
The algorithm of Section 3.3 requires a method of discretising samples from continuous distributions. (For simplicity, a fixed bin width can be assumed for each dimension of a multivariate distribution.) The following observations offer some insight for what makes a good bin width in this context. In the limit of the bin width going to zero, the binned empirical distribution after n samples would have n non-empty bins each containing one sample. Although the identity of those bins would vary across samplers, these empirical distributions would be indistinguishable in terms of both entropy and (3); so each sampler would be allocated the same sample size. In the opposite limit of the bin width becoming arbitrarily large, all samples of the same dimension would fall into the same bin. For fixed dimension problems, this would mean all sample sizes would be equal, and otherwise in transdimensional problems the strategies would simplify to working with marginal distributions of the dimension, which reduces the potential diversity of sample sizes. So a good bin width would lie well within these two extremes, ideally maximising the resulting differences in sample sizes. That is, a good bin width should distinguish well the varying complexity of the different samplers. Further to these observations, the next section suggests a novel maximum likelihood approach for determining an optimal number of bins, which in principle could be deployed adaptively as the algorithm progresses.

Maximum likelihood bin width estimation for Bayesian histograms
Consider a regular histogram of K equal width bins on the interval [a, b], and let p = (p 1 , . . . , p K ) be the bin probabilities. The Bayesian formulation of this histogram (?) treats the probabilities p as unknown, and a conjugate Dirichlet prior distribution based on a Lebesgue base measure with confidence level α presumes p ∼ Dirichlet(α(b − a)/K · 1 T ). For n samples, the marginal likelihood of observing bin counts n 1 , . . . , n K under this model is Using standard optimisation techniques, identifying the pair (K,α) that jointly maximise (15) suggests thatK serves as a good number of bins for a regular histogram of the observed data.

Alternative strategies
To calibrate the performance of the proposed method, some variations of the strategy for selecting sample sizes are considered. This section considers some alternative measures of the Monte Carlo error of a sampler, to be used in place of the divergence estimatesê

Chi-squared goodness of fit statistic
In the context of particle filters, Fox (2003) proposed a method for choosing the number of samples n required from a single sampler to guarantee that, under a chi square approximation, with a desired probability (1 − δ) the Kullback-Leibler divergence between the binned empirical and true distributions does not exceed a certain threshold ε. This was achieved by noting an identity between 2n times this divergence and the likelihood ratio statistic for testing the true distribution against the empirical distribution, assuming the true distribution had the same number of bins, K, as the observed empirical distribution. Since the likelihood ratio statistic should approximately follow a chi-squared distribution with K − 1 degrees of freedom, this suggested a sample size of where χ 2 K−1,1−δ is the 1 − δ quantile of that distribution. Adapting this idea to the algorithm of Section 3.3 simply requires a rearrangement of (16) to give the approximate error as a function of sample size, This error estimate can be substituted directly into the algorithm in place of the Monte Carlo divergence error estimateê (j) n (j) to provide an alternative scheme for choosing sample sizes when using loss function L max . The same (arbitrary) value of δ must be used for each rival sampler, and here this was specified as δ = 0.05 although the results are robust to different choices.
By the central limit theorem, the chi-squared distribution quantiles grow near-linearly with the degrees of freedom parameter for K > 30 (Fisher, 1959), so it should be noted that (17), which depends only on the number of bins, has much similarity, and almost equivalence, with the Miller-Madow estimate of entropy error cited in Section 2.2. By the reasoning given in Section 2.2, use of this error function should show some similarity in performance with the proposed method, but be less robust to distinguishing differences in distributions beyond the number of non-empty bins.
Recall from Section 3.2 that the sequential allocation strategy for minimising the loss function L ave requires an estimate of the expected reduction in error which would be achieved from obtaining another observation from a sampler. Since this error criterion depends entirely upon the number of non-empty bins K, in this case an estimate is required for the probability of the new observation falling into a new bin. A simple empirical estimate of the probability of falling into a new bin is provided by the proportion of samples after the first one that have fallen into new bins, given by (K − 1)/(n − 1). Note that this estimate will naturally carry positive bias, since discovery of new bins should decrease over time, and so a sliding window of this quantity might be more appropriate in some contexts.

Reference points
As a convergence diagnostic for transdimensional samplers, Sisson and Fan (2007) proposed running replicate sampling chains for the same target distribution, and comparing the variability across the chains of the empirical distributions of a distance-based function of interest. The method requires that the target π be a probability distribution for a point process, and maps multidimensional sampled tuples of events from π to a fixed-dimension space. Specifically, a set of reference points V are chosen, and for any sampled tuple of events the distance from each reference point to the closest event in the tuple is calculated. Thus π is summarised by a |V|-dimensional distribution, where |V| is the number of reference points in V.
One example considered in Sisson and Fan (2007) is a Bayesian continuous-time changepoint analysis of a changing regression model with an unknown number of changepoint locations. A variation of this example is analysed in Section 5.1.2 in this article, where instead the analysis will be for the canonical problem of detecting changes in the piecewise constant intensity function λ(t), t ∈ [0, 1] of an inhomogeneous Poisson process (see Raftery and Akman, 1986, and the subsequent literature). For Poisson process data, there are two natural functions of interest which could be evaluated at each reference point. The first is the distance to the nearest changepoint, the second is the intensity level. Both will be considered in Section 5.1.2. Note that in Sisson and Fan (2007), reference points are selected from random components from an initial sample from the single target distribution. Here, since there are multiple target distributions, a grid of one hundred uniformly spaced points across the domain [0, 1] is used.
The convergence diagnostic of Sisson and Fan (2007) did not formally provide a method for calibrating error or selecting sample size. Here, to compare the performance of the proposed sample size algorithm of Section 3.3, the sum across the reference points of the Monte Carlo variances of either of these functions of interest is used as the error criterion in the algorithm.

Ad hoc strategies
To demonstrate the value of the sophisticated sample size selection strategies given above, two simple strategies which have similar motivation but are otherwise ad hoc are included in the numerical comparisons of Section 5. These strategies are now briefly explained.

Extent
The extent of a distribution is the exponential of its entropy, and was introduced as a measure of spread by Campbell (1966). A simple strategy might be to choose sample size proportional to the estimated square of the extent of π, exp{2H(π)}. Note that the Gaussian distribution N(µ, σ 2 ), has an extent which is directly proportional to the standard deviation σ, and so in the univariate Gaussian example which will be considered in Section 5.1.1, this sample allocation strategy will be approximately equivalent to the optimal strategy when minimising the maximum Monte Carlo error of the sample means (cf. Section 3.1). Robert and Casella (2005) present a class of convergence tests for monitoring the stationarity of the output of a sampler from a single run which operate by splitting the current sample (x 1 , x 2 , . . . , x n ) in two and quantifying the difference between the empirical distributions of the first half of the sample (x 1 , x 2 , . . . , x ⌊n/2⌋ ), and the second half (x ⌊n/2⌋+1 , x ⌊n/2⌋+2 , . . . , x n ). For univariate samplers the Kolmogorov-Smirnov test, for example, is used to obtain a p-value as a measure of evidence that the second half of the sample is different from the first, and hence neither half is adequately representative of the target. The test statistics which are used condition on the sample size, and so the sole purpose of these procedures is to investigate how well the sampler is mixing and exploring the target distribution.

Jensen-Shannon divergence
To adapt these ideas to the current context, any mixing issues can first be discounted by splitting the sample for each target in half by allocating the samples into two groups alternately, so that the distribution of, say, (x 1 , x 3 , . . . , x n−1 ) can be compared with the distribution of (x 2 , x 4 , . . . , x n ). This method of splitting up the sample is also computationally much simpler in a streaming context, as incrementing the sample size n does not change the required groupings of the existing samples. Letπ odd n andπ even n be the respective empirical distributions of these two subsamples. A crude variation on using the Monte Carlo divergence error criteria of (5) is to estimate the error of the sampler by the Jensen-Shannon divergence ofπ odd n andπ even n , e JSD,n = JSD(π odd n ,π even n ) = H(π n ) − H(π odd n ) + H(π even n ) 2 .
If sufficiently many samples have been taken forπ n to be a good representation of the target distribution, then both halves of the sample should also provide reasonable approximations of the target and therefore have low divergence between one another. As in Section 2.2.1, calculation of (18) during sampling can be updated at each iteration very quickly. Let i ′ be the bin in which the nth observation falls. Then, for example, updating the first term of (18) simply requires H(π n ) = H(π n−1 ) − n i ′ log(n i ′ ) + (n i ′ − 1) log(n i ′ − 1) + log(n) − log(n − 1).

Examples
The methodology from this article is demonstrated on three different data problems. The first two examples assume only two or three data processes respectively, to allow a detailed examination of how the allocation strategies differ. Then finally a larger scale example with 400 data processes is considered, derived from the IEEE VAST 2008 Challenge concerning communication network anomaly detection.

Small scale examples
Two straightforward, synthetic examples are now considered. The first is a univariate problem of fixed dimension with two Gaussian target distributions, and the second is a transdimensional problem of unbounded dimension, concerning the changepoints in the piecewise constant intensity functions of three inhomogeneous Poisson processes. In both examples, it is assumed that a priori nothing is known about the target distributions and that computational limitations determine that only a fixed total number of samples can be obtained from them overall, which will correspond to an average of 50,000 samples per target distribution. Both loss functions from Section 3.1 are considered, measuring either the maximum error or average error across the target samplers. For each loss function, the following sample size allocation strategies are considered: 1. "Fixed" -The default strategy, 50,000 samples are obtained from each sampler.
2. Dynamically, aiming to minimise the expected loss, with sampling error estimated using the following methods: (a) "Grassberger" -Monte Carlo divergence error estimation from Section 2.2; (b) "Fox" -the χ 2 goodness of fit statistic of Section 4.1; (c) "Sisson" (only for the transdimensional example) -the Monte Carlo variances of one of the two candidate fixed dimension functions from Section 4.2 evaluated at 100 equally spaced reference points ("Sisson-i", for the intensity function, "Sisson-n" for the distance to nearest changepoint function); (d) "Extent" and "JSD" -two ad hoc criteria from Section 4.3.
Each sample size allocation strategy is evaluated over a large number of replications M, where M = 1 million or M = 100,000 respectively in the two examples. Good performance of a sample allocation strategy is measured by the chosen loss function when applied to the realised Monte Carlo divergence error e KL for each sampler. Good estimates of the true values of e KL are obtained by calculating the Jensen-Shannon divergence of the Monte Carlo empirical distributions obtained from the M runs (cf. Section 2.1).
Note that in all simulations, the same random number generating seeds are used for all strategies, so that all strategies are making decisions based on exactly the same samples.
This choice corresponds to an interior range of plus or minus five times the largest of the standard deviations of the two targets, divided into 100 evenly spaced bins, along with two extra bins for the extreme tails. Results were robust to allowing wider ranges or more bins. For further validation, a simple experiment was conducted using the method from Section 3.4.1 on [−10, 10]: The maximum number of samples (100,000) were simulated from both π 1 and π 2 , leading to estimateŝ K = 92 andK = 76 respectively, suggesting 100 as a good number of bins for fitting these densities.
The varying sample sizes obtained from each target from the 1 million simulations using each of the sample allocation strategies listed above and the loss function L max are shown in Fig. 1. Tables 1 and 2 show the mean sample sizes and the implied Monte Carlo divergence error for each target distribution using each of the sample allocation strategies listed above; the two tables correspond to the two choices of loss function for combining errors. Table 1 gives the results under the loss function (12) which calculates the maximal error across samplers. Optimal performance would imply approximately equal Monte Carlo divergence errors for the two targets, and the proposed strategy based on Grassberger's entropy bias estimate is by far the closest to achieving this objective. Interestingly, note that under this best strategy, the average sample sizes are almost exactly in the ratio 1:2, the same ratio as the true standard deviations of the target distributions. Recall from Section 3.1, that such a ratio is optimal in another sense, minimising the average Monte Carlo errors of the sample means. This contrast highlights the importance of carefully specifying the desired error criterion as well as the correct loss function.
One of the two ad hoc strategies based on calculating the Jensen-Shannon divergence between the two halves of the sample is only slightly outperformed by the χ 2 goodness of fit method; however, note in Figure 1 the much higher variance of the sample sizes with the JSD method, which is indicative of an unreliable strategy. The other ad hoc method which takes sample sizes proportional to the extent of the empirical distributions is seen to overcompensate for the higher variance of the second Gaussian, and performs worse than even the default equal sample size strategy. For that strategy, note the sample sizes are almost exactly in the ratio 1:4, the same ratio as the true variances of the target distributions. Recall from Section 4.3.1 that such a strategy is approximately equivalent to minimising the maximum Monte Carlo errors of the sample means, which was noted in Section 3.1 to imply such an allocation ratio. Table 2 gives the results under the loss function L ave (11) which calculates the average error across samplers. The Monte Carlo divergence strategy based on Grassberger's entropy bias estimate performs best, although the the χ 2 goodness of fit method also performs very well here. The contrasting sample sizes between the loss functions L max and L ave for all dynamic allocation strategies are noteworthy, as remarked in Section 3.1.

Transdimensional mixture target distributions
For a more complex example, simulated data were generated from three inhomogeneous Poisson processes on [0, 1] with different piecewise constant intensity functions. In each case, prior beliefs for the intensity functions were specified by a homogeneous Poisson process prior distribution on the number and locations of the changepoints and independent, conjugate gamma priors on the intensity levels. The three rival target distributions for inference are the Bayesian marginal posterior distributions on the number and locations of the changepoints for each of the three processes. Each of the three Poisson processes had two changepoints, located at 1/3 and 2/3; the intensity levels of the three processes were respectively (200, 300, 400), (200, 350, 500), (200, 400, 600), so the processes differed only through magnitudes of intensity changes. To make the target distributions closer and therefore make the inferential problem harder, in each case the prior expectation for the number of changepoints was set to 1. For illustration of the differences in complexity of the resulting posterior distributions for the changepoints of the three processes, large sample esti-  mates of the true, discretised posterior distributions are shown in Fig. 2, based upon one trillion reversible jump Markov chain Monte Carlo samples (Green, 1995). In all cases there is insufficient information to strongly detect both changepoints, and so much of the mass of the posterior distributions is localised at a single changepoint at 1/2, the midpoint of the two true changepoints.
Additionally, Fig. 3 shows the posterior variance of two functions of interest identified in Section 4.2 for t ∈ [0, 1]: the distance to the nearest changepoint, and the intensity level.
To determine sample sizes, reversible jump Markov chain Monte Carlo simulation was used to sample from the marginal posterior distributions of the changepoints, and the chains were thinned with only one in every fifty iterations retained to give approximately independent samples. To discretise the distributions, the interval [0, 1] was divided into 50 equally sized bins; whilst for a single dimension this would be fewer bins than were used in the previous section, here the bins are applied to each dimension of a mixture model of unbounded dimension, meaning that actually a very large number of bins are visited; computational storage issues can begin to arise when using an even larger number of bins, simply through storing the frequency counts of the samples. Fig. 4 shows the distributions of sample sizes obtained from a selection of the strategies over M = 100,000 repetitions, and Tables 3 and 4 show results from the different strategies examined for these more complex transdimensional samplers. Performance is similar to the previous section, with the Grassberger entropy bias correction method performing best.
For this transdimensional sampling example, it also makes sense to consider the fixed dimen-   sion function of interest methods of Sisson and Fan (2007), using the mean intensity function of the Poisson process or the distance to nearest changepoint, each evaluated at 100 equally spaced grid points on [0, 1]. The Monte Carlo variances used in these strategies estimate the variances displayed in the plots of Fig. 3 at the reference points, divided by the current sample size. The performance of these fixed dimensional strategies is particularly poor under the loss function L max . Importantly, it should also be noted that the sample sizes and performance vary considerably depending upon which of the two arbitrary functions of the reference points are used.

IEEE VAST 2008 Challenge Data
This final example now illustrates how the method performs in the presence of a much larger number of target distributions, in the context of network security. The IEEE VAST 2008 Challenge data are synthetic but realistically generated records of mobile phone calls for a small community of 400 individuals over a ten day period. The data can be obtained from http://www.cs.umd.edu/hcil/VASTch The aim of the original challenge was to find anomalous behaviour within this social network, which might be indicative of malicious coordinated activity. One approach to this problem is to monitor the call patterns of individual users and detect any changes from their normal behaviour, with the idea that a smaller subset of anomalous individuals will then be investigated for community structure. In particular, this approach has been shown to be effective with these data when monitoring the event process of incoming call times for each individual (Heard et al., 2010). After correcting for diurnal effects on normal behaviour, this approach can be reduced to a changepoint analysis of the intensities of 400 Poisson processes of the same character as Section 5.1.2. For the focus of this article, it is of interest to see how such an approach could be made more feasible in real time by allocating computational resource between these 400 processes more efficiently.

Discussion
It was remarked in the review paper of Sisson (2005) on transdimensional samplers that "a more rigorous default assessment of sampler convergence" than the existing technology is required, and this has remained an open problem. This article is a first step towards establishing such a default method, proposing a framework and methodology which are rigorously motivated and fully general in their applicability to all distributional settings. The next step should be to explore how these ideas might extend to the common setting of Markov chain Monte Carlo sampling without thinning where the successive draws cannot be treated as approximately independent. However, this issue should not be seen as a decisive limitation of the proposed methodology, since although thinning was used in the Markov chain Monte Carlo examples of Sections 5.1.2 and 5.2 to obtain the next sample for use in calculating the convergence criteria, this would not prevent the full sample from being retained and utilised without thinning for the actual inference problem.
Another related problem which could be considered is that of importance sampling. If samples cannot be obtained directly from the target π but instead from some importance distribution with the same support, then it would be useful to understand how these error estimates and sample size strategies can be extended to the case where the empirical distribution of the samples has associated weights. In addressing the revised question of how large an importance sample should be, there should be an interesting trade-off between the inherent complexity of the target distributions, which has been the subject of this article, and how well the importance distributions match those targets.