diff --git a/07_project_stats.py b/07_project_stats.py
new file mode 100644
index 0000000..9b78b25
--- /dev/null
+++ b/07_project_stats.py
@@ -0,0 +1,77 @@
+import numpy as np
+from qiskit import (
+    QuantumCircuit,
+    execute,
+    Aer)
+from qiskit.visualization import plot_histogram, plot_bloch_vector
+from matplotlib import pyplot as plt
+from collections import defaultdict
+
+# Use Aer's qasm_simulator
+simulator = Aer.get_backend('qasm_simulator')
+
+
+def create_bloch_vector(theta, phi):
+    x = np.sin(theta) * np.cos(phi)
+    y = np.sin(theta) * np.sin(phi)
+    z = np.cos(theta)
+    return [x, y, z]
+
+
+def perform_exp(iterations, case, draw=False):
+    simulator = Aer.get_backend('qasm_simulator')
+    all_counts = defaultdict(int)
+
+    for i in range(iterations):
+        qc = QuantumCircuit(2, 2)
+
+        # uniform distribution on the sphere
+        t = round(np.random.uniform(0, 1), 10)
+        theta = np.arccos(1 - 2 * t)
+        phi = round(np.random.uniform(0, 2 * np.pi), 10)
+
+        CASE_ANGLES = {
+            'identical': (theta, phi),
+            'orthogonal': (theta + np.pi, phi)
+        }
+
+        qc.r(theta, phi, 0)
+        theta1, phi1 = CASE_ANGLES[case]
+        qc.r(theta1, phi1, 1)
+
+        # print(theta, phi, phi1)
+
+        if draw:
+            q0 = create_bloch_vector(theta, phi)
+            plot_bloch_vector(q0)
+            plt.show()
+            q1 = create_bloch_vector(theta1, phi1)
+            plot_bloch_vector(q1)
+            plt.show()
+
+        # Measure in the Bell's basis
+        qc.cx(0, 1)
+        qc.h(0)
+        qc.measure([0, 1], [0, 1])
+
+        # print(qc)
+
+        # execute the job
+        job = execute(qc, simulator, shots=1)
+        result = job.result()
+        counts = result.get_counts(qc)
+
+        for k, v in counts.items():
+            all_counts[k] += v
+
+    # gather all counts
+    print("\nTotal counts for experiment '{}' are:\n".format(case))
+    for k, v in sorted(all_counts.items(), key=lambda x: x[0]):
+        print("{}: {}".format(k, v))
+
+
+if __name__ == "__main__":
+    # draw produces a Bloch sphere but does not work in this notebook
+
+    perform_exp(iterations=100, case='identical', draw=False)  # same state
+    perform_exp(iterations=100, case='orthogonal', draw=False)  # different orthogonal states