| 1 | package de.uka.ipd.sdq.statistics; |
| 2 | |
| 3 | import java.util.ArrayList; |
| 4 | import java.util.LinkedList; |
| 5 | import java.util.List; |
| 6 | |
| 7 | import org.apache.log4j.BasicConfigurator; |
| 8 | import org.apache.log4j.Level; |
| 9 | import org.apache.log4j.Logger; |
| 10 | |
| 11 | import de.uka.ipd.sdq.statistics.independence.IIndependenceTest; |
| 12 | import de.uka.ipd.sdq.statistics.independence.RunUpTest; |
| 13 | |
| 14 | /** |
| 15 | * Implements a batch means procedure based on phi-mixing conditions as |
| 16 | * described in [1]. Appropriate batch sizes and the number of batches are |
| 17 | * determined automatically. |
| 18 | * <p> |
| 19 | * The procedure utilizes an independence test in order to build a so-called |
| 20 | * "quasi independent" (QI) sample sequence. By default the {@link RunUpTest} |
| 21 | * will be used. "The aim of the QI method is to continue the simulation run |
| 22 | * until we have obtained a pre-specified number of essentially independent |
| 23 | * random samples by skipping highly correlated observations." [1] As soon as |
| 24 | * the QI sequence appears to be independent, the computed batches can be |
| 25 | * considered as valid. Samples in the QI sequence are only used to determine |
| 26 | * appropriate batch sizes. They are not used to compute the batch means! |
| 27 | * Instead, the batch means consist of all samples, regardless of statistical |
| 28 | * dependence. |
| 29 | * <p> |
| 30 | * [1] E. Chen, W. Kelton: A Stopping Procedure based on Phi-Mixing Conditions. |
| 31 | * Proceedings of the 2000 Winter Simulation Conference. |
| 32 | * |
| 33 | * @author Philipp Merkle |
| 34 | * |
| 35 | */ |
| 36 | public class PhiMixingBatchAlgorithm extends ABatchAlgorithm { |
| 37 | |
| 38 | private Logger logger; |
| 39 | |
| 40 | private IIndependenceTest independenceTest; |
| 41 | |
| 42 | private QuasiIndependentSampleSequence qiSamples; |
| 43 | |
| 44 | /** iteration number k */ |
| 45 | private int k = 0; |
| 46 | |
| 47 | /** |
| 48 | * number of quasi independent samples needed until the next iteration can |
| 49 | * start |
| 50 | */ |
| 51 | private int demandedQiSamples; |
| 52 | |
| 53 | /** total number of samples observed */ |
| 54 | private int totalSamples; |
| 55 | |
| 56 | private State state; |
| 57 | |
| 58 | public PhiMixingBatchAlgorithm() { |
| 59 | this(new RunUpTest()); |
| 60 | } |
| 61 | |
| 62 | public PhiMixingBatchAlgorithm(IIndependenceTest test) { |
| 63 | this.independenceTest = test; |
| 64 | |
| 65 | initialize(); |
| 66 | } |
| 67 | |
| 68 | private void initialize() { |
| 69 | logger = Logger.getLogger("de.uka.ipd.sdq.statistics.PhiMixingBatchAlgorithm.log"); |
| 70 | |
| 71 | qiSamples = new QuasiIndependentSampleSequence(); |
| 72 | |
| 73 | state = new Iteration0(); |
| 74 | } |
| 75 | |
| 76 | @Override |
| 77 | public void offerSample(double value) { |
| 78 | state.offerSample(value); |
| 79 | } |
| 80 | |
| 81 | private void setState(State state) { |
| 82 | this.state = state; |
| 83 | } |
| 84 | |
| 85 | /** |
| 86 | * Merges adjacent batches by taking their averages. |
| 87 | * |
| 88 | * @param batches |
| 89 | */ |
| 90 | private void mergeAdjacentBatches(List<Batch> batches) { |
| 91 | assert (batches.size() % 2 == 0); |
| 92 | for (int i = 0; i < batches.size() / 2; i++) { |
| 93 | Double mergedSum = batches.get(2 * i).getSum() |
| 94 | + batches.get(2 * i + 1).getSum(); |
| 95 | Integer mergedCount = batches.get(2 * i).getSize() |
| 96 | + batches.get(2 * i + 1).getSize(); |
| 97 | batches.set(i, new Batch(mergedSum, mergedCount)); |
| 98 | } |
| 99 | batches.subList(batches.size() / 2, batches.size()).clear(); |
| 100 | } |
| 101 | |
| 102 | /** |
| 103 | * Computes batch means from the buffered data. Only applicable if the quasi |
| 104 | * independent sample sequence has lag 1 since otherwise the sequence does |
| 105 | * not represent all seen samples but a subset. |
| 106 | * |
| 107 | * @param number |
| 108 | * the number of batches to compute |
| 109 | */ |
| 110 | private void computeBatches(int number) { |
| 111 | // TODO Ensure buffer.size() / batches is integer |
| 112 | assert (qiSamples.size() % number == 0); |
| 113 | int batchSize = qiSamples.size() / number; |
| 114 | |
| 115 | for (int i = 0; i < number; i++) { |
| 116 | double batchSum = 0.0; |
| 117 | int batchCount = 0; |
| 118 | for (int j = 0; j < batchSize; j++) { |
| 119 | int idx = i * batchSize + j; |
| 120 | batchSum += qiSamples.get(idx); |
| 121 | ++batchCount; |
| 122 | } |
| 123 | |
| 124 | batches.add(new Batch(batchSum, batchCount)); |
| 125 | } |
| 126 | } |
| 127 | |
| 128 | private abstract class State { |
| 129 | |
| 130 | private static final int INITIALIZE = 1; |
| 131 | |
| 132 | private static final int READY = 2; |
| 133 | |
| 134 | private static final int FINISH = 1; |
| 135 | |
| 136 | private int internalState = INITIALIZE; |
| 137 | |
| 138 | abstract void start(); |
| 139 | |
| 140 | public void offerSample(double sample) { |
| 141 | if (internalState == INITIALIZE) { |
| 142 | start(); |
| 143 | internalState = READY; |
| 144 | } |
| 145 | totalSamples++; |
| 146 | } |
| 147 | |
| 148 | abstract void end(); |
| 149 | |
| 150 | protected void changeState(State state) { |
| 151 | internalState = FINISH; |
| 152 | end(); |
| 153 | setState(state); |
| 154 | } |
| 155 | |
| 156 | protected void independenceConditionalStateChange(int lag, |
| 157 | State independent, State notIndependent) { |
| 158 | if (independenceTest.testIndependence(qiSamples.lagSublist(lag))) { |
| 159 | logger.info("Independence test passed. Batch means are valid."); |
| 160 | if (batches.isEmpty()) |
| 161 | computeBatches(4); |
| 162 | setValid(true); |
| 163 | changeState(independent); |
| 164 | } else { |
| 165 | logger.info("Independence test not passed. Need more samples."); |
| 166 | changeState(notIndependent); |
| 167 | } |
| 168 | } |
| 169 | |
| 170 | protected void logBatchStatus() { |
| 171 | if (batches.isEmpty()) { |
| 172 | logger.info("No batches determinined so far. Total samples " |
| 173 | + "seen: " + totalSamples); |
| 174 | } else { |
| 175 | logger.info("There are " + batches.size() + " batches of size " |
| 176 | + batches.get(0).getSize() + ". Total samples seen: " |
| 177 | + totalSamples); |
| 178 | } |
| 179 | } |
| 180 | |
| 181 | } |
| 182 | |
| 183 | private class Iteration0 extends State { |
| 184 | |
| 185 | public void start() { |
| 186 | k = 0; |
| 187 | logger.info("Start of iteration 0."); |
| 188 | demandedQiSamples += 4000; |
| 189 | } |
| 190 | |
| 191 | @Override |
| 192 | public void offerSample(double sample) { |
| 193 | super.offerSample(sample); |
| 194 | qiSamples.offerSample(sample); |
| 195 | if (demandedQiSamples == 0) { |
| 196 | super.independenceConditionalStateChange(1, |
| 197 | new IncreaseBatches(), new Iteration1a()); |
| 198 | } |
| 199 | } |
| 200 | |
| 201 | public void end() { |
| 202 | logger.info("End of iteration 0."); |
| 203 | super.logBatchStatus(); |
| 204 | } |
| 205 | |
| 206 | } |
| 207 | |
| 208 | private class Iteration1a extends State { |
| 209 | |
| 210 | public void start() { |
| 211 | logger.info("Start of iteration 1a"); |
| 212 | k = 1; |
| 213 | demandedQiSamples += 4000; |
| 214 | } |
| 215 | |
| 216 | @Override |
| 217 | public void offerSample(double sample) { |
| 218 | super.offerSample(sample); |
| 219 | qiSamples.offerSample(sample); |
| 220 | if (demandedQiSamples == 0) { |
| 221 | super.independenceConditionalStateChange(2, |
| 222 | new IncreaseBatches(), new Iteration1b()); |
| 223 | } |
| 224 | } |
| 225 | |
| 226 | public void end() { |
| 227 | logger.info("End of iteration 1a."); |
| 228 | super.logBatchStatus(); |
| 229 | } |
| 230 | |
| 231 | } |
| 232 | |
| 233 | private class Iteration1b extends State { |
| 234 | |
| 235 | public void start() { |
| 236 | logger.info("Start of iteration 1b"); |
| 237 | demandedQiSamples += 4000; |
| 238 | } |
| 239 | |
| 240 | @Override |
| 241 | public void offerSample(double sample) { |
| 242 | super.offerSample(sample); |
| 243 | qiSamples.offerSample(sample); |
| 244 | if (demandedQiSamples == 0) { |
| 245 | super.independenceConditionalStateChange(3, |
| 246 | new IncreaseBatches(), new IterationKa()); |
| 247 | } |
| 248 | } |
| 249 | |
| 250 | public void end() { |
| 251 | logger.info("End of iteration 1b"); |
| 252 | super.logBatchStatus(); |
| 253 | } |
| 254 | |
| 255 | } |
| 256 | |
| 257 | private class IterationKa extends State { |
| 258 | |
| 259 | public void start() { |
| 260 | ++k; |
| 261 | logger.info("Start of iteration " + k + "a."); |
| 262 | if (batches.isEmpty()) { |
| 263 | computeBatches(3); |
| 264 | } else { |
| 265 | assert (batches.size() == 3); |
| 266 | } |
| 267 | // prepare 4th batch mean |
| 268 | batches.add(new Batch()); |
| 269 | // discard every 2nd sample |
| 270 | qiSamples.discardAndRearrange(); |
| 271 | // obtain 2000 quasi independent samples... |
| 272 | demandedQiSamples = 2000; |
| 273 | // ...having lag 2^(k-1) |
| 274 | qiSamples.setLag((int) Math.pow(2.0, k - 1)); |
| 275 | } |
| 276 | |
| 277 | @Override |
| 278 | public void offerSample(double sample) { |
| 279 | super.offerSample(sample); |
| 280 | qiSamples.offerSample(sample); |
| 281 | |
| 282 | // batch means for iterations 0, 1(a) and 1(b) are calculated |
| 283 | // directly from the buffer. For iterations k >= 2 we do no longer |
| 284 | // store every observation in buffer because of a given lag. Hence |
| 285 | // we need to adjust the current batch mean every observation. |
| 286 | batches.get(batches.size() - 1).addSample(sample); |
| 287 | |
| 288 | if (demandedQiSamples == 0) { |
| 289 | super.independenceConditionalStateChange(2, |
| 290 | new IncreaseBatches(), new IterationKb()); |
| 291 | } |
| 292 | } |
| 293 | |
| 294 | public void end() { |
| 295 | logger.info("End of iteration " + k + "a."); |
| 296 | super.logBatchStatus(); |
| 297 | } |
| 298 | |
| 299 | } |
| 300 | |
| 301 | private class IterationKb extends State { |
| 302 | |
| 303 | public void start() { |
| 304 | logger.info("Start of iteration " + k + "b."); |
| 305 | assert (batches.size() == 4); |
| 306 | // reduce the number of batches from 4 to 2 by |
| 307 | // taking average of adjacent batch means |
| 308 | mergeAdjacentBatches(batches); |
| 309 | // obtain 4000 quasi independent samples... |
| 310 | demandedQiSamples = 4000; |
| 311 | // ...having lag 2^(k-1) |
| 312 | qiSamples.setLag((int) Math.pow(2.0, k - 1)); |
| 313 | // prepare 3rd batch mean |
| 314 | batches.add(new Batch()); |
| 315 | } |
| 316 | |
| 317 | @Override |
| 318 | public void offerSample(double sample) { |
| 319 | super.offerSample(sample); |
| 320 | qiSamples.offerSample(sample); |
| 321 | |
| 322 | // batch means for iterations 0, 1(a) and 1(b) are calculated |
| 323 | // directly from the buffer. For iterations k >= 2 we do no longer |
| 324 | // store every observation in buffer because of a given lag. Hence |
| 325 | // we need to adjust the current batch mean every observation. |
| 326 | batches.get(batches.size() - 1).addSample(sample); |
| 327 | |
| 328 | if (demandedQiSamples == 0) { |
| 329 | super.independenceConditionalStateChange(3, |
| 330 | new IncreaseBatches(), new IterationKa()); |
| 331 | } |
| 332 | } |
| 333 | |
| 334 | public void end() { |
| 335 | logger.info("End of iteration " + k + "b."); |
| 336 | super.logBatchStatus(); |
| 337 | } |
| 338 | |
| 339 | } |
| 340 | |
| 341 | private class IncreaseBatches extends State { |
| 342 | |
| 343 | private int demandedSamples; |
| 344 | |
| 345 | public void start() { |
| 346 | logger.info("Start increasing batches."); |
| 347 | |
| 348 | if (batches.size() == 30) { |
| 349 | mergeAdjacentBatches(batches); |
| 350 | } |
| 351 | |
| 352 | batches.add(new Batch()); |
| 353 | int batchSize = batches.get(0).getSize(); |
| 354 | demandedSamples = batchSize; |
| 355 | } |
| 356 | |
| 357 | @Override |
| 358 | public void offerSample(double sample) { |
| 359 | super.offerSample(sample); |
| 360 | batches.get(batches.size() - 1).addSample(sample); |
| 361 | --demandedSamples; |
| 362 | |
| 363 | if (demandedSamples == 0) { |
| 364 | changeState(new IncreaseBatches()); |
| 365 | } |
| 366 | } |
| 367 | |
| 368 | public void end() { |
| 369 | logger.info("Stop increasing batches."); |
| 370 | super.logBatchStatus(); |
| 371 | } |
| 372 | |
| 373 | } |
| 374 | |
| 375 | private class QuasiIndependentSampleSequence { |
| 376 | |
| 377 | /** |
| 378 | * buffer to store the quasi independent samples (i.e. a subset of the |
| 379 | * observation sequence) |
| 380 | */ |
| 381 | private List<Double> buffer; |
| 382 | |
| 383 | /** number of forthcoming observations to ignore until a sample is taken */ |
| 384 | private int lag; |
| 385 | |
| 386 | /** number of observations ignored */ |
| 387 | private int ignoredSamples; |
| 388 | |
| 389 | public QuasiIndependentSampleSequence() { |
| 390 | buffer = new LinkedList<Double>(); |
| 391 | lag = 1; |
| 392 | } |
| 393 | |
| 394 | public void setLag(int lag) { |
| 395 | assert (lag >= 1); |
| 396 | this.lag = lag; |
| 397 | } |
| 398 | |
| 399 | public void offerSample(double sample) { |
| 400 | if (ignoredSamples == lag - 1) { |
| 401 | buffer.add(sample); |
| 402 | |
| 403 | --demandedQiSamples; |
| 404 | ignoredSamples = 0; |
| 405 | } else { |
| 406 | ++ignoredSamples; |
| 407 | assert (ignoredSamples <= ignoredSamples); |
| 408 | } |
| 409 | } |
| 410 | |
| 411 | public double get(int index) { |
| 412 | return buffer.get(index); |
| 413 | } |
| 414 | |
| 415 | public int size() { |
| 416 | return buffer.size(); |
| 417 | } |
| 418 | |
| 419 | public List<Double> lagSublist(int lag) { |
| 420 | ArrayList<Double> res = new ArrayList<Double>(); |
| 421 | for (int i = 0; i < buffer.size(); i += lag) { |
| 422 | res.add(buffer.get(i)); |
| 423 | } |
| 424 | return res; |
| 425 | } |
| 426 | |
| 427 | /** |
| 428 | * Discards every 2nd sample and rearranges the remaining samples in the |
| 429 | * first half of the buffer. |
| 430 | * |
| 431 | * @param buffer |
| 432 | */ |
| 433 | public void discardAndRearrange() { |
| 434 | for (int i = 0; i < buffer.size() / 2; i++) { |
| 435 | buffer.set(i, buffer.get(2 * i)); |
| 436 | } |
| 437 | buffer.subList(buffer.size() / 2, buffer.size()).clear(); |
| 438 | } |
| 439 | |
| 440 | } |
| 441 | |
| 442 | } |