text:
# Using Hierarchical Data Mining to Characterize Performance of Wireless System Configurations

## 1 Introduction

Data mining is becoming increasingly relevant in simulation methodology and computational science \[RG01\]. It entails the ‘non-trivial process of identifying valid, novel, potentially useful, and ultimately understandable patterns in data’ \[FPSS96\]. Data mining can be used in both predictive (e.g., quantitative assessment of factors on some performance metric) and descriptive (e.g., summarization and system characterization) settings. Our goal in this paper is to demonstrate a hierarchical data mining framework applied to the problem of characterizing wireless system performance.

This work is done in the context of the S4W problem solving environment \[VHW02\]—‘Site-Specific System Simulator for Wireless System Design’. S4W provides site-specific (deterministic) electromagnetic propagation models as well as stochastic wireless system models for predicting the performance of wireless systems in specific environments, such as office buildings. S4W is also designed to support the inclusion of new models into the system, visualization of results produced by the models, integration of optimization loops around the models, validation of models by comparison with field measurements, and management of the results produced by a large series of experiments. In this paper, we study the effect of configuration parameters on the bit error probability (BEP) of a system simulated in S4W.

The approach we take is to accumulate a performance database of simulation runs that sweep over a targeted space of system configurations. This database is then mined to obtain regions of the configuration space that exhibit acceptable average performance. Exploiting prior knowledge about the underlying simulation, organizing the computational steps in data mining, and interpreting the results at every stage, are important research issues. In addition, we bring out the often prevailing tension between making statistically meaningful conclusions and the assumptions required for efficient and effective data mining algorithms. This interplay leads to a novel set of problems that we address in the context of the wireless systems performance domain.

Data mining algorithms work in a variety of ways but, for the purposes of this paper, it is helpful to think of them as performing systematic aggregation and redescription of data into higher-level objects. Our work can be viewed as employing three such layers of aggregation: points, buckets, and regions. Points (configurations) are records in the performance database. These records contain configuration parameters as well as unbiased estimates of bit error probabilities that we use as performance metrics. Buckets represent averages of points. We use buckets to reduce data dimensionality to two, which is the most convenient number of dimensions for visualization. Finally, buckets are aggregated into 2D regions of constrained shape. We find regions of buckets where we are most confident that the configurations exhibit acceptable average performance. The shapes of these regions illustrate the nature of the joint influence of the two selected configuration parameters on the configuration performance. Specific region attributes, such as region width, provide estimates for the thresholds of sensitivity of configurations to variations in parameter values.

### 1.1 Reader’s Guide

Our major contribution is the development of a statistical framework for assessing wireless system performance using data mining techniques. The following section outlines wireless systems performance simulation methodology and develops a statistical framework for spatial aggregation of simulation results. Section 3 demonstrates a substantial subset of this framework in the context of a performance study of WCDMA (wideband code division multiple access \[HT00\]) systems that employ two-branch STTD (space-time transmit diversity \[Ala98\]) techniques and 1/2 rate convolutional coding (forward error correction codes \[HT00\]). We study the effect of power imbalance between the branches on the BEP of the system across a wide range of average signal-to-noise ratios (SNRs). Section 4 extends the statistical framework to support computation of optimized regions of the bucket space. Such regions are computed by a well-known data mining algorithm \[FMMT01, YFM97\]. Section 5 applies these concepts to the example in Section 3. Section 6 summarizes present findings and outlines directions for future research.

## 2 The Statistics of Aggregation and the Aggregation of Statistics

Temporal variations in wireless channels have been extensively studied in the literature \[ZT02\]. The present work uses a Monte Carlo simulation of WCDMA wireless systems to study the effect of these variations. The simulation traces a number of frames of random information bits through the encoding filters, the channel (a Rayleigh fading linear filter \[Has93\]), and the decoding filters. The inputs are hardware parameters, average SNR, channel impulse response, and the number of frames to simulate. The output is the bit error rate—the ratio of the number of information bits decoded in error to the total number of information bits simulated. Simulations of this kind statistically model channel variations due to changes in the environment and device movement across a small geographical area (*small-scale fading* \[Has93\]). We refer to this kind of channel variation as *temporal variation* because a system is simulated over a period of time. Further, we say that a given list of inputs to the WCDMA simulation is a *configuration* or a *point* in the configuration space.

*Spatial variations* are due to changes in system configurations. We use this term to describe two quite different phenomena: changes in the average SNR and channel impulse response due to *large-scale fading* \[Has93\] and variations of hardware parameters. A typical approach to the analysis of spatial variations is to run several temporal variation simulations (i.e., compute bit error rates — BERs — at several points within a given area of interest) and plot 1D or 2D slices of the configuration space, as shown in Figure 1. In this paper, we augment this approach with statistically meaningful aggregation of performance estimates across several points. The result of this aggregation is a space of buckets, each bucket representing the aggregation of a number of points. Moving up one level of aggregation in this manner allows us to bring data mining algorithms to operate at the level of buckets. The space of buckets mined by the data mining algorithm is then visualized using color maps. The color of each bucket is the confidence that the points (configurations) that map to this bucket exhibit acceptable average performance.

### 2.1 The First Level of Aggregation: Points

Table 1 summarizes some of the syntactic conventions used in this paper. Mathematically, we can think of the WCDMA simulation as estimating the mean $`E[x_k]`$ of a random variable $`x_k`$ with some (unknown) distribution \[JBS92\] ($`x_k`$ is one when the information bit is decoded in error or zero when it is decoded correctly). Each BER $`\widehat{x}_{kj}`$, $`1jn_k`$, output by the simulation is an unbiased estimate of the BEP $`E[x_k]`$ of the simulated configuration $`c_k`$. Instead of building a detailed stochastic model of the simulation (analytically, from the distribution of $`x_k`$), we choose to work with the simpler distribution of the BER $`\widehat{x}_{kj}`$, referred to henceforth as just $`b_k`$. Thus, each sample from the distribution of $`b_k`$ is realized by simulating a number of frames and obtaining an estimate of $`E[x_k]`$. The distribution of $`b_k`$ is approximately Gaussian due to the Central Limit Theorem. Technically, we assume that the number of frames per estimate $`\widehat{x}_{kj}`$ is ‘large enough’ so that the Lindeberg condition is satisfied, that the variance of $`\widehat{x}_{kj}`$ is finite, and that $`\{\widehat{x}_{kj}\}_{j=1}^{n_k}`$ are i.i.d. We say that $`E[b_k]=E[E[x_k]]`$ is the *expected BEP* of configuration $`c_k`$ under Rayleigh fading.

### 2.2 The Second Level of Aggregation: Buckets

Let us now aggregate several points (i.e., random variables) into one bucket. The purpose of this aggregation is to reduce data dimensionality to a size that is easy to visualize, usually one or two dimensions. The basic idea is to linearly average all points that map to the same bucket but we must do so carefully, in order to preserve a meaningful statistical interpretation. Let $`\{b_k\}_{k=1}^n`$ be Gaussian random variables with means $`\{E[b_k]\}_{k=1}^n`$ and variances $`\{\sigma _k^2\}_{k=1}^n`$. As in the previous paragraph, let each such variable $`b_k`$ be the estimated BEP of some configuration $`c_k`$, $`1kn`$. For bucket $`C`$, define a *bucket* (mixture) *random variable* $`B`$ as the convex combination

$$B=\underset{k=1}{\overset{n}{}}p_kb_k,$$

where the $`p_k0`$ and $`_{k=1}^np_k=1`$. It is convenient to make $`\{p_k\}_{k=1}^n`$ the probabilities of occurrence of the configurations $`\{c_k\}_{k=1}^n`$ in the dataset being analyzed. This setup underlines the dependence of the outputs on the distribution of the inputs and frees the user from having to provide values for the constants $`\{p_k\}_{k=1}^n`$. It is well known that, as long as $`\{b_k\}_{k=1}^n`$ are *mutually independent* and Gaussian with means $`\{E[b_k]\}_{k=1}^n`$ and variances $`\{\sigma _k^2\}_{k=1}^n`$, $`B`$ is Gaussian with mean $`E[B]=_{k=1}^np_kE[b_k]`$ and variance $`\mathrm{\Sigma }^2=_{k=1}^np_k^2\sigma _k^2`$ \[CB02\]. The expected value $`E[B]`$ of the random variable $`B`$ can be viewed as the expected BEP of bucket $`C=\{c_1,c_2,\mathrm{},c_n\}`$ in a Rayleigh fading environment, conditional on the (discrete) distribution of the configurations in $`C`$.

The values $`\{p_k\}_{k=1}^n`$ are what the statisticians call *prior probabilities*. For most purposes of this paper, we simply estimate $`\{p_k\}_{k=1}^n`$ from available data. These values are explicitly or implicitly constructed during experiment design and we assume that they remain constant during experiment analysis. However, one can collect additional data as long as doing so does not change $`\{p_k\}_{k=1}^n`$. Prior probabilities can come from a number of sources: channel sounding measurements, propagation simulations, hardware and budget constraints, or even educated guesses by wireless system designers. The rest of the paper silently assumes that the values $`\{p_k\}_{k=1}^n`$ have been established beforehand. It is important to remember that even though the prior probabilities are for the most part transparent to the analysis presented here, they nonetheless always exist and all conclusions of data analysis are made conditional on the prior probabilities.

