first commit

.idea/.gitignore

@@ -0,0 +1,8 @@
+# Default ignored files
+/shelf/
+/workspace.xml
+# Editor-based HTTP Client requests
+/httpRequests/
+# Datasource local storage ignored files
+/dataSources/
+/dataSources.local.xml

.idea/inspectionProfiles/profiles_settings.xml

@@ -0,0 +1,6 @@
+<component name="InspectionProjectProfileManager">
+  <settings>
+    <option name="USE_PROJECT_PROFILE" value="false" />
+    <version value="1.0" />
+  </settings>
+</component>

.idea/misc.xml

@@ -0,0 +1,7 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<project version="4">
+  <component name="Black">
+    <option name="sdkName" value="Python 3.10 (STA3)" />
+  </component>
+  <component name="ProjectRootManager" version="2" project-jdk-name="Python 3.10 (STA3)" project-jdk-type="Python SDK" />
+</project>

.idea/modules.xml

@@ -0,0 +1,8 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<project version="4">
+  <component name="ProjectModuleManager">
+    <modules>
+      <module fileurl="file://$PROJECT_DIR$/.idea/pythonProject1.iml" filepath="$PROJECT_DIR$/.idea/pythonProject1.iml" />
+    </modules>
+  </component>
+</project>

.idea/pythonProject1.iml

@@ -0,0 +1,10 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<module type="PYTHON_MODULE" version="4">
+  <component name="NewModuleRootManager">
+    <content url="file://$MODULE_DIR$">
+      <excludeFolder url="file://$MODULE_DIR$/venv" />
+    </content>
+    <orderEntry type="jdk" jdkName="Python 3.10 (STA3)" jdkType="Python SDK" />
+    <orderEntry type="sourceFolder" forTests="false" />
+  </component>
+</module>

.idea/ruff.xml

@@ -0,0 +1,6 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<project version="4">
+  <component name="RuffConfigService">
+    <option name="globalRuffExecutablePath" value="$USER_HOME$/.local/bin/ruff" />
+  </component>
+</project>

.idea/vcs.xml

@@ -0,0 +1,6 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<project version="4">
+  <component name="VcsDirectoryMappings">
+    <mapping directory="$PROJECT_DIR$" vcs="Git" />
+  </component>
+</project>

EVM.py

@@ -0,0 +1,28 @@
+import random
+
+
+N = 100000
+TA = 3000
+TB = 6000
+T1 = 1200
+T2 = 1500
+k = 4548
+
+A = random.randint(0, N)
+B = N - A
+Na = k*10-B/2
+Nb = (k-0.1*A)/0.05
+
+print(f"A: {A}, B: {B}")
+print(f"A: {Na}, B: {Nb}")
+
+# pravděpodobnost rozpadu jednotlivých prvků v intervalu T1,T2
+RA = (T2-T1)/TA  # 0.1
+RB = (T2-T1)/TB  # 0.05
+
+# počet rozpadlých atomů v intervalu
+NB = 2*k/3  # 3032
+NA = k-NB  # 1516
+
+print(RA, RB)
+print(NB, k-NB)

FFT.py

