Added autocorrelation function and covariance code

doc/src/Statistics/Statistics.do.txt

@@ -1758,6 +1758,72 @@ the true $\angle\theta\rangle$. As final result for the observable one quotes $\
 !ec
 
 
+!split
+===== Autocorrelation function =====
+!bc pycod
+# Importing various packages
+from math import exp, sqrt
+from random import random, seed
+import numpy as np
+import matplotlib.pyplot as plt
+
+def autocovariance(x, n, k, mean_x):
+    sum = 0.0
+    for i in range(0, n-k):
+        sum += (x[(i+k)]-mean_x)*(x[i]-mean_x)
+    return  sum/n
+
+n = 1000
+x=np.random.normal(size=n)
+autocor = np.zeros(n)
+figaxis = np.zeros(n)
+mean_x=np.mean(x)
+var_x = np.var(x)
+print(mean_x, var_x)
+for i in range (0, n):
+    figaxis[i] = i
+    autocor[i]=(autocovariance(x, n, i, mean_x))/var_x    
+
+plt.plot(figaxis, autocor, "r-")
+plt.axis([0,n,-0.1, 1.0])
+plt.xlabel(r'$i$')
+plt.ylabel(r'$\gamma_i$')
+plt.title(r'Autocorrelation function')
+plt.show()
+
+
+!ec
+
+!split
+=====  Covariance =====
+!bc pycod
+# Importing various packages
+from math import exp, sqrt
+from random import random, seed
+import numpy as np
+import matplotlib.pyplot as plt
+
+def covariance(x, y, n):
+    sum = 0.0
+    mean_x = np.mean(x)
+    mean_y = np.mean(y)
+    for i in range(0, n):
+        sum += (x[(i)]-mean_x)*(y[i]-mean_y)
+    return  sum/n
+
+n = 10
+
+x=np.random.normal(size=n)
+y = 4+3*x+np.random.normal(size=n)
+covxy = covariance(x,y,n)
+print(covxy)
+z = np.vstack((x, y))
+c = np.cov(z.T)
+
+print(c)
+
+!ec
+