This discussion of $`\{p_k\}_{k=1}^n`$ can be interpreted as a deferral of the exact definition of $`B`$ until experiment setup, or as parameterization of the analysis procedure. A natural question is whether or not this level of parameterization is sufficient. It is sufficient for the purposes of this paper but, strictly speaking, the interrelations between $`\{b_k\}_{k=1}^n`$ should also be defined during experiment setup. Mutual independence of $`\{b_k\}_{k=1}^n`$ is a simplifying assumption and it might be desirable to model interactions between $`\{b_k\}_{k=1}^n`$ in practice. This implies adding covariance terms to $`\mathrm{\Sigma }^2`$ and re-thinking the distribution of $`B`$. Such analysis is necessarily specific to a particular experiment. For the sake of simplicity, the rest of this paper assumes mutual independence of variables in a given bucket.

### 2.3 Confidence Estimation

Point and bucket estimates of the expected BEP are meaningful performance metrics for wireless systems. Let us also estimate our confidence in these estimates. Confidence analysis enables wireless system designers to make more practical claims than point estimates alone. A statement of the form ‘this configuration will exhibit acceptable performance in 95% of the cases’ is often preferable to a statement of the form ‘the expected BEP of this configuration is approximately $`5\times 10^4`$’. More precisely, we say that configuration $`c_k`$ *exhibits acceptable performance* when the expected BEP $`E[b_k]`$ of configuration $`c_k`$ is below some fixed threshold $`T`$. This statement is conditional on the temporal simulation assumptions, i.e., Rayleigh fading. Standard values for $`T`$ are $`10^3`$ for voice quality systems and $`10^6`$ for data quality systems. Likewise, we say that bucket $`C`$ (a subspace of configurations) *exhibits acceptable average performance* when the expected BEP $`E[B]`$ of bucket $`C`$ is below some fixed threshold $`T`$. This statement is conditional on both the temporal simulation assumptions and the distribution of configurations $`\{c_k\}_{k=1}^n`$ in the bucket (the prior probabilities).

The confidence that configuration $`c_k`$ (resp. bucket $`C`$) exhibits acceptable (average) performance is $`P(E[b_k]<T)`$ (resp. $`P(E[B]<T)`$). Since $`b_k`$ and $`B`$ are Gaussian, these probabilities can be estimated as

$$P(E[b_k]<T)F_{n_k1}\left(\frac{T\widehat{b}_k}{\widehat{\sigma }_k/\sqrt{n_k}}\right),P(E[B]<T)F_{N1}\left(\frac{T\widehat{B}}{\widehat{\mathrm{\Sigma }}/\sqrt{N}}\right),$$

where $`F_{N1}()`$ is the CDF of the Student $`t`$ distribution with $`N1`$ degrees of freedom and $`n_k`$ and $`N`$ are the sample sizes for configuration $`c_k`$ and bucket $`C`$, respectively. For configuration $`c_k`$,

$$\widehat{b}_k=\frac{1}{n_k}\underset{j=1}{\overset{n_k}{}}\widehat{x}_{kj},\widehat{\sigma }_k^2=\frac{1}{(n_k1)}\underset{j=1}{\overset{n_k}{}}(\widehat{x}_{kj}\widehat{b}_k)^2,$$

where $`\widehat{b}_k`$ and $`\widehat{\sigma }_k^2`$ are the estimates of the expected BEP and the BEP variance at point $`c_k`$, $`n_k2`$ is sample size, and $`\{\widehat{x}_{kj}\}_{j=1}^{n_k}`$ are sample values. For bucket $`C`$, we substitute point estimates into $`E[B]=_{k=1}^np_kE[b_k]`$ and $`\mathrm{\Sigma }^2=_{k=1}^np_k^2\sigma _k^2`$ to obtain

$$\widehat{B}=\underset{k=1}{\overset{n}{}}\widehat{p}_k\widehat{b}_k,\widehat{\mathrm{\Sigma }}^2=\underset{k=1}{\overset{n}{}}\widehat{p}_k^2\widehat{\sigma }_k^2,$$

where $`\widehat{B}`$ and $`\widehat{\mathrm{\Sigma }}^2`$ are the estimates of the expected BEP and the BEP variance at bucket $`C`$, and $`\{\widehat{p}_k\}_{k=1}^n`$ are the prior probabilities estimated from the dataset as $`\widehat{p}_k=n_k/_{i=1}^nn_i`$. Observe that

$$\widehat{B}=\underset{k=1}{\overset{n}{}}\widehat{p}_k\widehat{b}_k=\frac{1}{_{k=1}^nn_k}\underset{k=1}{\overset{n}{}}\underset{j=1}{\overset{n_k}{}}\widehat{x}_{kj}$$

is exactly the sample mean of all observations in the bucket, but $`\widehat{\mathrm{\Sigma }}^2`$ is *not* the variance of this sample. This is the case because $`\{\{\widehat{x}_{kj}\}_{j=1}^{n_k}\}_{k=1}^n`$ are not i.i.d. samples from the mixture distribution of $`B`$—they are samples from the constituent distributions of $`\{b_k\}_{k=1}^n`$.

## 3 Extended Example

Let us now apply the techniques developed so far to analyze the performance of a space of configurations. The wireless systems under consideration employ WCDMA technology with two-branch STTD and $`1/2`$ rate convolutional coding. We require that the transmitter has two antennas (branches) separated by a distance large enough for their signals to be uncorrelated, but small enough for the mean path losses and impulse responses of their channels to be approximately equal at receiver locations of interest. We assume Rayleigh flat fading channels, which is reasonable for indoor applications in the ISM and UNII carrier frequency bands (2.4 and 5.2 GHz, respectively). The goal is to study the effect of power imbalance between the branches on the BEP of the configurations across a wide range of average SNRs.

This section presents a number of plots that summarize simulated BERs. We also outline the process of statistically significant sampling of the configuration space. The next section develops a data mining methodology that solves a practically important problem: given a dataset similar to the one presented next, find a region of the configuration space where we can confidently claim that configurations will exhibit acceptable (average) performance.

Let us begin with an initial sample of the configuration space, as shown in Figure 2 (top). This figure shows the simulated BER as a 2D function $`\widehat{f}(S_1,S_2)`$ of the average branch bit energy-to-noise ratios (SNRs) $`S_1`$ and $`S_2`$, in dB. The parallel simulation ran for three days on 120 machines (AMD Athlon 1.0 GHz) at a speed of approximately 2.5 points per machine per day. 10000 frames, or 800000 information bits, were simulated for each of the 820 points $`S_2=3,4,\mathrm{},42`$; $`S_1=3,4,\mathrm{},S_2`$. Since $`\widehat{f}(S_1,S_2)`$ is symmetric \[DDJS02\], we show $`M=1600`$ points $`\{c_k\}_{k=1}^M`$ for a full cross-product of $`S_1`$ and $`S_2`$.

Wireless system designers are more accustomed to 1D slices of the configuration space, e.g., the ones shown in Figure 3. Define the *branch power imbalance factor*

$$\alpha =10^{0.1|S_1S_2|},$$

where $`S_1`$ and $`S_2`$ are the average SNRs of the branches, in dB. (This definition applies as long as the mean path losses of the branches are equal.) By definition, $`0\alpha 1`$, where zero corresponds to a total malfunction of one of the branches and one corresponds to a perfect balance of branch powers. The graphs in Figure 3 were obtained by fixing $`\alpha `$ and varying the *effective SNR*

$$S=10\mathrm{log}_{10}\left((10^{0.1S_1}+10^{0.1S_2})/2\right),$$

in dB (top), and fixing the effective SNR and varying $`\alpha `$ (bottom). (Note that fixing the effective SNR is equivalent to fixing total transmitter power.) The sample of configurations came from the dataset shown in Figure 2 (bottom), described in detail later. However, this sample does not contain the exact points for typical slices, so we used a fitted surface—Figure 4—to approximate the BERs for the slices in Figure 3. We choose to work with the axes $`S_1,S_2`$ in Figure 2 because it simplifies the discussion later.

