# Coins pooling or not
Student ID: F74082125
Student Name: 張議隆
More information in [hackmd website](https://hackmd.io/@bebo1010/coin_pooling)
:::spoiler TOC
[TOC]
:::
## How to run my code
Most commands are already written in Makefile

SEPERATE_CODE=Seperate_experiment
COMBINE_CODE=Combined_experiment

setup:
        export GSL_RNG_SEED=74082125
run_seperate:
        gcc -o $(SEPERATE_CODE)  $(SEPERATE_CODE).c  GSLfun.c -lgsl -lgslcblas -lm
        ./$(SEPERATE_CODE)

run_combine:
        gcc -o $(COMBINE_CODE)  $(COMBINE_CODE).c  GSLfun.c -lgsl -lgslcblas -lm
        ./$(COMBINE_CODE)

clean:
        rm $(SEPERATE_CODE)
        rm $(COMBINE_CODE)

To set the RNG seed, run `make setup` or `make` in short.
This command may not always work(at least in my enviorment). If `make setup` did not set the RNG seed, run `export GSL_RNG_SEED=74082125` manually to set the seed.
Then you could run `make run_seperate` and `run_combine` to get the simulation result.
If the simulation took too long to generate the result, modify line 23 to a smaller number.

const unsigned long long int sample_count = 1000 * (1 << sample_power);

## Seperate trials
### Possible outcomes of 4 trials(s1)
Table entry$(s < 4)$ taken directly from given pdf file
$C^{i - \frac{1}{2}}_{i} = \dfrac{(i-\frac{1}{2})(i-\frac{3}{2}) \cdots \frac{1}{2}}{i!}$
| s | 0 | 1 | 2 | 3 | 4 |
| -------- | -------- | -------- | -------- | -------- | -------- |
| $C^{s - 0.5}_{s}$ | 1 | $\frac{1}{2}$ | $\frac{3}{8}$ | $\frac{5}{16}$ | $\frac{35}{128}$ |
| $P[s]$ | $\frac{35}{128}$ | $\frac{5}{32}$ | $\frac{9}{64}$ | $\frac{5}{32}$ | $\frac{35}{128}$ |
| $128*P[s]$ | $35$ | $20$ | $18$ | $20$ | $35$ |

### Possible outcomes of 5 trials(s2)
Extended from the table above, and re-calculate $P[s]$

| s | 0 | 1 | 2 | 3 | 4 | 5 |
| -------- | -------- | -------- | -------- | -------- | -------- | -------- |
| $C^{s - 0.5}_{s}$ | 1 | $\frac{1}{2}$ | $\frac{3}{8}$ | $\frac{5}{16}$ | $\frac{35}{128}$ | $\frac{63}{256}$ |
| $P[s]$ | $\frac{63}{256}$ | $\frac{35}{256}$ | $\frac{15}{128}$ | $\frac{15}{128}$ | $\frac{35}{256}$ | $\frac{63}{256}$ |
| $256*P[s]$ | $63$ | $35$ | $30$ | $30$ | $35$ | $63$ |

### Probability table of s1, s2
| | $s_1 = 0$ | 1 | 2 | 3 | 4 |
| --------- | ---------- | - | - | - | - |
| $s_2 = 0$ | $\frac{63}{256} * \frac{35}{128}$ | $\frac{63}{256} * \frac{5}{32}$ | $\frac{63}{256} * \frac{9}{64}$ | $\frac{63}{256} * \frac{5}{32}$ | $\frac{63}{256} * \frac{35}{128}$ |
| 1 | $\frac{35}{256} * \frac{35}{128}$ | $\frac{35}{256} * \frac{5}{32}$ | $\frac{35}{256} * \frac{9}{64}$ | $\frac{35}{256} * \frac{5}{32}$ | $\frac{35}{256} * \frac{35}{128}$ |
| 2 | $\frac{15}{128} * \frac{35}{128}$ | $\frac{15}{128} * \frac{5}{32}$ | $\frac{15}{128} * \frac{9}{64}$ | $\frac{15}{128} * \frac{5}{32}$ | $\frac{15}{128} * \frac{35}{128}$ |
| 3 | $\frac{15}{128} * \frac{35}{128}$ | $\frac{15}{128} * \frac{5}{32}$ | $\frac{15}{128} * \frac{9}{64}$ | $\frac{15}{128} * \frac{5}{32}$ | $\frac{15}{128} * \frac{35}{128}$ |
| 4 | $\frac{35}{256} * \frac{35}{128}$ | $\frac{35}{256} * \frac{5}{32}$ | $\frac{35}{256} * \frac{9}{64}$ | $\frac{35}{256} * \frac{5}{32}$ | $\frac{35}{256} * \frac{35}{128}$ |
| 5 | $\frac{63}{256} * \frac{35}{128}$ | $\frac{63}{256} * \frac{5}{32}$ | $\frac{63}{256} * \frac{9}{64}$ | $\frac{63}{256} * \frac{5}{32}$ | $\frac{63}{256} * \frac{35}{128}$ |

### Normalized form of above table
the below comes from above table $* 2^{15} = 32768$

| | $s_1 = 0$ | 1 | 2 | 3 | 4 |
| --------- | ---------- | - | - | - | - |
| $s_2 = 0$ | $2205$ | $1260$ | $1134$ | $1260$ | $2205$ |
| 1 | $1225$ | $700$ | $630$ | $700$ | $1225$ |
| 2 | $1050$ | $600$ | $540$ | $600$ | $1050$ |
| 3 | $1050$ | $600$ | $540$ | $600$ | $1050$ |
| 4 | $1225$ | $700$ | $630$ | $700$ | $1225$ |
| 5 | $2205$ | $1260$ | $1134$ | $1260$ | $2205$ |

### Simulation Output
:::info

gcc -o Seperate_experiment  Seperate_experiment.c  GSLfun.c -lgsl -lgslcblas -lm
./Seperate_experiment

:::
total sample points: 32768000
Using default random seed
0 2203936 1223708 1050776 1051424 1224938 2206024
1 1257968 701588 599610 601165 699895 1259920
2 1136410 630719 540751 539500 630802 1133816
3 1259445 699621 600675 600773 699740 1259132
4 2204491 1223057 1050313 1050340 1224293 2203170

## Combined trials

### Possible outcomes of 9 trials

| s | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 |
| - | - | - | - | - | - | - | - | - | - | - |
| $C^{s - 0.5}_{s}$ | 1 | $\frac{1}{2}$ | $\frac{3}{8}$ | $\frac{5}{16}$ | $\frac{35}{128}$ | $\frac{63}{256}$ | $\frac{231}{1024}$ | $\frac{429}{2048}$ | $\frac{6435}{32768}$ | $\frac{12155}{65536}$ |
| $P[s]$ | $\frac{12155}{65536}$ | $\frac{6435}{65536}$ | $\frac{19305}{262144}$ | $\frac{2145}{32768}$ | $\frac{2205}{32768}$ | $\frac{2205}{32768}$ | $\frac{2145}{32768}$ | $\frac{19305}{262144}$ | $\frac{6435}{65536}$ | $\frac{12155}{65536}$ |
| $262144*P[s]$ | $48620$ | $25740$ | $19305$ | $17160$ | $17640$ | $17640$ | $17160$ | $19305$ | $25740$ | $48620$ |

### Combination table
Combination table are calculated using python code
:::spoiler Python Code

# combinations.py
from math import factorial

def nCk(n,k): 
  return factorial(n) / (factorial(k) * factorial(n - k))

S1_element = 4
S2_element = 5
max_number = 3

tuple_list = list()
combination_list = list()
for i in range(max(max_number-S2_element, 0),min(S1_element+1, max_number + 1)):
    tuple_list.append((i, max_number - i))
print(tuple_list)

for item in tuple_list:
    combination_list.append((int(nCk(S1_element, item[0])), int(nCk(S2_element, item[1]))))
print(combination_list)

ratio_list = list()
for item in combination_list:
   ratio_list.append(item[0] * item[1])
ratio_sum = sum(ratio_list)
print(ratio_list)
print(ratio_sum)

:::

| s | combinations($s_1, s_2$) | ratio | probability |
| - | ------------ | ------ | ----------- |
| 0 | $(0, 0)$ | $1$ | |
| 1 | $(0, 1), (1, 0)$ | $5 : 4$ | $\frac{5}{9},\frac{4}{9}$ |
| 2 | $(0, 2), (1, 1), (2, 0)$ | $10 : 20 : 6$ | $\frac{5}{18}, \frac{5}{9}, \frac{1}{6}$ |
| 3 | $(0, 3), (1, 2), (2, 1), (3, 0)$ | $10 : 40 : 30: 4$ | $\frac{5}{42}, \frac{10}{21}, \frac{15}{42}, \frac{1}{21}$ |
| 4 | $(0, 4), (1, 3), (2, 2), (3, 1), (4, 0)$ | $5 : 40 : 60 : 20 : 1$ | $\frac{5}{126}, \frac{20}{63}, \frac{10}{21}, \frac{10}{63}, \frac{1}{126}$ |
| 5 | $(0, 5), (1, 4), (2, 3), (3, 2), (4, 1)$ | $1 : 20 : 60 : 40 : 5$ | $\frac{1}{126}, \frac{10}{63}, \frac{10}{21}, \frac{20}{63}, \frac{5}{126}$ |
| 6 | $(1, 5), (2, 4), (3, 3), (4, 2)$ | $4 : 30 : 40 : 10$ | $\frac{1}{21}, \frac{15}{42}, \frac{10}{21}, \frac{5}{42}$ |
| 7 | $(2, 5), (3, 4), (4, 3)$ | $6 : 20 : 10$ | $\frac{1}{6}, \frac{5}{9}, \frac{5}{18}$ |
| 8 | $(3, 5), (4, 4)$ | $4 : 5$ | $\frac{4}{9},\frac{5}{9}$ |
| 9 | $(4, 5)$ | $1$ | |

### Probability Table
The following table was generated using the code below, some part was from modified version of code above.
:::spoiler Code

# probability.py
from math import factorial

def nCk(n,k): 
  return factorial(n) / (factorial(k) * factorial(n - k))

S1_element = 4
S2_element = 5
max_length = S1_element + S2_element

tuple_matrix = list()
ratio_matrix = list()
sum_list = list()
for current_index in range(max_length + 1):
  tuple_list = list()
  combination_list = list()
  for i in range(max(current_index-S2_element, 0),min(S1_element+1, current_index + 1)):
      tuple_list.append((i, current_index - i))
  print(tuple_list)

  for item in tuple_list:
      combination_list.append((int(nCk(S1_element, item[0])), int(nCk(S2_element, item[1]))))
  print(combination_list)

  ratio_list = list()
  for item in combination_list:
    ratio_list.append(item[0] * item[1])
  ratio_sum = sum(ratio_list)
  print(ratio_list)
  print(ratio_sum)

  tuple_matrix.append(tuple_list)
  ratio_matrix.append(ratio_list)
  sum_list.append(ratio_sum)

# I am too lazy to write a code to generate list of probability
import numpy as np
# P_list = [6435, 3432, 2772, 2520, 2450, 2520, 2772, 3432, 6435]
P_list = [48620, 25740, 19305, 17160, 17640, 17640, 17160, 19305, 25740, 48620]

P_matrix = np.zeros((S1_element + 1, S2_element + 1))

for item_list in tuple_matrix:
    for item, index in zip(item_list, range(len(item_list))):
        first_index = item[0]
        second_index = item[1]
        print(first_index, second_index)

        probability = P_list[first_index + second_index]
        ratio = ratio_matrix[first_index + second_index][index]
        denominator = sum_list[first_index + second_index]
        P_matrix[first_index][second_index] = probability * ratio / denominator

# print(P_matrix) 
for item_list in P_matrix:
    print("|", end="")
    for item in item_list:
        print(item, "|", end="")
    print()

:::
| | $s_2 = 0$ | 1 | 2 | 3 | 4 | 5 |
| --------- | ---------- | - | - | - | - | - |
| $s_1 = 0$ |48620.0 |14300.0 |5362.5 |2042.86 |700.0 |140.0 |
| 1 |11440.0 |10725.0 |8171.43 |5600.0 |2800.0 |817.14 |
| 2 |3217.5 |6128.57 |8400.0 |8400.0 |6128.57 |3217.5 |
| 3 |817.14 |2800.0 |5600.0 |8171.43 |10725.0 |11440.0 |
| 4 |140.0 |700.0 |2042.86 |5362.5 |14300.0 |48620.0 |

### Simulation Output
:::warning

gcc -o Combined_experiment  Combined_experiment.c  GSLfun.c -lgsl -lgslcblas -lm
./Combined_experiment

This could take a while......
:::

total sample points: 262144000
Using default random seed
0 48622130 14298062 5716439 2200221 698170 139973
1 11438522 11441474 8802979 5596545 2799594 880475
2 3430081 6598164 8400003 8403540 6604028 3428588
3 879872 2800758 5600501 8800023 11443106 11438750
4 140201 701729 2203215 5719598 14292911 48624348

###### tags: `1112_courses` `probability models`