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 | } |