What can be gathered from Figure 2 (top)? The deep valley along the diagonal is due to the fact that, provided that the effective SNR is fixed, we expect the BEP to be smallest when the branch power is balanced ($`S_1=S_2`$, $`\alpha =1`$\[DDJS02\]. Somewhat less expected were (a) the wide *tolerance region* where $`|S_1S_2|`$ is large (up to 12 dB) but the BER is still small, (b) a very sharp decline in performance at the edge of the tolerance region, and (c) a region of high local variability in the upper part of the diagonal. The surface is truncated at

$$\underset{1kM}{\mathrm{min}}\{\widehat{b}_k\}=3.75\times 10^6$$

because smaller estimates of the (expected) BEP require an enormous computation time due to the convergence properties of Monte Carlo Estimation (more on this below).

### 3.1 Statistically Significant Sampling Methodology

The initial sample looks reasonable and uncovers interesting trends in system performance, but it does not contain enough information to make statistically significant claims. Estimating the probability that a configuration exhibits acceptable average performance requires several samples per point $`c_k`$. The simulation is computationally expensive and different regions of the configuration space exhibit different variability. Therefore, we must define tight stopping criteria for sampling. Figure 2 (bottom) shows the output obtained with the following (per point $`c_k`$) stopping criteria. The criteria are designed to achieve high estimation accuracy.

1. Sampling $`\{\widehat{x}_{kj}\}`$ stops when the relative error in the estimate $`\widehat{b}_k`$ of the expected BEP $`E[b_k]`$ is smaller than the *relative accuracy threshold* $`\beta =0.1`$ times the current estimate $`\widehat{b}_k`$, at a $`\gamma =0.9`$ confidence level, i.e., when

$$P(|E[b_k]\widehat{b}_k|<\beta \widehat{b}_k)\gamma .$$

We required $`n_k2`$ samples to obtain an estimate $`\widehat{\sigma }_k^2`$ of the BEP variance $`\sigma _k^2`$. Notice that the target is the relative error, not the absolute error, because the range of $`\{\widehat{b}_k\}_{k=1}^M`$ in the configuration space spans four orders of magnitude. Therefore, absolute error measures are misleading.
2. Sampling $`\{\widehat{x}_{kj}\}`$ also stops when we can say, with confidence $`\gamma =0.9`$, that the expected BEP $`E[b_k]`$ is below the *sampling threshold* $`t=10^4`$, i.e., when

$$P(E[b_k]<t)\gamma .$$

This work considers voice quality applications, so the exact value of the expected BEP is irrelevant as long as it is smaller than the performance threshold $`T=10^3`$. The sampling threshold $`t`$ was set to an order of magnitude below the performance threshold $`T`$ to avoid large approximation error of a fitted surface near $`T`$.
3. Finally, sampling $`\{\widehat{x}_{kj}\}`$ stops when more than 50 samples of 10000 frames each are required to satisfy either of the previous rules. This rule fired in 5% of the cases, all at the boundary of the tolerance region and most in mid diagonal.

Altogether, 5154 samples were collected for an average of 6.3 samples per point. Needless to say, the computational expense of such sampling remains too high for practical applications. While a large number of samples is typically desirable (for validation purposes), we will show that our data mining framework makes very effective use of data and thus requires fewer samples in practice. Let us now look at the data in more detail.

### 3.2 Results of Statistically Significant Sampling

It is also likely that the samples output by the WCDMA simulation are approximately Gaussian distributed. Intuitively, we are simulating a large number of information bits (800000) per BEP estimate $`\widehat{x}_{kj}`$, so the Lindeberg condition for the Central Limit Theorem should hold. Figure 5 shows empirical evidence that this is the case. We have arbitrarily chosen one point among those with 20–30 sample values $`\{\widehat{x}_{kj}\}_{j=1}^{n_k}`$ and plotted the empirical CDF of this sample against that of the Gaussian distribution with the mean equal to sample mean $`\widehat{b}_k`$ and the variance equal to sample variance $`\widehat{\sigma }_k^2`$. The curves are close to each other and the Shapiro-Wilk test yields $`W=0.98`$ ($`0W1`$) and $`p`$-value of $`0.88`$. Other points also demonstrate similar curves and high values of $`W`$, but $`p`$-values vary significantly. This dataset contains sufficient samples to estimate $`\{E[b_k]\}_{k=1}^M`$ with high relative accuracy, but 6.3 samples per point are insufficient to formally justify a Gaussian assumption.

It is also instructive to see some measure of how the sample variance is distributed across the configuration space. Figures 6 and 7 show sample sizes and sample standard deviation-to-mean ratios for the samples in Figure 2 (recall that we prefer relative measures because the range of $`\{\widehat{b}_k\}_{k=1}^M`$ is large). Both figures indicate high variance around the boundary of the tolerance region. This is not surprising because the edges of the tolerance region are relatively steep. Figure 7 also shows relatively high variance at some points inside the tolerance region. This is because the simulation achieved the sampling threshold $`t=10^4`$ and stopped before it achieved the relative accuracy threshold $`\beta =0.1`$. Knowing this, one would expect a larger relative variance in the tolerance region. Let us examine why this is not the case.

We treat the BEP as a continuous Gaussian random variable $`b_k`$, but all sample values $`\{\widehat{x}_{kj}\}_{j=1}^{n_k}`$ are discrete—they are ratios of two integers, the number of errors and the number of bits simulated. The simulation may not detect any bit errors when the expected BEP $`E[b_k]`$ is relatively small (e.g., one error in the number of bits simulated). Since no channel is perfect, zero is too optimistic an estimate for the expected BEP. Instead, we conservatively assume that at least three bit errors have been detected. This is why the smallest estimate $`\widehat{b}_k`$ of $`E[b_k]`$ is $`3/800000=3.75\times 10^6`$. However, using any constant cutoff prevents us from estimating the variance $`\sigma _k^2`$. We would need to simulate a large number of frames to estimate $`\sigma _k^2`$ when the expected BEP is small. Instead, we can empirically show that the probability that the expected BEP is smaller than the performance threshold $`T=10^3`$ is close to one. Let $`\widehat{b}_k=3.75\times 10^6`$ be the sample mean, $`n_k=2`$ be the sample size, and $`\sigma _k^2`$ be the BEP variance at point $`c_k`$ where two independent simulations detected three or fewer bit errors each. Sampling $`\{\widehat{x}_{kj}\}`$ will stop because sample variance is zero, so the first stopping rule applies.

We need to show that sampling can indeed stop, i.e., that the probability that the expected BEP is below the performance threshold $`T`$ is

$$P(E[b_k]<T)F_{n_k1}\left(\frac{T\widehat{b}_k}{\widehat{\sigma }_k/\sqrt{n_k}}\right)0.995.$$

This statement can only be false when $`(T\widehat{b}_k)\sqrt{n_k}/\widehat{\sigma }_k64`$, or $`\widehat{\sigma }_k2.2\times 10^5`$, almost an order of magnitude bigger than the conservative estimate $`\widehat{b}_k`$ of the expected BEP $`E[b_k]`$. This is unlikely because Figure 7 (bottom) shows that the sample standard deviation rarely exceeds the sample mean even by half an order of magnitude. In other words, we do not have accurate estimates for variance $`\sigma _k^2`$ in the tolerance region. However, we can still reasonably conclude that configurations exhibit acceptable performance in this region.

## 4 The Third Level of Aggregation: Regions

Consider a set of buckets $`\{C_k\}_{k=1}^M`$ with corresponding random variables $`\{B_k\}_{k=1}^M`$. Given a number of sample values, the framework developed in Section 2 allows us to estimate the probabilities $`\{P(E[B_k]<T)\}_{k=1}^M`$ that buckets $`\{C_k\}_{k=1}^M`$ exhibit acceptable average performance. (All arguments about buckets equally apply to points because a point is a special case of a bucket.) This section is concerned with finding an optimal subset of random variables from among $`\{B_k\}_{k=1}^M`$. This optimal subset corresponds to an optimal region of a 2D bucket space. We would like to find a sufficiently large admissible region $`R_m`$ such that we are sufficiently confident that buckets in $`R_m`$ exhibit acceptable average performance.

There are many ways to define admissibility and we are interested in adopting a definition that is both meaningful in the wireless domain and permits effective data mining algorithms. Among a space of such admissible regions, we can define different optimality criteria and data mining then reduces to searching within this space. In this paper, a region $`R_m`$ is admissible when it has a particular type of shape. We will explore three different criteria for the mining of optimal regions; the algorithms and these criteria are based on the work of Fukuda et al. \[YFM97\] and have been adapted to the problem of mining simulation data in this paper.

Additional notation relating buckets to regions is introduced in Table 2. Let $`X`$ and $`Y`$ be two discrete parameters to the temporal (e.g., WCDMA) simulations such that $`X`$ and $`Y`$ partition the point space into disjoint buckets $`\{C_k\}_{k=1}^M`$. More precisely, let $`X,Y`$ have ordinal domains $`D_X,D_Y`$, let $`|D_X|=M_X,|D_Y|=M_Y,|D_X||D_Y|=M`$, and assume that the map $`\rho :D_X\times D_Y\{C_k\}_{k=1}^M`$ is bijective. In other words, $`X`$ and $`Y`$ define a discrete 2D space of buckets. Since the domains of $`X`$ and $`Y`$ are ordinal, this space is easily visualized as a 2D color map or a 3D perspective plot.

For example, the average SNRs $`S_1`$ and $`S_2`$ in the previous section partition the space of configurations into buckets. Both $`S_1`$ and $`S_2`$ vary from 3 to 42 in steps of 1 (in dB), so $`M_X=M_Y=40`$ and $`M=40\times 40=1600`$ (recall, from Section 3, that only 820 of these points were simulated and the remaining ones were symmetrically reflected). Furthermore, the domains of $`S_1`$ and $`S_2`$ are ordinal because the values of $`S_1`$ and $`S_2`$ are directly related to the powers of the transmitter antennas. In this case, the buckets are simply the points in the space of configurations. In general, buckets can be convex combinations of points, as detailed in Section 2. Recall that we defined the color of a bucket as the probability that the bucket exhibits acceptable average performance. Figure 8 shows these ‘colors’ as a perspective plot for the STTD example.

### 4.1 Region Shape

Consider regions (subsets) of buckets in the bucket space. If the shape of these regions is unconstrained, there are $`2^M`$ possible regions $`\{R_m\}_{m=1}^{2^M}`$. Let region $`R_m`$, $`1m2^M`$, consist of buckets $`\{C_{\kappa (m,i)}\}_{i=1}^{\eta (m)}`$, where $`\eta (m)`$, $`1m2^M`$, is a mapping from region number $`m`$ to the number of buckets in this region, and $`\kappa (m,i)`$, $`1m2^M`$, $`1i\eta (m)`$, is a mapping from region number $`m`$ and bucket number $`i`$ within region $`R_m`$ to bucket number $`k`$, $`1kM`$, that we use to subscript buckets $`\{C_k\}_{k=1}^M`$. The exact definitions of $`\eta (m)`$ and $`\kappa (m,i)`$ are not important as long as they generate all possible regions (subsets) $`\{R_m\}_{m=1}^{2^M}`$.

The shape of admissible regions should be constrained because unconstrained regions are hard to interpret and tend to overfit the training data. Besides, the problem of selecting an optimal unconstrained region is computationally intractable—all $`2^M`$ possible regions must be considered, where $`M=1600`$ in the STTD example. The region shape can be constrained in a number of different ways (rectangular, x-monotone, etc.). Our restrictions on region shape are discussed next.

Without loss of generality, assume that $`D_X=\{1,2,\mathrm{},M_X\}`$ and $`D_Y=\{1,2,\mathrm{},M_Y\}`$. Intuitively, region $`R_m`$ is rectilinear when its intersection with any horizontal or vertical line is connected. More formally, region $`R_m`$ is *rectilinear* if and only if whenever buckets $`C_{\kappa (m,i)}`$ at $`(x_{\kappa (m,i)},y_{\kappa (m,i)})`$ and $`C_{\kappa (m,j)}`$ at $`(x_{\kappa (m,j)},y_{\kappa (m,j)})`$ are both in $`R_m`$, then (a) $`\rho (r,s)=C_{\kappa (m,i)}`$ and $`\rho (r,t)=C_{\kappa (m,j)}`$ imply buckets $`\rho (r,u)`$ are also in $`R_m`$ for all $`u[s,t]`$, and (b) $`\rho (r,t)=C_{\kappa (m,i)}`$ and $`\rho (s,t)=C_{\kappa (m,j)}`$ imply buckets $`\rho (u,t)`$ are also in $`R_m`$ for all $`u[r,s]`$. Here $`[a,b]`$ means all integers between the integers $`a`$, $`b`$, inclusive. We use Manhattan geometry to define connectedness. Region $`R_m`$ is *connected* if and only if for every pair of buckets $`C_{\kappa (m,i)}`$ and $`C_{\kappa (m,j)}`$ in $`R_m`$ there exists a sequence of buckets

$$C_{\kappa (m,i)}=C_{\kappa (m,l_1)},C_{\kappa (m,l_2)},\mathrm{},C_{\kappa (m,l_n)}=C_{\kappa (m,j)}$$

in $`R_m`$ such that for every $`1k<n`$

$$\rho ^1(C_{\kappa (m,l_k)})\rho ^1(C_{\kappa (m,l_{k+1})})_1=1.$$

Furthermore, we say that region $`R_m`$ is *admissible* if it is both rectilinear and connected.

This definition of admissibility can be viewed as a relaxed definition of convexity. Geometrically, it is easy to see that region $`R_m`$ is admissible if and only if, when we look at $`R_m`$ from left to right, its upper boundary first increases and then decreases monotonically (a pseudoconcave function), and its lower boundary first decreases and then increases monotonically (a pseudoconvex function). In other words, the region boundary need not be strictly convex or strictly concave, but it must be pseudoconvex or pseudoconcave. Admissible regions are informally summarized in Figure 9. All admissible regions are composed of regions of four primitive types: W (region gets wider from left to right), N (region gets narrower), U (region slants up), and D (region slants down). Twelve combinations of these types yield all types of admissible regions: W, WU, WUN, WD, WDN, WN, UN, DN, U, D, N, and the empty region.

Our choice of connected rectilinear regions is due to primarily heuristic considerations. These considerations are commonly applicable, but must be re-evaluated for each study. Both the connectedness and the rectilinearity restrictions can be justified for the STTD example (see next section). In general, it is easy to justify connectedness, but hard to justify rectilinearity. We advocate the use of connected rectilinear regions primarily because this shape is resistant to noise in the sample, not because we can analytically show that the region boundary is rectilinear. In data mining, the choice of region shape is most commonly dictated by the desired tradeoff between bias and variance \[HTF01\]. Regions with flexible shape exhibit small bias (they can fit any data) but high variance (they can be overly sensitive to a particular dataset). Regions with rigid shape exhibit high bias but small variance. Connected rectilinear regions provide a reasonable tradeoff between bias and variance for many applications.

### 4.2 Evaluating Regions

Another prerequisite to finding regions with the desired properties is a definition of region ‘goodness’. Let us map bucket confidence $`P(E[B_{\kappa (m,i)}]<T)`$ to a discrete range $`[0\mathrm{}1000]`$ and define the *hit of bucket* $`C_{\kappa (m,i)}`$ as

$$h_{\kappa (m,i)}=1000P(E[B_{\kappa (m,i)}]<T)+0.5$$

($`X`$ denotes the largest integer that does not exceed $`X`$), the *support of bucket* $`C_{\kappa (m,i)}`$ as

$$s_{\kappa (m,i)}=1000$$

(this constant was chosen to make the discretization error reasonably small), the *hit of region* $`R_m`$ as

$$H_m=\underset{i=1}{\overset{\eta (m)}{}}h_{\kappa (m,i)},$$

and the *support of region* $`R_m`$ as

$$S_m=\underset{i=1}{\overset{\eta (m)}{}}s_{\kappa (m,i)}=1000\eta (m).$$

The key to efficient computation of optimized-confidence and optimized-support admissible regions is the definition of region confidence as

$$\mathrm{\Theta }_m=H_m/S_m,$$

where $`H_m`$ is the hit and $`S_m`$ is the support of region $`R_m`$. Let us explore the implications of these definitions in more detail.

#### 4.2.1 Model-Based and Model-Free Analyses

Suppose, $`n_{\kappa (m,i)}=6`$ samples have been collected for bucket $`C_{\kappa (m,i)}`$ that consists of a single point. Let the sample mean be $`\widehat{B}_{\kappa (m,i)}=5\times 10^4`$ and the sample standard deviation be $`\widehat{\mathrm{\Sigma }}_{\kappa (m,i)}=8.87\times 10^4`$. Furthermore, suppose that five of these samples have the BER below $`10^3`$ and one has the BER above $`10^3`$. Then,

$$P(E[B_{\kappa (m,i)}]<T)F_5\left(\frac{10^35\times 10^4}{8.87\times 10^4/\sqrt{6}}\right)0.887.$$

A purely model-free approach would interpret the above simulation results as ‘bucket $`C_{\kappa (m,i)}`$ will exhibit acceptable average performance in 5 out of 6 cases.’ A strongly model-based approach would interpret the simulation results as ‘we are 88.7% confident that bucket $`C_{\kappa (m,i)}`$ exhibits acceptable average performance.’ Our interpretation lies between the model-based approach and a model-free approach and posits that ‘bucket $`C_{\kappa (m,i)}`$ will exhibit acceptable average performance in 887 out of 1000 cases.’ These interpretations provide confidence estimates under different simplifying assumptions.

The model-free interpretation does not take either sample variance or sample distribution into account. This interpretation is only reliable for a sufficiently large number of samples, which is a luxury in our application. Our middle-ground interpretation explicitly accounts for sample variance and sample distribution. When sample size is small, our interpretation provides a statistically valid estimate of confidence that the bucket exhibits acceptable average performance. For a single bucket, this interpretation is as good as a strongly model-based interpretation, modulo a reasonably small discretization error. However, our interpretation diverges from the model-based interpretation at the region level.

A strongly model-based analysis procedure would define a region random variable

$$Q_m=\frac{1}{W_m}\underset{i=1}{\overset{\eta (m)}{}}w_{\kappa (m,i)}B_{\kappa (m,i)},$$

where $`\{B_{\kappa (m,i)}\}_{i=1}^{\eta (m)}`$ are bucket random variables, $`\{w_{\kappa (m,i)}\}_{i=1}^{\eta (m)}`$ are *a priori* (positive) constant weights, and

$$W_m=\underset{i=1}{\overset{\eta (m)}{}}w_{\kappa (m,i)}$$

is a normalization factor that maps these weights to probabilities of bucket occurrence in the region. A procedure similar to that in Section 2 would then be used to estimate $`P(E[Q_m]<T)`$ for a threshold $`T`$. The result of this calculation can be interpreted as the probability that region $`R_m`$ exhibits acceptable average performance, conditional on the temporal simulation assumptions, the bucketing prior probabilities, and the region prior probabilities. However, as we shall see later, this definition of region confidence violates a property that permits an efficient data mining algorithm.

We think of region confidence in terms of average bucket confidence over the whole region, namely,

$$\mathrm{\Theta }_m\frac{1}{\eta (m)}\underset{i=1}{\overset{\eta (m)}{}}P(E[B_{\kappa (m,i)}]<T).$$

(If region size $`\eta (m)`$ is large enough, we can reasonably expect the discretization errors to cancel each other.) This interpretation of $`\mathrm{\Theta }_m`$ does not correspond to the strongly model-based probability that region $`R_m`$ exhibits acceptable average performance. Instead, we define a region random variable $`P_m`$ as the probability that *any* bucket $`C_{\kappa (m,i)}`$ in region $`R_m`$ exhibits acceptable average performance. Then, we estimate the expected value $`E[P_m]`$ across the region $`R_m`$ by the sample mean $`\widehat{P}_m\mathrm{\Theta }_m`$ of estimates of bucket confidences $`\{P(E[B_{\kappa (m,i)}]<T)\}_{i=1}^{\eta (m)}`$.

How do these two definitions relate to each other? It is easy to show that they are equivalent only under very restrictive assumptions. Basically, we are assuming that the buckets are mutually independent, that population variance is small, and that the region is consistent, i.e., ‘good’ and ‘bad’ buckets are never mixed in the same region. Let bucket random variables $`\{B_{\kappa (m,i)}\}_{i=1}^{\eta (m)}`$ be mutually independent, let the estimates $`\{\widehat{\mathrm{\Sigma }}_{\kappa (m,i)}^2\}_{i=1}^{\eta (m)}`$ of bucket variances be approximately equal to zero, and let the estimates $`\{\widehat{B}_{\kappa (m,i)}\}_{i=1}^{\eta (m)}`$ of bucket expected BEPs be either all greater than the performance threshold $`T`$ or all smaller than the performance threshold $`T`$ (i.e., all $`\{T\widehat{B}_{\kappa (m,i)}\}_{i=1}^{\eta (m)}`$ have the same sign). Then, bucket confidences

$$P(E[B_{\kappa (m,i)}]<T)F_{n_{\kappa (m,i)}1}\left(\frac{T\widehat{B}_{\kappa (m,i)}}{\widehat{\mathrm{\Sigma }}_{\kappa (m,i)}/\sqrt{n_{\kappa (m,i)}}}\right),$$

$`1i\eta (m)`$, will be either all approximately equal to zero ($`\widehat{B}_{\kappa (m,i)}>T`$), or all approximately equal to one ($`\widehat{B}_{\kappa (m,i)}<T`$). Therefore, region confidence $`\mathrm{\Theta }_m`$ will be approximately equal to zero or one. Likewise, the strongly model-based region confidence

$$P(E[Q_m]<T)F_{\eta (m)1}\left(\frac{T\widehat{Q}_m}{\widehat{\mathrm{\Psi }}_m/\sqrt{\eta (m)}}\right)$$

will be approximately equal to zero or one because the estimate $`\widehat{\mathrm{\Psi }}_m^2`$ of region variance is (see Section 2)

$$\widehat{\mathrm{\Psi }}_m^2=\frac{1}{W_m^2}\underset{i=1}{\overset{\eta (m)}{}}w_{\kappa (m,i)}^2\widehat{\mathrm{\Sigma }}_{\kappa (m,i)}^20.$$

The sign of $`T\widehat{Q}_m`$ determines whether $`P(E[Q_m]<T)`$ is approximately equal to zero or one. After a minor rearrangement of terms,

$$T\widehat{Q}_m=\frac{1}{W_m}\underset{i=1}{\overset{\eta (m)}{}}w_{\kappa (m,i)}(T\widehat{B}_{\kappa (m,i)}).$$

We assumed that $`\{T\widehat{B}_{\kappa (m,i)}\}_{i=1}^{\eta (m)}`$ have the same sign, so we have shown that $`P(E[Q_m]<T)\mathrm{\Theta }_m`$. The equality is asymptotically exact as all variance estimates $`\{\widehat{\mathrm{\Sigma }}_{\kappa (m,i)}^2\}_{i=1}^{\eta (m)}`$ approach zero. This argument applies regardless of the distributions of $`\{B_{\kappa (m,i)}\}_{i=1}^{\eta (m)}`$, as long as these random variables are mutually independent.

### 4.3 Optimized Regions

We now pursue the definition of optimized regions. Given a slope $`\tau `$, $`0\tau 1`$, define the *gain* of region $`R_m`$, $`1m2^M`$, as

$$G(R_m,\tau )=H_m\tau S_m,$$

where $`H_m`$ is the region hit and $`S_m`$ is the region support. Let an *optimized-gain admissible region* $`R_\tau `$ with respect to slope $`\tau `$, $`0\tau 1`$, be an admissible region with the maximum gain $`G(R_\tau ,\tau )`$ over all admissible regions (this region need not be unique). Optimized-gain admissible regions are easy to define, compute, and analyze, but hard to interpret. Common practice is to define optimized-confidence and optimized-support admissible regions. Admissible region $`R_{}`$ is an *optimized-confidence admissible region* with respect to a given support threshold $`1000\eta `$, $`0\eta M`$, if $`R_{}`$ has the maximum confidence $`\mathrm{\Theta }_{}=H_{}/S_{}`$ among all admissible regions with support of at least $`1000\eta `$. Likewise, admissible region $`R_{}`$ is an *optimized-support admissible region* with respect to a given confidence threshold $`\theta `$, $`0\theta 1`$, if $`R_{}`$ has the maximum support $`S_{}=1000\eta ()`$ among all admissible regions with confidence of at least $`\theta `$. In other words, we can either fix the region confidence $`\theta `$ and find the largest region $`R_{}`$ with confidence of at least $`\theta `$, or we can fix the minimum region size (support) $`1000\eta `$ and find the most confident region $`R_{}`$ with support of at least $`1000\eta `$.

Observe that $`\tau `$ in the definition of an optimized-gain admissible region is the relative importance of support vs. that of confidence. We can find a small region with high confidence or a large region with small confidence, but both objectives cannot be maximized simultaneously. Increasing $`\tau `$ will increase the confidence of the optimized-gain admissible region, but decrease its support. Likewise, decreasing $`\tau `$ will decrease the confidence of the optimized-gain admissible region, but increase its support. Therefore, we can find approximate optimized-confidence and optimized-support admissible regions by a binary search for the value of $`\tau `$ where the respective threshold is barely satisfied. The search can stop at a given level of precision $`\mathrm{\Delta }\tau `$, where the lower bound on $`\mathrm{\Delta }\tau `$ can be found in \[YFM97\] (they show that the number of steps in this search is logarithmic in the support $`1000M`$ of the bucket space). This algorithm is approximate because an optimized-confidence (resp. optimized-support) admissible region need not be an optimized-gain admissible region for any value of $`\tau `$. Yoda et al. \[YFM97\] argue that this approximation is reasonable for large datasets.

Let us revisit the definition of region ‘goodness’. Geometrically, the buckets with the same value of $`X`$ are the columns and the buckets with the same value of $`Y`$ are the rows. An optimized-gain admissible region can be computed in $`O(M_XM_Y^2)`$ time by a set of rules of the following form. Recall that a region of type W gets wider from left to right (see Figure 9). Let $`R_W(m,[s,t])`$ be the region of type W with maximum gain $`f_W(m,[s,t])`$ over all admissible regions of type W that end in column $`m`$ and span rows $`s`$ through $`t`$ in this column. Then, either (a) $`m`$ is the first column of $`R_W(m,[s,t])`$, or (b) $`R_W(m,[s,t])`$ includes the region $`R_W(m1,[s^{},t^{}])`$ with the maximum gain $`f_W(m,[s^{},t^{}])`$ over all admissible regions that end in column $`m1`$ and span rows $`s^{}s`$ through $`t^{}t`$ in this column. \[YFM97\] keeps the regions with maximum gain for every region type and every triple $`(m,[s,t])`$ in a dynamic programming table. These locally maximal regions then grow according to a set of rules that compute an optimized-gain admissible region. This efficient greedy algorithm for computing optimized-gain admissible regions depends on the property of the gain function that we refer to as monotonicity. Let $`0\tau 1`$ be a slope and $`R_m^{}`$ and $`R_{m^{\prime \prime }}`$ be two admissible regions with gains $`G(R_m^{},\tau )G(R_{m^{\prime \prime }},\tau )`$. The gain function $`G(R_m,\tau )`$ is *monotonic* if for any region $`R_k`$ disjoint with both $`R_m^{}`$ and $`R_{m^{\prime \prime }}`$

$$G(R_m^{}R_k,\tau )G(R_{m^{\prime \prime }}R_k,\tau ),$$

where the union of regions is defined in the obvious way. It is easy to see that our gain function

$$G(R_m,\tau )=H_m\tau S_m=\underset{i=1}{\overset{\eta (m)}{}}1000P(E[B_{\kappa (m,i)}]<T)+0.51000\tau \eta (m)$$

is monotonic because it is additive. However, a strongly model-based gain function

$$G^{(M)}(R_m,\tau )=P(E[Q_m]<T)\tau \eta (m)/M$$

is not monotonic even if we assume independence of bucket random variables $`\{B_{\kappa (m,i)}\}_{i=1}^{\eta (m)}`$ that make up $`Q_m`$. To the best of our knowledge, only monotonic gain functions are known to result in practical algorithms for computing optimized-gain admissible regions.

What happens when no estimates of mean and/or variance are available for some bucket $`C_k`$? The answer to this question depends on problem-specific considerations. As was demonstrated in Section 3, it is sometimes possible to provide conservative estimates for these values. For example, we have empirically shown that the expected BEPs of some configurations $`\{c_k\}`$ are smaller than $`T=10^3`$ with confidence $`P(E[b_k]<T)0.995`$. Likewise, we know that as the effective SNR approaches negative infinity (in dB), the BEP approaches 0.5, which is the probability of correctly guessing the value of a random bit when the transmitter is turned off. Thus, we can let $`P(E[b_k]<T)=0`$ for points with sufficiently small effective SNRs and a reasonable performance threshold $`T`$. If no such estimates are available, we can simply omit the missing buckets from the probability computation. This must be done with care because such buckets will contribute nothing to the confidence of the region. This fact can be used to reduce the computational expense of sampling.

This section has highlighted the sometimes contradictory objectives that aggregation must satisfy: permit valid statistical interpretations and afford structure that can be exploited by data mining algorithms. Our approach has been a judicious mix of concepts from both statistics and data mining. We showed that our formulation of the data mining problem lies between the completely model-free approach and the strongly model-based approach. The next section applies the data mining methodology described here to the example in Section 3.

## 5 Optimized-Support Regions for the STTD Example

This section continues the example in Section 3. First, we show that optimized-gain regions are both rectilinear and connected for this example. It immediately follows that optimized-support and optimized-confidence regions are also admissible. An optimized-support admissible region is presented next. We show that the elaborate region mining setup leads to simple engineering interpretations. Finally, we look at the performance of data mining when the number of samples is small. Three-fold cross-validation shows that data mining performs well under these circumstances.

### 5.1 Justification of Data Mining for the STTD Example

Let the average SNRs $`S_1=X`$ and $`S_2=Y`$ partition the space of configurations in Figure 2 into disjoint points (buckets) $`\{c_k\}_{k=1}^M`$, $`1M1600`$. We now give an intuitive argument to justify the suitability of the data mining algorithm for the STTD study. Without loss of generality, consider only the points with $`XY`$, i.e., $`S_1S_2`$. It is easy to extend all arguments to $`X>Y`$, but this adds little to the discussion.

Let $`c_1`$ at $`(x_1,y_1)`$ and $`c_2`$ at $`(x_2,y_1)`$, $`x_1<x_2<y_1`$, be two points in an optimized-gain region (of arbitrary shape) for some slope $`0<\tau <1`$ (see Figure 10). This means that the confidences of these points are one, and thus the expected BEPs of these points are smaller than the performance threshold $`T`$. When $`x_1,x_2<y_1`$ and $`y_1`$ is fixed, the BEP is a monotonically decreasing function of $`x`$—increasing $`x`$ decreases the power imbalance and increases the effective SNR, so the BEP must decrease. Therefore, the expected BEP of any point $`c_u`$ at $`(x_u,y_1)`$, $`x_1<x_u<x_2`$, is below the performance threshold $`T`$. Thus, the confidences of points $`\{c_u\}`$ are one and these points must also be in the optimized-gain region $`R_\tau `$. Three more symmetric arguments of this kind show that optimized-gain regions are rectilinear.

Likewise, let $`c_1`$ at $`(x_1,y_1)`$ and $`c_2`$ at $`(x_2,y_2)`$, $`x_1<x_2<y_1<y_2`$, be two points in an optimized-gain rectilinear region (refer to Figure 10). Since $`c_1`$ is in the optimized-gain region and $`x_1<x_2`$, the point at $`(x_2,y_1)`$ is also in this region because it has a smaller BEP than $`c_1`$. Since the optimized-gain region is rectilinear, there is a horizontal path from $`(x_1,y_1)`$ to $`(x_2,y_1)`$ and a vertical path from $`(x_2,y_1)`$ to $`(x_2,y_2)`$. Thus, there is a Manhattan path from $`(x_1,y_1)`$ to $`(x_2,y_2)`$. Arguments of this kind show that optimized-gain rectilinear regions must be connected as long as they are ‘wide enough’.

To summarize, we have shown that optimized-gain (and thus optimized-support and optimized-confidence) regions are admissible. The data mining algorithm described in Section 4, which results in optimal admissible regions, is thus appropriate for the STTD example. We now show and interpret data mining results.

### 5.2 Optimized-Support Admissible Regions

Figure 11 shows an optimized-support admissible region for the confidence threshold $`\theta =0.99`$. Intuitively, this is the largest admissible region where we can claim, with confidence of at least $`0.99`$, that configurations exhibit acceptable performance. This claim is conditional on temporal simulation assumptions and on mutual independence of configurations in the region. The shape of this region confirms that, under a fixed effective SNR, the BEP is minimal when the average SNRs of the two branches are equal. The width of this region shows the largest acceptable power imbalance. For this example, the system tolerates power imbalance of up to 12 dB. However, the width of the optimized region is not uniform. The region is narrower for small effective SNRs and wider for large effective SNRs. This means that configurations with low effective SNRs are more sensitive to power imbalance than configurations with high effective SNRs. None of these observations are news to an informed reader. The contribution of data mining in this context is not qualitative discoveries; it is statistically significant quantitative results.

Let us see how the data mining algorithm performs when data is scarce. The initial sample of the configuration space in Figure 2 (top) contains one sample value per bucket. The statistically significant sample in Figure 2 (bottom) contains at least two additional sample values per bucket (recall that we required at least two sample values to estimate bucket variance $`\sigma _k^2`$). Therefore, three-fold cross-validation is the most elaborate cross-validation procedure that this dataset affords. The regions in the top left, bottom left, and top right of Figure 12 have been computed for sample values in Figure 2 (top) and the first two sample values per bucket in Figure 2 (bottom). Each of these regions has been computed with two out of the three sample values per bucket. The region in the lower right has been computed with all data in Figure 2 (bottom). All four regions are optimized-support admissible regions with the confidence threshold $`\theta =0.95`$. The regions are overlaid on top of the color-coded bucket confidence values. Red (dark) corresponds to low confidence and white (light) corresponds to high confidence that configuration $`c_k`$ exhibits acceptable average performance w.r.t. the voice quality threshold $`T=10^3`$.

The regions in Figure 12 are identical except for the lower left corner. This is not surprising because this part of the configuration space exhibits high relative variance. Also, the data is symmetric but the regions are asymmetric in the lower-left corner. Recall that optimized-gain admissible regions, and thus optimized-support admissible regions, are not unique. The ties in region gains are broken arbitrarily. Therefore, region asymmetry is an additional indicator of region instability.

Figure 12 also shows that additional data improves image contrast but does not significantly affect region shape. Collecting additional sample values separates the points into ones with low confidence and ones with high confidence. A curious side effect occurs when the difference in confidence estimates of low-confidence points falls below the discretization error (1/1000). In this case, the ‘confidence slack’ $`1\theta `$ is allocated to arbitrary points with low confidence. One way to correct this situation is to raise the confidence threshold $`\theta `$—after all, more accurate data should afford stronger claims. Another alternative is to lower the discretization threshold. In general, optimized regions work best when the data is noisy. A contour plot will suffice when the data is highly accurate.

It can also be seen that the high contrast created by the sharp edge of the tolerance region is advantageous to data mining. The region is stable where the contrast is high. When the image is blurred, data mining tries to avoid the questionable boundary points.

To summarize, this section has demonstrated that optimized-gain regions are rectilinear and connected for a non-trivial space of wireless system configurations. We have also shown that optimized-support admissible regions are easy to interpret. Finally, we have shown that data mining works well when sample sizes are small.

## 6 Discussion and Future Work

We have demonstrated a hierarchical formulation of data mining suitable for assessing performance of wireless system configurations. WCDMA simulation results are systematically aggregated and redescribed, leading to intuitive regions that allow the engineer to evaluate wireless system configuration parameters. We have shown that the assumptions about region shape and properties made by data mining algorithms can be valid in the wireless design context; the patterns mined hence lead to explainable and statistically valid design conclusions. As a methodology, data mining is thus shown to be extremely powerful when coupled with statistically meaningful performance evaluation.

This work is the first (known to the authors) application of data mining methodology to solve problems in wireless system design. Therefore, a large number of extensions are possible and called for. We outline possible extensions at the three levels of aggregation: points, buckets, and regions.

At the point level, it may be advantageous to model temporal simulations more precisely. This paper assumes a ‘large enough’ number of frames per simulation and works with the distribution of estimated BEPs. We have shown reasonable analytical and empirical evidence that this distribution is Gaussian. The advantage of this problem formulation is the independence of spatial aggregation from the assumptions of temporal simulation. This helps introduce wireless engineers to the methodology of data mining for studying design problems. However, a stronger model of temporal simulation (e.g., Markov chains in \[Wan95\]) may yield appreciable gains in software performance. This direction is worth pursuing because few research groups have access to parallel computing facilities of the scale used in this work. For instance, the initial sample of the configuration space in Figure 2 (top) would take one year of computation time on a modern workstation. The study presented in this paper would clearly be impossible without significant computational power.

Aggregation of points into buckets is the least developed part of this work. Suppose that we would like to simulate the effects of interference on configuration performance. Assume that the distribution of the average strengths of the interfering signals is known *a priori* (e.g., estimated by ray tracing). We can either make this distribution known to the temporal simulation, or, alternatively, run several temporal simulations for different strengths of interfering signals. The former is more accurate and computationally more efficient, but the latter is more generic and simpler to implement. Bucketing of simulation results with varying simulation parameters is intended to approximate the performance of a single device under varying conditions. This paper does not employ such bucketing but instead builds all the necessary kinds of parameter variation into the temporal simulation (which can be argued to be the right way to do it). However, bucketing may be necessary when one has to work with a given dataset (e.g., measurements). Bucket space can be viewed as a configuration space for a more complex temporal simulation. Therefore, an in-depth treatment of bucketing is orthogonal to the primary topic of this paper, which is data mining.

Significant work remains to be done at the region level as well. For instance, the assumption of small variance could conceivably be relaxed. One can also pursue the relatively difficult task of incorporating strongly model-based prior knowledge into the data mining algorithm, or the somewhat easier task of applying different kinds of region mining algorithms to problems in wireless system design.

Defining additional case studies is another obvious direction for future work. We have studied a relatively small part of the parameter space of modern wireless systems. More studies of this type must be performed to highlight the merits and the shortcomings of data mining in this domain.

Finally, the strict staging of data collection and data mining can be relaxed. One can fruitfully interleave the two activities and have the results of data mining drive subsequent data collection. In data-scarce domains, it would be advantageous to focus the data collection effort on only those regions deemed most important to support a particular data mining objective. Methodologies for closing-the-loop in this manner are becoming increasingly prevalent \[RBK02\]. This will also help define alternative criteria for evaluating experiment designs and layouts.


format: text

text:
# Spin reorientation transition in the incommensurate stripe-ordered phase of La3/2Sr1/2NiO4

## Abstract

The spin ordering of La3/2Sr1/2NiO4 was investigated by magnetization measurements, and by unpolarized- and polarized-neutron diffraction. Spin ordering with an incommensurability $`\epsilon `$ $``$ 0.445 is observed below $`T_{\mathrm{SO}}`$ $``$ 80 K. On cooling, a spin reorientation is observed at 57 $`\pm `$ 1 K, with the spin axes rotating from 52 $`\pm `$ 4 to 78 $`\pm `$ 3. This is the first time a spin reorientation has been observed in a La2-xSrxNiO4+δ compound having incommensurate stripe order.

preprint: ATB-1

Incommensurate magnetic fluctuations observed in superconducting La2-xSrxCuO4(LSCO)super  are widely believed to be associated with charge stripe correlations.tranquada-Nature-1995  Assuming this to be the case, it has been concluded that charge stripes in LSCO lie parallel to the Cu–O bonds in the CuO planes for $`x`$ 0.053,fujita-PRB-2002  but are aligned at $`45^{}`$ to the Cu–O bond for $`x`$ 0.07. fujita-PRB-2002 ; superstripe  Given that the critical doping level for superconductivivty in LSCO is $`x`$ = 0.053–0.056,fujita-PRB-2002  these observations suggest a strong link between spin-charge stripe ordering and high temperature superconductivity.

Diagonal charge stripes have also been observed in La2-xSrxNiO4+δ(LSNO),tranquada-PRL-1993  and have been studied by neutronneutron ; yoshizawa-PRB-2000 ; lee-PRB-2001 ; kajimoto-PRB-2001 ; kajimoto  and x-rayx-ray ; pash-PRL-2000  diffraction for doping levels in the range 0.135 $`x`$ 0.5. The charge stripes are found to be ordered over a long range, with correlation lengths in excess of 100 Å for certain levels of dopingyoshizawa-PRB-2000 ; pash-PRL-2000 . LSNO with $`x`$ = 1/3 has stripe order that is particularly stable owing to a commensurability effect that pins the charges to the lattice.ramirez-PRL-1996 ; yoshizawa-PRB-2000 ; kajimoto-PRB-2001  Having charge stripes that are static and well correlated, LSNO is a good system for studying the basic properties of spin-charge stripes.

Recently, Lee et al.lee-PRB-2001  measured the direction of the ordered Ni moments in the spin-ordered phase of compounds with $`x`$ = 0.275 and 1/3. The angle between the spin axes and the stripe direction was found to be larger at $`x`$ = 1/3 than $`x`$ = 0.275. Furthermore, a transition was observed at $``$50 K in the $`x`$ = 1/3 sample such that on cooling the spins rotate away from the stripe direction and a further localization of the charge ordering occurs. Lee et al. concluded from this that an intimate connection exists between the local spin order and the stability of the charge order.

Early work on La3/2Sr1/2NiO4+δ by electron diffractionchen-PRL-1993  suggested that the case of half-doping is another very stable charge-ordering phase. This charge ordering has a 2-D checkerboard-like pattern as shown in figure 1(a). Recently, Kajimoto et al.kajimoto  studied the spin and charge order in a single crystal of La3/2Sr1/2NiO4 by neutron diffraction. They observed checkerboard ordering below $``$480 K, but they also found an incommensurate charge-stripe ordering phase below $`T_{\mathrm{ICO}}`$$``$ 180 K coexisting with the checkerboard ordering. Spin ordering was reported to be found only in association with the incommensurate component of the charge ordering.

Here we describe a study of spin-ordering in La3/2Sr1/2NiO4+δ employing magnetometry and polarized- and unpolarized-neutron diffraction. Our results are consistent with those of Kajimoto et al.,kajimoto  but in addition we observed a spin reorientation transition similar to, but more pronounced than that observed in La5/3Sr1/3NiO4. This is the first time a spin reorientation associated with an incommensurate spin ordering has been observed in a LSNO material. The spin reorientation is signalled by a sharp anomaly in the magnetization at $``$57 K, and an analysis of polarized-neutron diffraction data shows that on cooling through 57 K the Ni spins rotate away from the stripe direction.

Two single crystals of La3/2Sr1/2NiO4+δ were used for these experiments. The crystal used for the magnetization measurements had dimensions $``$$`7\times 5\times 2`$ mm3 and that for neutron diffraction was a rod of 7–8 mm diameter and $``$35 mm in length. Both these samples were grown from the same starting powder by the floating-zone technique.Prab  The oxygen excess was determined by thermogravimetric analysis to be $`\delta `$ = 0.02 $`\pm `$ 0.01.

Magnetization data were collected with a SQUID magnetometer (Quantum Design) with the field applied parallel to the $`ab`$ plane of the crystal. We carried out DC measurements by cooling in an applied field of 500 Oe (FC) and also by cooling in zero field and measuring while warming in a field of 500 Oe (ZFC). Figure 2 shows the temperature variation of the FC and ZFC magnetization. A small change in slope at 90 $`\pm `$ 5 K marks the spin ordering temperature, in agreement with previous observations.yoshizawa-PRB-2000  A clearer anomaly is seen at $`T_{\mathrm{SR}}`$ = 57 K, where the FC and ZFC curves sharply separate. Below 5 K both the FC and ZFC curves begin to rise for reasons that are not yet understood, but could be due to a small amount of paramagnetic impurity in the crystal. We measured the FC magnetization parallel to the crystal $`c`$ direction (not shown) and found no anomaly around $`T_{\mathrm{SR}}`$, thus indicating that the transition involves the in-plane components of the spins.

In what follows we describe neutron diffraction studies aimed at identifying the nature of the transition at $`T_{\mathrm{SR}}`$. Referred to the tetragonal unit cell of LSNO, with unit cell parameters $`a`$ $``$ 3.8 Å, $`c`$ $``$ 12.7 Å, charge ordering peaks due to structural distortions are located at ($`h\pm \epsilon `$, $`k\pm \epsilon `$, $`l`$) positions in reciprocal space. The spin order peaks are located at ($`h`$+1/2$`\pm \epsilon `$/2, $`k`$+1/2$`\pm \epsilon `$/2, $`l`$), the strongest being when $`l`$ is odd.lee-PRB-2001  To a first approximation the relationship between the incommensurability, $`\epsilon `$, and the hole concentration, $`n_h=x`$ \+ 2 $`\delta `$, is simply $`\epsilon `$ = $`n_h`$. In general, however, competition between apparently different wavevectors for the spin and charge orderingkajimoto-PRB-2001  and the separate issue of oxygen non-stoichiometryjest-PRB-1999 ; wochner-PRB-1998 ; chen-PRL-1993  causes small deviations from the $`\epsilon =n_h`$ rule.

The neutron experiments were performed on the triple-axis spectrometer IN20 at the Institut Laue-Langevin. The energies of the incident and elastically-scattered neutrons were selected by Bragg reflection from an array of either pyrolytic graphite (PG) or Heusler alloy crystals, for unpolarized- and polarized-neutron scattering respectively. A PG filter was present after the sample and before the analyzer to suppress higher-order harmonic scattering. Most of the data were obtained with initial and final neutron wavevectors of 2.66 Å-1, with some further data taken with a neutron wavevector of 3.54 Å-1. We mounted the sample with the  and  crystal directions in the horizontal scattering plane, and performed scans in reciprocal space either parallel to the ($`h`$$`h`$, 0) direction at constant $`l`$, or parallel to the (0, 0, $`l`$) direction at constant $`h`$.

We began by performing polarized-neutron diffraction measurements with the neutron polarization P parallel to the scattering vector Q, maintained by adjusting currents in a Helmholtz coil mounted around the sample. In this configuration a neutron’s spin is flipped during a magnetic scattering process, but remains unchanged when scattered by a non-magnetic process, e.g. from a distortion of the lattice. Thus, by measuring the spin- flip (SF) and non-spin-flip (NSF) channels one can identify whether the origin of the scattering is magnetic or not.

We observed charge-order diffraction in the vicinity of the (1.5, 1.5, $`l`$) line in reciprocal space. At $`T`$ = 10 K the peak shape observed in ($`h`$$`h`$, 0) scans was consistent with that observed by Kajimoto et al., with a central peak at $`h`$ = 1.5 and two satellites at $`h`$ 1.45 and 1.55. We observed no variation of this scattering with $`l`$ over a single Brillouin zone, showing that both the commensurate and incommensurate components of charge ordering are completely two-dimensional.

We observed spin-order diffraction peaks at ($`h`$+1/2$`\pm \epsilon `$/2, $`h`$+1/2$`\pm \epsilon `$/2, $`l`$) positions with $`l`$ odd and $`\epsilon `$ = 0.445 $`\pm `$ 0.005. Figure 3 shows a scan parallel to (0, 0, $`l`$). The widths of the peaks in this scan convert to a correlation length along the $`c`$ axis of 16.4 $`\pm `$ 0.3 Å. For comparison, the in-plane correlation length in the direction perpendicular to the stripes is 78.3 $`\pm `$ 1.3 Å. There were no peaks at $`l`$ = even positions, as can be seen in figure 3. There were also no spin-order peaks at ($`h\pm `$1/4, $`h\pm `$1/4, l) positions, confirming the lack of any commensurate Fourier component in the spin ordering, in agreement with the observations of Kajimoto et al..kajimoto

Figure 4a show the temperature dependence of the integrated intensity of two magnetic reflections, Q1 = (0.275,0.275,5) and Q2 = (0725,0.725,1), measured in scans parallel to $`(h,h,0)`$. The spin ordering transition is seen to be rather sluggish, but a sharpening of the peaks below $`T_{\mathrm{SO}}`$ $``$ 80 K (not shown) suggests that this is where long-range ordering extending over many unit cells sets in. However short-range correlations persist above 100 K. On cooling below $`T_{\mathrm{SO}}`$, the peak at Q1 increases monotonically, while that at Q2 first increases, then reaches a maximum at $``$60 K, then decreases again. This anomalous behaviour correlates with the transition observed in the magnetization, Fig. 2. To understand this behaviour we must recall that magnetic neutron diffraction is sensitive to spin components perpendicular to Q. As shown in figure 1(b), Q1 is aproximately parallel to $`(0,0,l)`$ while Q2 is approximately parallel to $`(h,h,0)`$. Therefore, to a good approximation, the scattering at Q1 is from the total in-plane spin moment, irrespective of its direction, while that at Q2 comes mainly from the spin components parallel to the stripes and along the $`c`$ axis. The data in fig. 4a therefore imply that below $``$57 K the spins rotate away from the stripe direction.

To fully analyze the direction of the spins over this temperature range we varied the direction of the neutron polarization P relative to the scattering vector Q. Three configuration were used: 1) P $``$ Q, 2) P $``$ Q but in the scattering plane, and 3) P $``$ Q but out of the scattering plane. (Recall that  and  are the crystal directions in the scattering plane.) As mentioned before, neutrons scatter via magnetic interactions from spin components perpendicular to Q, and SF scattering is due to spin components perpendicular to P. Therefore, analysis of configurations 1), 2) and 3) leads to the direction of the ordered moment.111Actually, 2 out of the 3 configurations are sufficient. Table 1 summarizes the relations between the observed intensities at Q1 and Q2 and the ordered spin components parallel ($`\mu `$) and perpendicular ($`\mu `$) to the stripe direction (see figure 1(a)), and parallel to the $`c`$ axis ($`\mu _c`$) .