@@ -0,0 +1,47 @@
+# DFT vs FFT - porovnani casu
+from cmath import exp, pi
+import numpy as np
+import time
+
+
+def fft(f):
+    N = len(f)
+
+    if N <= 1:
+        return f
+
+    even = fft(f[0::2])
+    odd = fft(f[1::2])
+
+    T = [exp(-2j * pi * k / N) * odd[k] for k in range(N // 2)]
+
+    return [even[k] + T[k] for k in range(N // 2)] + [even[k] - T[k] for k in range(N // 2)]
+
+
+def dft(f):
+    N = len(f)
+    result = [0] * N
+
+    for i in range(N):
+        for j in range(N):
+            result[i] += f[j] * exp(-2j * pi * i * j / N)
+    return result
+
+
+print("i,  2**i      DFT,      FFT,    numpy")
+for i in range(10, 21):
+    f = [1] * 2**i
+    start_time_D = end_time_D = 0
+    if i < 15:
+        start_time_D = time.time()
+        dft(f)
+        end_time_D = time.time()
+
+    start_time_F = time.time()
+    fft(f)
+    end_time_F = time.time()
+
+    start_time_N = time.time()
+    np.fft.fft(f)
+    end_time_N = time.time()
+    print(i, 2**i, end_time_D - start_time_D, end_time_F - start_time_F, end_time_N - start_time_N)

Fourier.py

@@ -0,0 +1,35 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+# Délka jedné periody
+T = 2
+
+# Počet period před a po intervalu
+num_periods = 3
+
+# Definice intervalu t pro funkci <0, 1/2)
+t1 = np.linspace(0, 0.5, 100)
+
+# Definice intervalu t pro funkci <1/2, 2)
+t2 = np.linspace(0.5, 2, 100)
+
+# Vytvořte pole hodnot f(t) pro první interval
+f1 = 1 + np.exp(-t1)
+
+# Vytvořte pole hodnot f(t) pro druhý interval a posuňte ho tak, aby měl hodnotu 2 v bodech vzdálených o násobek periody T = 2
+f2 = 1 + np.exp(-t2 - T * num_periods)
+
+# Vykreslete první část funkce
+plt.plot(t1, f1, label='1+exp(-t)', color='blue')
+
+# Vykreslete druhou část funkce
+plt.plot(t2, f2, label='1+exp(-t) shifted', color='blue')
+
+# Nastavte osy a přidejte popisky
+plt.xlabel('t')
+plt.ylabel('f(t)')
+plt.title('Graf funkce f(t) s posunutím')
+plt.legend()
+
+# Zobrazte graf
+plt.show()

ITDT/Gabor.py

@@ -0,0 +1,99 @@
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy import ndimage as ndi
+
+from skimage.io import imread
+from skimage.util import img_as_float
+from skimage.filters import gabor_kernel
+
+
+def compute_feats(image, kernels):
+    feats = np.zeros((len(kernels), 2), dtype=np.double)
+    for k, kernel in enumerate(kernels):
+        filtered = ndi.convolve(image, kernel, mode='wrap')
+        feats[k, 0] = filtered.mean()
+        feats[k, 1] = filtered.var()
+    return feats
+
+
+def match(feats, ref_feats):
+    min_error = np.inf
+    min_i = None
+    for i in range(ref_feats.shape[0]):
+        error = np.sum((feats - ref_feats[i, :])**2)
+        if error < min_error:
+            min_error = error
+            min_i = i
+    return min_i
+
+
+sigma = 10  # 1 - 10
+frequency = 0.35  # 0.05 - 0.35
+
+
+kernels = []
+for theta in range(4):
+    theta = theta / 4. * np.pi
+    kernel = np.real(gabor_kernel(frequency, theta=theta, n_stds=sigma))
+    kernels.append(kernel)
+
+
+shrink = (slice(0, None, 3), slice(0, None, 3))
+disc = img_as_float(imread("discord-logo.png", as_gray=True))[shrink]
+feng = img_as_float(imread("Feng_wind.jpg", as_gray=True))[shrink]
+vazka = img_as_float(imread("vazka.png", as_gray=True))[shrink]
+image_names = ('Discord', 'Znak', 'Vážka')
+images = (disc, feng, vazka)
+
+# prepare reference features
+ref_feats = np.zeros((3, len(kernels), 2), dtype=np.double)
+ref_feats[0, :, :] = compute_feats(disc, kernels)
+ref_feats[1, :, :] = compute_feats(feng, kernels)
+ref_feats[2, :, :] = compute_feats(vazka, kernels)
+
+def power(image, kernel):
+    # Normalize images for better comparison.
+    image = (image - image.mean()) / image.std()
+    return np.sqrt(ndi.convolve(image, np.real(kernel), mode='wrap')**2 +
+                   ndi.convolve(image, np.imag(kernel), mode='wrap')**2)
+
+# Plot a selection of the filter bank kernels and their responses.
+results = []
+kernel_params = []
+for theta in range(0, 10):
+    theta = theta / np.pi
+    kernel = gabor_kernel(frequency, theta=theta)
+    params = f"Θ={18*(theta*np.pi):.0f}°"
+    kernel_params.append(params)
+    # Save kernel and the power image for each image
+    results.append((kernel, [power(img, kernel) for img in images]))
+
+fig, axes = plt.subplots(nrows=len(results)+1, ncols=len(image_names)+1, figsize=(4, 9))
+plt.gray()
+
+fig.suptitle(f'Gaborovy filtry\n pro λ={frequency}, ϕ={0}, σ={sigma} a γ={1}.', fontsize=12)
+
+axes[0][0].axis('off')
+
+# Plot original images
+for label, img, ax in zip(image_names, images, axes[0][1:]):
+    ax.imshow(img)
+    ax.set_title(label, fontsize=9)
+    ax.axis('off')
+
+for label, (kernel, powers), ax_row in zip(kernel_params, results, axes[1:]):
+    # Plot Gabor kernel
+    ax = ax_row[0]
+    ax.imshow(np.real(kernel))
+    ax.set_ylabel(label, fontsize=7)
+    ax.set_xticks([])
+    ax.set_yticks([])
+
+    # Plot Gabor responses with the contrast normalized for each filter
+    vmin = np.min(powers)
+    vmax = np.max(powers)
+    for patch, ax in zip(powers, ax_row[1:]):
+        ax.imshow(patch, vmin=vmin, vmax=vmax)
+        ax.axis('off')
+
+plt.show()

ITDT/fourier/fourier.py

@@ -0,0 +1,62 @@
+import numpy as np
+import sympy as sp
+import matplotlib.pyplot as plt
+
+t = sp.symbols('t')
+T = 2  # Perioda funkce
+to = 1/2
+omega = 2 * sp.pi / T
+
+f = 1 + sp.exp(-t)
+
+def fc(x):
+    return 1 + sp.exp(-x)
+
+
+def get_tabulka():
+    print(f"a_0 = {a_0: .2f}")
+    print(f"c_0 = {((1 / T) * sp.integrate(f, (t, 0, to))):.2f}")
+    print(f"A_0 = {(a_0 / 2):.2f}")
+
+    pa_coeffs = [((2 / T) * sp.integrate(f * sp.cos(n * omega * t), (t, 0, to)).evalf(), n) for n in range(1, n_max + 1)]
+    pb_coeffs = [((2 / T) * sp.integrate(f * sp.sin(n * omega * t), (t, 0, to)).evalf(), n) for n in range(1, n_max + 1)]
+    pc_coeffs = [((1 / T) * sp.integrate(f * sp.exp(1j * n * omega * t), (t, 0, to)).evalf(), n) for n in range(1, n_max + 1)]
+    for a, b, c in zip(pa_coeffs, pb_coeffs, pc_coeffs):
+        print(f"n= {a[1]}")
+        print(f"a= {a[0]:.2f}")
+        print(f"b= {b[0]:.2f}")
+        print(f"A= {(a[0]**2+b[0]**2)**(1/2):.2f}")
+        print(f"c_n= {np.absolute(c[0]):.2f}")
+        print(f"phi= {np.angle(float(a[0])+float(b[0])*1j, deg=True):.2f}")
+
+def calculate_coefficients(n_max):
+    a_0 = (2 / T) * sp.integrate(f, (t, 0, to))
+    a_coeffs = [((2 / T) * sp.integrate(f * sp.cos(n * omega * t), (t, 0, to)), n) for n in range(1, n_max + 1)]
+    b_coeffs = [((2 / T) * sp.integrate(f * sp.sin(n * omega * t), (t, 0, to)), n) for n in range(1, n_max + 1)]
+    return a_0, a_coeffs, b_coeffs
+
+def fourier_series(x, a_0, a_coeffs, b_coeffs):
+    series = a_0 / 2
+    for a_n, n in a_coeffs:
+        series += a_n * sp.cos(n * omega * x)
+    for b_n, n in b_coeffs:
+        series += b_n * sp.sin(n * omega * x)
+    return series
+
+
+for n_max in [5, 10, 20, 30]:
+    a_0, a_coeffs, b_coeffs = calculate_coefficients(n_max)
+
+    # Definujte interval pro vykreslení
+    t_values = np.linspace(-4, 6, 200)
+
+    # Výpočet Fourierovy řady a vykreslení
+    fourier = [fourier_series(x, a_0, a_coeffs, b_coeffs) for x in t_values]
+    original = [fc(x%T) if (x%T) < to else 0 for x in t_values]
+    plt.plot(t_values, original, label="f(t)")
+    plt.plot(t_values, fourier, label="FŘ")
+    plt.title(f"n={n_max}")
+    plt.legend()
+    plt.xlabel("t")
+    plt.ylabel("f(t)")
+    plt.show()

ITDT/fourier/sinus_cosinus.py

@@ -0,0 +1,78 @@
+import numpy as np
+import sympy as sp
+import matplotlib.pyplot as plt
+
+t = sp.symbols('t')
+T = 4  # Perioda funkce
+to = 1/2
+omega = 2 * sp.pi / T
+
+f = 1 + sp.exp(-t)
+
+def fc(x):
+    return 1 + sp.exp(-x)
+
+def get_tabulka():
+    print(f"a_0 = {a_0: .2f}")
+    print(f"c_0 = {((1 / T) * sp.integrate(f, (t, 0, to))):.2f}")
+    print(f"A_0 = {(a_0 / 2):.2f}")
+
+    pa_coeffs = [((2 / T) * sp.integrate(f * sp.cos(n * omega * t), (t, 0, to)).evalf(), n) for n in range(1, n_max + 1)]
+    pb_coeffs = [((2 / T) * sp.integrate(f * sp.sin(n * omega * t), (t, 0, to)).evalf(), n) for n in range(1, n_max + 1)]
+    pc_coeffs = [((1 / T) * sp.integrate(f * sp.exp(1j * n * omega * t), (t, 0, to)).evalf(), n) for n in range(1, n_max + 1)]
+    for a, b, c in zip(pa_coeffs, pb_coeffs, pc_coeffs):
+        print(f"n= {a[1]}")
+        print(f"a= {a[0]:.2f}")
+        print(f"b= {b[0]:.2f}")
+        print(f"A= {(a[0]**2+b[0]**2)**(1/2):.2f}")
+        print(f"c_n= {np.absolute(c[0]):.2f}")
+        print(f"phi= {np.angle(float(a[0])+float(b[0])*1j, deg=True):.2f}")
+
+def calculate_coefficients(n_max):
+    a_0 = (2 / T) * sp.integrate(f, (t, 0, to))
+    a_coeffs = [((2 / T) * sp.integrate(f * sp.cos(n * omega * t), (t, 0, to)), n) for n in range(1, n_max + 1)]
+    b_coeffs = [((2 / T) * sp.integrate(f * sp.sin(n * omega * t), (t, 0, to)), n) for n in range(1, n_max + 1)]
+    return a_0, a_coeffs, b_coeffs
+
+def cosinus_series(x, a_0, a_coeffs):
+    cosinus = a_0 / 2
+    for a_n, n in a_coeffs:
+        cosinus += a_n * sp.cos(omega * n * x)
+    return cosinus
+
+def sinus_series(x, b_coeffs):
+    sinus = 0
+    for b_n, n in b_coeffs:
+        sinus += b_n * sp.sin(omega * n * x)
+    return sinus
+
+# Předpočítat koeficienty
+n_maxs = [5, 10, 20, 30]
+
+for n_max in n_maxs:
+    a_0, a_coeffs, b_coeffs = calculate_coefficients(n_max)
+
+    # Definujte interval pro vykreslení
+    t_values = np.linspace(-4, 6, 200)
+
+    cosinus = [cosinus_series(x, a_0, a_coeffs) for x in t_values]
+    cosinus_scaled = [2 * y for y in cosinus]
+    original = [fc(-x % (2*2)) if (-x % (2*2)) < to else fc(x % (2*2)) if (x % (2*2)) < to else 0 for x in t_values]
+    plt.plot(t_values, original, label="f(t)")
+    plt.plot(t_values, cosinus_scaled, label="CŘ")
+    plt.title(f"n={n_max}")
+    plt.legend()
+    plt.xlabel("t")
+    plt.ylabel("f(t)")
+    plt.show()
+
+    sinus = [sinus_series(x, b_coeffs) for x in t_values]
+    sinus_scaled = [2 * y for y in sinus]
+    original = [-fc(-x % (2*2)) if (-x % (2*2)) < to else fc(x % (2*2)) if (x % (2*2)) < to else 0 for x in t_values]
+    plt.plot(t_values, original, label="f(t)")
+    plt.plot(t_values, sinus_scaled, label="SŘ")
+    plt.title(f"n={n_max}")
+    plt.legend()
+    plt.xlabel("t")
+    plt.ylabel("f(t)")
+    plt.show()

ITDT/jjj.py

@@ -0,0 +1,31 @@
+import sympy as sp
+
+# Symbolická proměnná t
+t = sp.symbols('t')
+
+
+# Definujte funkci f(t) podle vašich potřeb
+def f(t):
+    return t if 0 < t < 1 else 0
+
+# Zde nastavte hodnotu n pro konkrétní číslo harmonické složky
+n = 1
+
+# Perioda a úhlová frekvence
+T = 2  # Perioda funkce
+omega = 2 * sp.pi / T  # Úhlová frekvence
+
+# Výpočet Fourierových koeficientů
+a_0 = (2 / T) * sp.integrate(f(t), (t, 0, T))
+a_n = (2 / T) * sp.integrate(f(t) * sp.cos(n * omega * t), (t, 0, T))
+b_n = (2 / T) * sp.integrate(f(t) * sp.sin(n * omega * t), (t, 0, T))
+
+
+# Výpočet konkrétního koeficientu c_n na základě a_n a b_n
+c_n = a_n + b_n
+
+# Tisk koeficientů
+print(f'a_0: {a_0}')
+print(f'a_{n}: {a_n.subs(n, 1)}')
+print(f'b_{n}: {b_n.subs(n, 1)}')
+print(f'c_{n}: {c_n.subs(n, 1)}')

convolution1D.py

@@ -0,0 +1,39 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+#  vytvoreni signalu
+t = np.linspace(0, 2, 200)
+clean_signal = np.sin(2 * np.pi * 1 * t)
+
+#  zasumeni signalu
+noise = 0.1 * np.random.randn(len(t))
+noisy_signal = clean_signal + noise
+
+
+#  1D konvoluce
+def convolution(signal, clean):
+  c = []
+  for i, bv in enumerate(clean):
+    for j, value in enumerate(signal):
+      if i+j < len(c):
+        c[i+j] += value*bv
+      else:
+        c.append(value*bv)
+  return c
+
+
+#  vytvoření filtru na signál
+def filter_signal(signal, number):
+  return convolution(signal, [1/number for _ in range(number)])
+
+
+#  první obrázek - vytvořený signál
+plt.plot(noisy_signal)
+plt.show()
+#  druhý obrázek - signál po filtru n=5
+plt.plot(filter_signal(noisy_signal, 5))
+plt.show()
+#  třetí obrázek - signál po filtru n=20
+plt.plot(filter_signal(noisy_signal, 20))
+plt.show()
+

convolution2D.py

@@ -0,0 +1,40 @@
+from PIL import Image
+
+
+def convolution(input_image, image_filter):
+    width, height = input_image.size
+    output_image = Image.new("RGB", (width, height))
+
+    for x in range(width):
+        for y in range(height):
+            red, green, blue = 0, 0, 0
+            div = 0
+            filter_from = -int(len(image_filter[0])/2)
+            filter_to = len(image_filter[0]) + filter_from
+            for i in range(filter_from, filter_to):
+                for j in range(filter_from, filter_to):
+                    if x + i > 0 and x + i < width-1 and y + j > 0 and y + j < height-1:
+                        neighbor_pixel = input_image.getpixel((x + i, y + j))
+                        red += image_filter[i + 1][j + 1] * neighbor_pixel[0]
+                        green += image_filter[i + 1][j + 1] * neighbor_pixel[1]
+                        blue += image_filter[i + 1][j + 1] * neighbor_pixel[2]
+                        div += image_filter[i + 1][j + 1]
+            if div > 0:
+                red = int(red/div)
+                blue = int(blue/div)
+                green = int(green/div)
+            output_image.putpixel((x, y), (red, green, blue))
+
+    output_image.save("moucha3.jpg")
+    output_image.close()
+    input_image.close()
+
+# moucha3.jpg
+odsumeni = [[1, 2, 1],  [2, 5, 2],  [1, 2, 1]]
+# moucha1.jpg
+sum = [[2, 4, 5, 4, 2], [4, 9, 12, 9, 4], [5, 12, 15, 12, 5], [4, 9, 12, 9, 4], [2, 4, 6, 4, 2]]
+# moucha2.jpg
+hrany = [[0, 1, 0],  [1, -4, 1],  [0, 1, 0]]
+# convolution(Image.open("pop.jpg"), sum)
+convolution(Image.open("moucha.jpg"), odsumeni)
+

fractal.py

@@ -0,0 +1,95 @@
+from cmath import sin
+
+import matplotlib.pyplot as plt
+
+
+def plot_logistic_map(r, x0, num_iterations):
+    x = [x0]
+    for _ in range(num_iterations):
+        x0 = logistic_map(r, x0)
+        x.append(x0)
+
+    plt.plot(range(num_iterations + 1), x, marker='o', linestyle='-')
+    plt.xlabel("Národy")
+    plt.ylabel("Hodnota x")
+    plt.title(f"Logistická mapa (r = {r})")
+    plt.show()
+
+
+def logistic_map(r, x):
+    # return r * x * (1 - x) * (- x + 1)
+    return r * x * (1 - x)
+
+
+def bifurcation_diagram(r_values, x0, num_iterations, num_transient, num_points):
+    all_x = []
+    all_r = []
+
+    for r in r_values:
+        x = x0
+        transient = []
+        stable = []
+
+        for i in range(num_iterations):
+            x = logistic_map(r, x)
+
+            if i >= num_transient:
+                transient.append(r)
+                stable.append(x)
+
+        all_x.extend(stable)
+        all_r.extend(transient)
+
+    plt.scatter(all_r, all_x, s=1, marker=".")
+    plt.xlabel("r")
+    plt.ylabel("Hodnota x")
+    plt.title("Bifurkační diagram")
+    plt.show()
+
+
+def create_diagrams():
+    # Parametry
+    r = 6
+    r_values = [i / 1000 for i in range(1000, 6000)]  # Rozsah hodnot r
+    x0 = 0.2  # Počáteční hodnota
+    num_iterations = 1000
+    num_transient = 100  # Počet přechodných iterací
+    num_points = 10000  # Počet bodů v diagramu
+
+    bifurcation_diagram(r_values, x0, num_iterations, num_transient, num_points)
+    plot_logistic_map(r, x0, 100)
+
+
+def midpoint(point1, point2):
+    return ((point1[0] + point2[0]) / 2, (point1[1] + point2[1]) / 2)
+
+
+def sierpinski_triangle(ax, depth, p1, p2, p3):
+    if depth == 0:
+        ax.plot(*p1, 'ro')
+        ax.plot(*p2, 'ro')
+        ax.plot(*p3, 'ro')
+    else:
+        p12 = midpoint(p1, p2)
+        p23 = midpoint(p2, p3)
+        p31 = midpoint(p3, p1)
+
+        sierpinski_triangle(ax, depth - 1, p1, p12, p31)
+        sierpinski_triangle(ax, depth - 1, p12, p2, p23)
+        sierpinski_triangle(ax, depth - 1, p31, p23, p3)
+
+
+def main():
+    fig, ax = plt.subplots()
+    p1 = (0, 0)
+    p2 = (1, 0)
+    p3 = (0.5, (3 ** 0.5) / 2)
+    depth = 5  # Změňte tuto hodnotu pro různé úrovně fraktálu
+
+    sierpinski_triangle(ax, depth, p1, p2, p3)
+    ax.set_aspect('equal', adjustable='box')
+    plt.show()
+
+
+if __name__ == "__main__":
+    main()

main.py

@@ -0,0 +1,34 @@
+import numpy as np
+import pymc as pm
+import arviz as az
+from matplotlib import pyplot as plt
+
+
+x = np.random.normal(3, 2, size=30)
+# plt.hist(x)
+# plt.show()
+
+with pm.Model() as model:
+    mu = pm.Uniform('mu', lower=-5, upper=5)
+    sig = pm.HalfNormal('sig', sigma=10)
+    data = pm.Normal('data', mu=mu, sigma=sig, observed=x)
+
+pm.model_to_graphviz(model)
+
+with model:
+    trace = pm.sample(5000, tune=1000)
+
+# plt.plot(trace.posterior.data_vars["mu"][0, 1:100])
+
+# az.plot_trace(trace)
+summary = az.summary(trace)
+
+with model:
+    post = pm.sample_posterior_predictive(trace)
+
+plt.hist(post.posterior_predictive.data.to_numpy().flatten(), density=True, bins=50)
+plt.hist(x, alpha=0.4, density=True, bins=6)
+
+az.plot_ppc(post)
+
+plt.show()