# Comparing multiple frequency distributions from data and from simulation

I am attempting to analyze the shuffling in the card based online game Magic Arena, which I suspect has a specific bug that causes the distribution of outputs to be biased. Is the approach I have designed appropriate, and am I executing it correctly?

What I have:

• Usable data from an estimated 100000 to 300000 games (I haven't run the aggregations to find out yet), gathered client-side by a mechanism I am reasonably certain has no relevant self-selection bias.
• A specific hypothesis for exactly what I think is wrong.
• A personal implementation of the hypothesized bug, which I can run arbitrarily many times to generate reference statistics.

What I don't have:

• Direct access to the actual game's shuffling. It is done by closed source software on Wizards of the Coast company servers, and only sends any output when someone actually plays the game.
• A theoretical derivation of the exact distribution that the hypothesized bug produces. I believe producing one would be technically possible but far too difficult and computationally expensive to be feasible.

A typical deck in this game has 60 cards, around 24 of which are of a type called lands, and those lands are of especially relevant significance. The bug's effect is dependent on position in the deck prior to shuffling, rather than card type, however. I am analyzing how many of the first ~24 cards in the deck are drawn in the opening hand of 7 cards, and separately how many of the last ~24 cards, which my simulation produced very different distributions for.

My plan:

1. Run my implementation of the bug one billion times, recording frequency distributions for the first and last 24 cards before shuffling showing up in the first 7 after shuffling.
2. Use a chi-square two sample test to compute two p-values - one for the first 24 cards in the deck, the other for the last 24, in both cases comparing data from the game vs data from my simulation. As I understand it, I need the two sample variation rather than the more common Pearson's version because my predicted distribution is itself generated by a random sample rather than derived theoretically.
3. Use Fisher's method to combine these p-values into one.
4. Compare the result with the chosen significance level of 0.05.

I have not yet aggregated the raw data into useful statistics, but I have run through the technique with made up data for an example.

Suppose I had data from 10000 games for the number of the first 24 cards that got shuffled into the opening hand, distributed as in the table below:

$$\begin{array} {|l|r|r|r|r|r|r|r|r|} \hline &\text{0 in hand}&\text{1 in hand}&\text{2 in hand}&\text{3 in hand}&\text{4 in hand}&\text{5 in hand}&\text{6 in hand}&\text{7 in hand}\\ \hline \text{Example}&102&793&1965&3190&2423&1266&241&20\\ \hline \text{Simulation}&9337686&68686547&201678618&306136866&259242574&122298732&29747381&2871596\\ \hline \end{array}$$

Following the instructions I have for the chi-square two sample test, I calculate the scaling constants:

$$K_1=\sqrt{1000000000/10000}=100\sqrt{10}\approx316.2278$$ $$K_2=\sqrt{10000/1000000000}=\sqrt{10}\div1000\approx0.003162278$$

Showing the whole calculation just for 0 in hand:

$$\frac{(K_1R_0-K_2S_0)^2}{R_0+S_0}=\frac{(100\sqrt{10}\times102-\sqrt{10}\div1000\times9337686)^2}{102+9337686}=\frac{10(10200-9337.686)^2}{9337788}\approx0.79632$$

Computing the sum for all columns, I get:

$$\chi^2=\sum_{i=0}^7\frac{(K_1R_i-K_2S_i)^2}{R_i+S_i}\approx49.8825$$

I have 8 bins and unequal sample sizes, so 8 degrees of freedom. The corresponding p-value according to Wolfram Alpha (assuming I interpreted what to enter correctly) is about $$4.3049\times10^{-8}$$.

Suppose when I go through this for the distributions of the last 24 cards I get a p-value of 0.67.

Fisher's method gives a combined test statistic of $$\chi^2=-2(\ln(4.3049\times10^{-8})+\ln(0.67))\approx34.723$$

I put this into Wolfram Alpha with 4 degrees of freedom (twice the number of p-values being combined), getting an overall p-value of $$5.2959487\times10^{-7}$$

As this is less than 0.05, I would reject my hypothesis for this example.

Is all of this correct, in both choice of testing methods and their execution, and if not what mistake(s) did I make?