We corrected the data to take into account the different background count rates in the SF and the NSF channels, and to correct for the imperfect spin polarization of neutron beam. The latter was calculated from the flipping ratio 18.5 $`\pm `$ 1.8 measured on a magnetic Bragg peak.

A full analysis of the intensities at the Q1 and Q2 points determined the spin orientations at $`T`$ = 10 K and $`T`$ = 65 K. These analyses showed that the non-magnetic (NM) and $`\mu _c`$ count rates were zero, to within error for both temperatures. Thus, it is reasonable to assume that the spins lie within the $`ab`$ plane over the whole temperature range, consistent with the lack of any anomaly observed at $`T_{\mathrm{SR}}`$ in the magnetization data measured with $`𝐇c`$. At the other temperatures we deduced the spin orientation within the $`ab`$ plane from one Q position ($`T`$= 55 K from Q1, $`T`$ = 50 and 75 K from Q2.) From these data we find that the spins reorientate from an angle of 52 $`\pm `$ 4 to the stripe direction at $`T`$ = 75 K, to 78 $`\pm `$ 3 at $`T`$ = 10 K, as shown in figure 4(b). The transition is seen to occur quite rapidly between 50 and 65 K, although the reorientation is not complete at 50 K.

How does the spin reorientation described here in the $`x`$ = 1/2 crystal compare with that reported by Lee et al.lee-PRB-2001  for $`x`$ = 1/3? In both compositions the spin reorientations occur at similar temperatures ($``$50 K for $`x`$ = 1/3, $``$57 K for $`x`$ = 1/2) and the spins rotate in the same sense, namely away from the stripe direction, on cooling. On the other hand, the angle through which the spins rotate is greater for $`x`$ = 1/2 ($``$26) than for $`x`$ = 1/3 ($``$13), and the final low temperature spin direction in $`x`$ = 1/2 is almost perpendicular to the stripes ($`\varphi `$ = 78) compared with 53 for $`x`$ = 1/3. Considering also the spin direction $`\varphi `$ = 27 observed for $`x`$ = 0.275,lee-PRB-2001  for which there is no reorientation transition, we find that as the doping level increases the angle from the stripe direction increases. However, the trend in the spin orientation also correlates with the increasing stability of the charge ordering as measured by the charge-ordering temperaturesyoshizawa-PRB-2000 ; lee-PRB-2001 ; kajimoto  ($`T_{\mathrm{CO}}`$($`x`$ = 0.275) $`200`$ K, $`T_{\mathrm{CO}}`$($`x`$ = 1/3) $`240`$ K, $`T_{\mathrm{CO}}`$($`x`$ = 1/2) $`480`$ K). Therefore, measurements on more compositions of La2-xSrxNiO4+δ are needed to establish whether it is the doping level or the stability of the charge order that is more important in determining the spin orientation.

