Add sheet8.py

Jaslo/sheet8.py

+from timeit import timeit
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from scipy.special import factorial
+
+from tqdm import tqdm
+
+# Problem 8.1
+# c)
+N_0 = 0
+N_try = int(1e6)
+zeta = 10
+
+N = np.zeros(N_try + 1, dtype=int)
+N[0] = N_0
+for i in range(N_try):
+    if np.random.random() < 0.5:
+        if np.random.random() < zeta / (N[i] + 1):
+            N[i+1] = N[i] + 1
+            continue
+    else:
+        if np.random.random() < N[i] / zeta:
+            N[i+1] = N[i] - 1
+            continue
+
+    N[i+1] = N[i]
+
+print(f"<N> = {np.mean(N)} ± {np.var(N)}")
+
+hist, bins = np.histogram(N, bins=np.max(N))
+
+x = bins[:-1]
+plt.plot(x, hist / np.max(hist), label="MC", zorder=3)
+
+x = np.linspace(x[0], x[-1])
+y = zeta**x * np.exp(-zeta) / factorial(x)
+plt.plot(x, y / np.max(y), label="Poisson")
+plt.xlabel(r"$N$")
+plt.ylabel(r"$p^\textrm{eq}(N)$ (normalised)")
+plt.legend()
+plt.savefig("problem8.1c.png")
+plt.show()
+
+# d)
+N = np.zeros(N_try + 1, dtype=int)
+N[0] = N_0
+for i in range(N_try):
+    if np.random.random() < 0.5:
+        if np.random.random() < zeta**3 / ((N[i] + 1) * (N[i] + 2) * (N[i] + 3)):
+            N[i+1] = N[i] + 3
+            continue
+    else:
+        if np.random.random() < N[i] * (N[i] + 1) * (N[i] + 2) / zeta**3 and N[i] > 2:
+            N[i+1] = N[i] - 3
+            continue
+    N[i+1] = N[i]
+
+print(f"<N> = {np.mean(N)} ± {np.var(N)}")
+
+hist, bins = np.histogram(N, bins=np.max(N))
+
+x = bins[:-1]
+plt.plot(x, hist / np.max(hist), label="MC")
+
+x = np.linspace(x[0], x[-1])
+y = zeta**x * np.exp(-zeta) / factorial(x)
+plt.plot(x, y / np.max(y), label="Poisson")
+plt.xlabel(r"$N$")
+plt.ylabel(r"$p^\textrm{eq}(N)$ (normalised)")
+plt.legend()
+plt.savefig("problem8.1d-3-10.png")
+plt.show()
+
+
+zeta = 100
+N = np.zeros(N_try + 1, dtype=int)
+N[0] = N_0
+for i in range(N_try):
+    if np.random.random() < 0.5:
+        if np.random.random() < zeta**3 / ((N[i] + 1) * (N[i] + 2) * (N[i] + 3)):
+            N[i+1] = N[i] + 3
+            continue
+    else:
+        if np.random.random() < N[i] * (N[i] + 1) * (N[i] + 2) / zeta**3 and N[i] > 2:
+            N[i+1] = N[i] - 3
+            continue
+    N[i+1] = N[i]
+
+print(f"<N> = {np.mean(N)} ± {np.var(N)}")
+
+hist, bins = np.histogram(N, bins=np.max(N))
+
+x = bins[:-1]
+plt.plot(x, hist / np.max(hist), label="MC")
+
+x = np.linspace(x[0], x[-1])
+y = zeta**x * np.exp(-zeta) / factorial(x)
+plt.plot(x, y / np.max(y), label="Poisson")
+plt.xlabel(r"$N$")
+plt.ylabel(r"$p^\textrm{eq}(N)$ (normalised)")
+plt.legend()
+plt.savefig("problem8.1d-3-100.png")
+plt.show()
+
+
+zeta = 10
+M = 3
+N = np.zeros(N_try + 1, dtype=int)
+N[0] = N_0
+for i in range(N_try):
+    m = np.random.randint(1, M+1)
+    if np.random.random() < 0.5:
+        if np.random.random() < zeta**m / np.multiply.accumulate(np.arange(N[i] + 1, N[i] + m + 1))[-1]:
+            N[i+1] = N[i] + m
+            continue
+    else:
+        if np.random.random() < np.multiply.accumulate(np.arange(N[i], N[i] + m))[-1] / zeta**m and N[i] >= m:
+            N[i+1] = N[i] - m
+            continue
+    N[i+1] = N[i]
+
+print(f"<N> = {np.mean(N)} ± {np.var(N)}")
+
+hist, bins = np.histogram(N, bins=np.max(N))
+
+x = bins[:-1]
+plt.plot(x, hist / np.max(hist), label="MC")
+
+x = np.linspace(x[0], x[-1])
+y = zeta**x * np.exp(-zeta) / factorial(x)
+plt.plot(x, y / np.max(y), label="Poisson")
+plt.xlabel(r"$N$")
+plt.ylabel(r"$p^\textrm{eq}(N)$ (normalised)")
+plt.legend()
+plt.savefig("problem8.1d-m-100.png")
+plt.show()
+
+# e)
+zeta = 10
+N = np.zeros(N_try + 1, dtype=int)
+N[0] = N_0
+for i in range(N_try):
+    if np.random.random() < 0.5:
+        if np.random.random() < zeta / (N[i] + 1):
+            N[i+1] = N[i] + 1
+            continue
+    else:
+        if np.random.random() < N[i] / zeta:
+            N[i+1] = N[i] - 1
+            continue
+
+    N[i+1] = N[i]
+
+k = np.arange(100, N_try - 20, 20, dtype=int)
+C = np.zeros((len(k), 21))
+l = np.arange(21, dtype=int)
+for i in tqdm(l):
+    for j in range(len(k)):
+        C[j,i] = np.mean(N[k[j] + i] * N[k[j]]) - np.mean(N)**2
+
+plt.plot(l, np.mean(C, axis=0))
+plt.xlabel(r"$l$")
+plt.ylabel(r"$C(l)$")
+plt.savefig("problem8.1e.png")
+plt.show()
+
+# f)
+def mc():
+    zeta = 1000
+    M = 3
+    N = np.zeros(N_try + 1, dtype=int)
+    N[0] = N_0
+    for i in range(N_try):
+        m = np.random.randint(1, M+1)
+        if np.random.random() < 0.5:
+            if np.random.random() < zeta**m / np.multiply.accumulate(np.arange(N[i] + 1, N[i] + m + 1))[-1]:
+                N[i+1] = N[i] + m
+                continue
+        else:
+            if np.random.random() < np.multiply.accumulate(np.arange(N[i], N[i] + m))[-1] / zeta**m and N[i] >= m:
+                N[i+1] = N[i] - m
+                continue
+        N[i+1] = N[i]
+
+def poisson():
+    zeta = 1000
+    l = np.exp(-zeta)
+    for _ in range(N_try):
+        p = 1
+        k = 0
+
+        while True:
+            p *= np.random.random()
+            k += 1
+
+            if p <= l:
+                break