The occurrence of spin reorientation transitions in the $`x`$ = 1/3 and $`x`$ = 1/2 but not in $`x`$ = 0.275 might suggest that these transitions are linked to the stability of the charge order. Irrespective of this, what the present study has shown is that the occurrence of a spin reorientation transition does not depend on the stripe order being commensurate with lattice. Indeed, the spin reorientation transition is even more striking in the $`x`$ = 1/2 composition, for which the spin order is incommensurate, than in $`x`$ = 1/3, for which it is commensurate.

These observations emphasize the importance of couplings between spin and charge degrees of freedom in stripe-ordered systems, and suggest that this coupling is particularly strong in La3/2Sr1/2NiO4.

This work was supported in part by the Engineering and Physical Sciences Research Council of Great Britain.


text:
# REFERENCES

\[

\]

Aji and Goldenfeld Reply: In a recent letter , we analyzed the critical dynamics of the superconducting to normal phase transition in zero magnetic field. We explained Monte Carlo (MC) simulation results in both strong and weak screening limits, where the dynamic critical exponent was found to be $`z_{MC}2.7`$ and $`z_{MC}1.5`$ respectively. These results, taken at face value were surprising, departing strongly from scaling expectations based on model A dynamics that $`z2`$. We showed that the simulations do not measure the true dynamic exponent $`z`$ and that in both the short and long range limits, provided the screening length is smaller than the system size, the dynamic exponent correctly inferred from the simulations is $`z2`$, thus removing the discrepancy.

In the preceding Comment, the author asserts that in our theory, the vorticity has the wrong scaling dimension, inconsistent with standard scaling theory and not supported by MC data. Second, he asserts, without justification, that the identification of MC time with real time is correct, and that our proposal that they are not the same is an influence of the discreteness of the lattice.

In fact, the non-equivalence of MC time and real time is a genuine effect, but arises due to an implementation of a dynamic MC algorithm that does not properly account for the scaling of $`J`$. Here, we show that our theory predicts a scale dependent coupling constant $`J(L)`$, and that the standard finite size scaling results are indeed satisfied, contrary to the statement in the Comment.

In terms of the vorticity, $`\stackrel{}{n}=\stackrel{}{}\times \stackrel{}{}\theta `$, the free energy is, $`F=\beta _{i,j}\stackrel{}{n_i}.\stackrel{}{n_j}G_{ij}[\lambda _0],`$ where the lattice Green function is

$$G_{ij}[\lambda _0]=J\frac{(2\pi )^2}{L^3}\underset{\stackrel{}{k}}{}\frac{\mathrm{exp}[i\stackrel{}{k}(\stackrel{}{r_i}\stackrel{}{r_j})]}{2_m^3[1\mathrm{cos}(k_m)]+\lambda _0^2}$$
(1)
Here $`J`$ is the coupling constant, $`\lambda _0`$ is the screening length and $`\beta =1/k_BT`$. In the long range case, $`\lambda _0\mathrm{}`$, the $`G_{ij}J/|\stackrel{}{r_i}\stackrel{}{r_j}|`$, while in the short range case, $`G_{ij}J\delta (\stackrel{}{r_i}\stackrel{}{r_j})`$. The free energy density scales as $`L^d`$ ($`F(L/\xi )^d`$), $`\theta `$, the phase, has a trivial scaling dimension of $`0`$ and $`n\xi ^2`$. In the long range case, since $`G_{ij}J/\xi `$, we get $`J\xi ^2/L^3`$ where $`\xi `$ is the correlation length ($`(L/\xi )^dL^{2d}\xi ^4J/\xi `$). In the short range case, where $`G_{ij}J/L^d`$, a similar analysis yields $`J\xi `$. In the dynamic MC simulations, the scaling of $`J`$ is not accounted for and leads to scaling dimensions of $`n`$ other than $`2`$. If one ignored the scaling of $`J`$, the above analysis would yield $`n\xi ^x`$ and $`\theta \xi ^{1x}`$, where $`x`$ is $`5/2`$ ($`(L/\xi )^dL^{2d}\xi ^{2x}/\xi `$) in the weak screening limit and $`x=3/2`$ in the strong screening limit. The superfluid density is obtained from, $`F=\beta 𝑑\stackrel{}{r}\rho |\stackrel{}{}\theta |^2.`$ Given the scaling dimension of $`\theta `$ above, and assuming that $`J`$ does not scale, $`\rho J\xi ^{2x2d}`$ ($`(L/\xi )^d\rho L^d\xi ^{22x}`$), which would disagree with standard scaling theory .

This discrepancy in the scaling dimensions arises only in MC simulations where a time dependent variable is measured, such as a two-time correlation function. The superfluid density in static MC simulations, for example, is obtained from the equal time current-current correlation function in the phase representation of the Villain model. Unlike the vorticity, the phase variable does not scale with system size ($`\theta L^0`$). This computation of the correlation function reproduces the standard scaling form for $`\rho `$.

The non-equivalence of real time and MC time arises from using the same value of $`J`$ for all lattice sizes $`L`$ , and as explained in our letter, is a consequence of the fact that in a single MC time step, regardless of $`L`$, an equal change must be made to the voltage pulse produced by the vortex loops and not the vortex loops themselves. In dynamic MC simulations, if we scaled $`J`$ appropriately, we would obtain $`z2`$, for strongly screened interactions. However, one usually does not know ahead of time what the scale dependence of $`J(L)`$ is.

In our theory, the scaling form for $`\rho `$ is calculated from $`\rho J(L)L^{2xd2}`$. In the case of short ranged interaction and the physically relevant case of weak screening ($`\lambda _0<L`$), $`x=3/2`$ and $`JL`$ ($`\xi L`$ at $`T_c`$), we get $`\rho L^{2d}`$ as expected. In the long range case, where $`\lambda _0`$ is set to infinity, the same argument yields $`\rho L^{2d}`$ since $`x=5/2`$ and $`JL^1`$ in this case. For the magnetic permeability, including the scaling of $`J`$ implies, $`\mu J^1L^{d2x}`$. With this correction, we obtain the correct scaling form, $`\mu L^{d4}`$, in both the strong and weak screening limit.

In summary, our analysis is fully consistent with standard results, and explains the surprising results of the MC simulations.

We acknowledge National Science Foundation grant NSF-DMR-99-70690.

Vivek Aji and Nigel Goldenfeld
University of Illinois at Urbana-Champaign
Department of Physics
1110 West Green Street
Urbana, Illinois 61801-3080