{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "source": [
        "# Κεφάλαιο 10 — Στοχαστικές διαδικασίες συνεχούς χρόνου (Ρυθμίσεις)\n",
        "\n",
        "Σε αυτό το notebook μελετάμε στοχαστικές διαδικασίες συνεχούς χρόνου, και\n",
        "πιο συγκεκριμένα Μαρκοβιανές αλυσίδες συνεχούς χρόνου με γεννήτορα\n",
        "πίνακα $Q = (q_{ij})$.\n",
        "\n",
        "Θα χρησιμοποιήσουμε:\n",
        "\n",
        "- την προσομοίωση τροχιών $X(t)$,\n",
        "- την ανάπτυξη σε σειρά\n",
        "  $$P(t) = \\sum_{n=0}^{\\infty} \\frac{t^n}{n!} Q^n,$$\n",
        "- την κατασκευή της ενσωματωμένης αλυσίδας μέσω των πινάκων\n",
        "  $$\\Lambda = (\\Lambda_i), \\quad \\Lambda_i = -q_{ii},$$\n",
        "  και του στοχαστικού πίνακα $R$,\n",
        "- τις εξισώσεις Chapman–Kolmogorov\n",
        "  $$P(t+s) = P(t)\\,P(s),$$\n",
        "- και τη στάσιμη κατανομή $\\pi$ με\n",
        "  $$\\pi Q = 0, \\quad \\sum_i \\pi_i = 1.$$\n",
        "\n",
        "Στο κελί αυτό κάνουμε μόνο τις βασικές ρυθμίσεις (βιβλιοθήκες, seed,\n",
        "μέγεθος διαγραμμάτων).\n"
      ],
      "metadata": {
        "id": "rv13RY80J5QG"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "from numpy.random import default_rng\n",
        "\n",
        "# τυχαίος γεννήτορας για αναπαραγωγιμότητα\n",
        "rng = default_rng(2025)\n",
        "\n",
        "# ρυθμίσεις matplotlib\n",
        "plt.rcParams[\"figure.figsize\"] = (8, 4)\n",
        "plt.rcParams[\"axes.grid\"] = True\n",
        "\n",
        "# Θα χρησιμοποιούμε αυτό τον Q σε όλο το κεφάλαιο:\n",
        "# 3-καταστατική CTMC με καταστάσεις 0,1,2\n",
        "Q = np.array([\n",
        "    [-1.5, 1.0 , 0.5 ],\n",
        "    [ 0.5, -1.0, 0.5 ],\n",
        "    [ 0.25, 0.25, -0.5]\n",
        "])\n",
        "\n",
        "print(\"Γεννήτορας Q:\")\n",
        "print(Q)\n",
        "print(\"Άθροισμα κάθε γραμμής (πρέπει να είναι 0):\", Q.sum(axis=1))\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "ovwzfhkOL01F",
        "outputId": "70736920-7cd7-4bfe-e3f2-5a5924c7e6ca"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Γεννήτορας Q:\n",
            "[[-1.5   1.    0.5 ]\n",
            " [ 0.5  -1.    0.5 ]\n",
            " [ 0.25  0.25 -0.5 ]]\n",
            "Άθροισμα κάθε γραμμής (πρέπει να είναι 0): [0. 0. 0.]\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# 10.1 Εισαγωγή στις στοχαστικές διαδικασίες συνεχούς χρόνου\n",
        "\n",
        "Θυμόμαστε ότι μια στοχαστική διαδικασία συνεχούς χρόνου είναι μια οικογένεια\n",
        "τυχαίων μεταβλητών\n",
        "$$\\{X(t) : t \\ge 0\\},$$\n",
        "ορισμένη σε μια κοινή πιθανωτική δομή, με τιμές σε ένα σύνολο καταστάσεων\n",
        "$S$.\n",
        "\n",
        "Μια Μαρκοβιανή αλυσίδα συνεχούς χρόνου με γεννήτορα $Q = (q_{ij})$ έχει την\n",
        "ιδιότητα:\n",
        "$$\n",
        "\\mathbb{P}(X(t_{n+1}) = i_{n+1} \\mid X(t_n) = i_n, \\dots, X(t_0) = i_0)\n",
        "= \\mathbb{P}(X(t_{n+1}) = i_{n+1} \\mid X(t_n) = i_n),\n",
        "$$\n",
        "και επιπλέον\n",
        "- όταν το σύστημα βρίσκεται στην κατάσταση $i$, ο χρόνος παραμονής εκεί\n",
        "  είναι εκθετική τ.μ. με παράμετρο\n",
        "  $$\\Lambda_i = -q_{ii},$$\n",
        "- ακολούθως γίνεται άλμα σε νέα κατάσταση $j \\ne i$ με πιθανότητα\n",
        "  \\begin{eqnarray*}\\mathbb{P}(X \\text{ πηγαίνει } i \\to j) =\n",
        "    \\frac{q_{ij}}{\\Lambda_i}.\\end{eqnarray*}\n",
        "\n",
        "Στο κελί αυτό:\n",
        "\n",
        "1. Ορίζουμε μια 3-καταστατική CTMC με γεννήτορα\n",
        "   $$\n",
        "   Q = \\begin{pmatrix}\n",
        "   -1.5 & 1   & 0.5 \\\\\n",
        "   0.5  & -1 & 0.5 \\\\\n",
        "   0.25 & 0.25 & -0.5\n",
        "   \\end{pmatrix}.\n",
        "   $$\n",
        "2. Προσομοιώνουμε μερικές τροχιές $X(t)$ στο διάστημα χρόνου $[0,T]$.\n",
        "3. Παρουσιάζουμε τις τροχιές ως βηματική συνάρτηση (piecewise constant),\n",
        "   ώστε να φανεί η δομή \"παραμονές + άλματα\".\n"
      ],
      "metadata": {
        "id": "T21e0ULVL741"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "def simulate_ctmc(Q, x0=0, T=10.0, rng=None):\n",
        "    \"\"\"\n",
        "    Προσομοίωση μιας CTMC με γεννήτορα Q.\n",
        "    Επιστρέφει πίνακες χρόνων και καταστάσεων (times, states).\n",
        "    \"\"\"\n",
        "    if rng is None:\n",
        "        rng = np.random.default_rng()\n",
        "\n",
        "    t = 0.0\n",
        "    times = [t]\n",
        "    states = [x0]\n",
        "    x = x0\n",
        "    n_states = Q.shape[0]\n",
        "\n",
        "    while t < T:\n",
        "        rate = -Q[x, x]  # Λ_x\n",
        "        if rate <= 0:\n",
        "            # απορροφητική ή χωρίς έξοδο: μένουμε μέχρι το T\n",
        "            times.append(T)\n",
        "            states.append(x)\n",
        "            break\n",
        "\n",
        "        # χρόνος αναμονής ~ Exp(rate)\n",
        "        wait = rng.exponential(1.0 / rate)\n",
        "\n",
        "        if t + wait > T:\n",
        "            # το επόμενο άλμα είναι εκτός του [0,T]\n",
        "            times.append(T)\n",
        "            states.append(x)\n",
        "            break\n",
        "\n",
        "        t = t + wait\n",
        "\n",
        "        # επιλέγουμε επόμενη κατάσταση με πιθανότητες q_{xj} / Λ_x\n",
        "        rates = Q[x].copy()\n",
        "        rates[x] = 0.0\n",
        "        probs = rates / rates.sum()\n",
        "        x = rng.choice(n_states, p=probs)\n",
        "\n",
        "        times.append(t)\n",
        "        states.append(x)\n",
        "\n",
        "    return np.array(times), np.array(states)\n",
        "\n",
        "# Παράδειγμα: 3 τροχιές στο [0, 10]\n",
        "rng = default_rng(123)\n",
        "T = 10.0\n",
        "\n",
        "fig, ax = plt.subplots()\n",
        "\n",
        "for k in range(3):\n",
        "    times, states = simulate_ctmc(Q, x0=0, T=T, rng=rng)\n",
        "    ax.step(times, states, where=\"post\", label=f\"Τροχιά {k+1}\")\n",
        "\n",
        "ax.set_xlabel(\"χρόνος $t$\")\n",
        "ax.set_ylabel(\"κατάσταση $X(t)$\")\n",
        "ax.set_yticks([0, 1, 2])\n",
        "ax.set_title(\"Παράδειγμα τροχιών CTMC\")\n",
        "ax.legend()\n",
        "plt.show()\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 416
        },
        "id": "LfVTe21rMPZP",
        "outputId": "eab636b6-73d4-4d4c-f951-bf43fbca520a"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 800x400 with 1 Axes>"
            ],
            "image/png": "[encoded data, 27716 characters]\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# 10.2 Ακολουθίες πινάκων και σύγκλιση\n",
        "\n",
        "Στις σημειώσεις ορίζεις τη νόρμα\n",
        "$$\n",
        "\\|M\\|_{u,v}\n",
        "= \\sup_{i \\in S} \\sum_{j \\in S} |M_{ij}|,\n",
        "$$\n",
        "και τον χώρο $M_{u,v}(S)$ των πινάκων για τους οποίους η νόρμα αυτή είναι\n",
        "πεπερασμένη. Με αυτή τη νόρμα:\n",
        "\n",
        "- $(M_{u,v}(S), \\|\\cdot\\|_{u,v})$ είναι διανυσματικός χώρος (και στην\n",
        "  πραγματικότητα πλήρης, άρα χώρος Banach),\n",
        "- για $A, B \\in M_{u,v}(S)$ ισχύει\n",
        "  $$\\|AB\\|_{u,v} \\le \\|A\\|_{u,v}\\,\\|B\\|_{u,v}.$$\n",
        "\n",
        "Για πεπερασμένο σύνολο καταστάσεων $S = \\{0,1,\\dots,n-1\\}$ η νόρμα γράφεται\n",
        "$$\n",
        "\\|M\\|_{u,v} = \\max_i \\sum_j |M_{ij}|\n",
        "$$\n",
        "δηλαδή μέγιστο άθροισμα απόλυτων τιμών γραμμής.\n",
        "\n",
        "Στο κελί αυτό:\n",
        "\n",
        "1. Ορίζουμε τη συνάρτηση `matrix_norm_uv`.\n",
        "2. Υπολογίζουμε την ακολουθία μερικών αθροισμάτων\n",
        "   $$S_n(t) = \\sum_{k=0}^n \\frac{t^k}{k!} Q^k,$$\n",
        "   που προσεγγίζει τον πίνακα μετάβασης $P(t)$.\n",
        "3. Δείχνουμε ότι η ακολουθία $(S_n(t))$ είναι Cauchy ως προς τη\n",
        "   $\\|\\cdot\\|_{u,v}$, δηλαδή $\\|S_n - S_{n-1}\\|_{u,v} \\to 0$.\n"
      ],
      "metadata": {
        "id": "lH-AoQrRMSTr"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def matrix_norm_uv(M):\n",
        "    \"\"\"\n",
        "    Νόρμα μέγιστου αθροίσματος γραμμής:\n",
        "    ||M||_{u,v} = max_i sum_j |M_{ij}|\n",
        "    \"\"\"\n",
        "    return np.max(np.sum(np.abs(M), axis=1))\n",
        "\n",
        "def exp_series_partial_sums(Q, t, N_max):\n",
        "    \"\"\"\n",
        "    Επιστρέφει τη λίστα [S_0, S_1, ..., S_N_max]\n",
        "    όπου\n",
        "        S_n(t) = sum_{k=0}^n t^k/k! Q^k.\n",
        "    \"\"\"\n",
        "    n = Q.shape[0]\n",
        "    S_list = []\n",
        "    term = np.eye(n)  # Q^0\n",
        "    S = np.eye(n)     # S_0 = I\n",
        "    S_list.append(S.copy())\n",
        "\n",
        "    for k in range(1, N_max + 1):\n",
        "        # term_k = term_{k-1} * (t/k) * Q\n",
        "        term = term @ Q * (t / k)\n",
        "        S = S + term\n",
        "        S_list.append(S.copy())\n",
        "\n",
        "    return S_list\n",
        "\n",
        "# Παράδειγμα για t=2\n",
        "t = 2.0\n",
        "N_max = 25\n",
        "S_list = exp_series_partial_sums(Q, t, N_max)\n",
        "\n",
        "# Υπολογίζουμε την \"διαφορά\" S_n - S_{n-1}\n",
        "diff_norms = [\n",
        "    matrix_norm_uv(S_list[k] - S_list[k-1]) for k in range(1, len(S_list))\n",
        "]\n",
        "\n",
        "# Εκτύπωση μερικών νόρμων\n",
        "for k in [1, 5, 10, 20, 25]:\n",
        "    if k <= N_max:\n",
        "        print(f\"n = {k}, ||S_n||_u,v = {matrix_norm_uv(S_list[k]):.4f}\")\n",
        "\n",
        "# Γράφημα της ||S_n - S_{n-1}||_u,v\n",
        "fig, ax = plt.subplots()\n",
        "ax.plot(range(1, N_max+1), diff_norms, marker=\"o\")\n",
        "ax.set_yscale(\"log\")\n",
        "ax.set_xlabel(\"δείκτης $n$\")\n",
        "ax.set_ylabel(r\"$\\|S_{n} - S_{n-1}\\|_{u,v}$\")\n",
        "ax.set_title(r\"Σύγκλιση της ακολουθίας $S_n(t)$ για $t = 2$\")\n",
        "plt.show()\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 505
        },
        "id": "YaX73P64L-ik",
        "outputId": "40017dc0-b645-414e-9dc7-f2963196b9f3"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "n = 1, ||S_n||_u,v = 5.0000\n",
            "n = 5, ||S_n||_u,v = 5.4667\n",
            "n = 10, ||S_n||_u,v = 1.0000\n",
            "n = 20, ||S_n||_u,v = 1.0000\n",
            "n = 25, ||S_n||_u,v = 1.0000\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 800x400 with 1 Axes>"
            ],
            "image/png": "[encoded data, 35936 characters]\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# 10.3 Ύπαρξη στοχαστικής διαδικασίας\n",
        "\n",
        "Στις σημειώσεις, στην ενότητα 10.3, δίνεται το Θεώρημα Ύπαρξης:\n",
        "\n",
        "- Δίνονται οι παράμετροι\n",
        "  $$\\Lambda_i \\ge 0$$\n",
        "  και ένας στοχαστικός πίνακας\n",
        "  $$R = (r_{ij}), \\quad \\sum_j r_{ij} = 1.$$\n",
        "- Υποθέτουμε ότι\n",
        "  $$\\sup_{i \\in S} \\Lambda_i < \\infty.$$\n",
        "- Τότε υπάρχει *μοναδική* στοχαστική διαδικασία $X(t)$ σε χρόνο συνεχές,\n",
        "  με Μαρκοβιανή ιδιότητα και κατάλληλες ιδιότητες (συνέχεια από δεξιά,\n",
        "  περιορισμένη πιθανότητα έκρηξης κ.λπ.).\n",
        "\n",
        "Για πεπερασμένο $S$ μπορούμε να ορίσουμε τον γεννήτορα\n",
        "$$\n",
        "Q = \\Lambda (R - I),\n",
        "$$\n",
        "όπου $\\Lambda$ είναι ο διαγώνιος πίνακας με στοιχεία $\\Lambda_i$.\n",
        "\n",
        "Σε αυτό το κελί:\n",
        "\n",
        "1. Ξεκινάμε από δεδομένα $\\Lambda$ και $R$ για 3 καταστάσεις.\n",
        "2. Ορίζουμε τον $Q = \\Lambda (R - I)$.\n",
        "3. Υπολογίζουμε τον πίνακα μετάβασης\n",
        "   $$P(t) \\approx \\sum_{n=0}^N \\frac{t^n}{n!} Q^n.$$\n",
        "4. Ελέγχουμε:\n",
        "   - ότι $P(t)$ είναι (περίπου) στοχαστικός (γραμμές που αθροίζουν στο 1),\n",
        "   - την ιδιότητα ημι-ομάδας $P(t+h) \\approx P(t)P(h)$.\n"
      ],
      "metadata": {
        "id": "wpd_bsvPMcjw"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Ορίζουμε τις παραμέτρους Λ και R\n",
        "Lambda_vals = np.array([1.5, 1.0, 0.5])  # Λ_0, Λ_1, Λ_2\n",
        "\n",
        "R = np.array([\n",
        "    [0.0 , 2/3, 1/3],\n",
        "    [0.5 , 0.0, 0.5],\n",
        "    [0.5 , 0.5, 0.0]\n",
        "])\n",
        "\n",
        "Lambda = np.diag(Lambda_vals)\n",
        "\n",
        "# Γεννήτορας από Q = Λ (R - I)\n",
        "Q_from_LR = Lambda @ (R - np.eye(3))\n",
        "\n",
        "print(\"Λ =\")\n",
        "print(Lambda)\n",
        "print(\"\\nR =\")\n",
        "print(R)\n",
        "print(\"\\nQ = Λ (R - I) =\")\n",
        "print(Q_from_LR)\n",
        "\n",
        "print(\"\\nΣύγκριση με τον Q που χρησιμοποιήσαμε ήδη:\")\n",
        "print(\"max |Q - Q_from_LR| =\",\n",
        "      np.max(np.abs(Q - Q_from_LR)))\n",
        "\n",
        "# Συνάρτηση matrix_exponential_series (θα τη χρησιμοποιήσουμε και σε 10.5, 10.6)\n",
        "def matrix_exponential_series(Q, t, N=40):\n",
        "    \"\"\"\n",
        "    Υπολογίζει προσεγγιστικά το exp(tQ) με ανάπτυξη σε σειρά μέχρι όρο N.\n",
        "    \"\"\"\n",
        "    n = Q.shape[0]\n",
        "    result = np.eye(n)\n",
        "    term = np.eye(n)\n",
        "\n",
        "    for k in range(1, N + 1):\n",
        "        term = term @ Q * (t / k)\n",
        "        result = result + term\n",
        "\n",
        "    return result\n",
        "\n",
        "# Έλεγχος στοχαστικότητας των P(t) για διάφορα t\n",
        "t_values = [0.5, 1.0, 2.0]\n",
        "N_series = 40\n",
        "\n",
        "for t in t_values:\n",
        "    P_t = matrix_exponential_series(Q_from_LR, t, N=N_series)\n",
        "    print(f\"\\n--- t = {t} ---\")\n",
        "    print(\"Άθροισμα γραμμής P(t):\", P_t.sum(axis=1))\n",
        "\n",
        "# Έλεγχος ημι-ομάδας: P(t + h) ≈ P(t) P(h)\n",
        "t, h = 1.0, 0.5\n",
        "P_t = matrix_exponential_series(Q_from_LR, t, N=N_series)\n",
        "P_h = matrix_exponential_series(Q_from_LR, h, N=N_series)\n",
        "P_t_plus_h = matrix_exponential_series(Q_from_LR, t + h, N=N_series)\n",
        "\n",
        "print(\"\\nΈλεγχος P(t+h) ≈ P(t)P(h):\")\n",
        "print(\"||P(t+h) - P(t)P(h)||_u,v =\",\n",
        "      matrix_norm_uv(P_t_plus_h - P_t @ P_h))\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "WH7WMgoWMgec",
        "outputId": "40f0bdd6-d53d-441c-8cba-8c7b90b38bb0"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Λ =\n",
            "[[1.5 0.  0. ]\n",
            " [0.  1.  0. ]\n",
            " [0.  0.  0.5]]\n",
            "\n",
            "R =\n",
            "[[0.         0.66666667 0.33333333]\n",
            " [0.5        0.         0.5       ]\n",
            " [0.5        0.5        0.        ]]\n",
            "\n",
            "Q = Λ (R - I) =\n",
            "[[-1.5   1.    0.5 ]\n",
            " [ 0.5  -1.    0.5 ]\n",
            " [ 0.25  0.25 -0.5 ]]\n",
            "\n",
            "Σύγκριση με τον Q που χρησιμοποιήσαμε ήδη:\n",
            "max |Q - Q_from_LR| = 0.0\n",
            "\n",
            "--- t = 0.5 ---\n",
            "Άθροισμα γραμμής P(t): [1. 1. 1.]\n",
            "\n",
            "--- t = 1.0 ---\n",
            "Άθροισμα γραμμής P(t): [1. 1. 1.]\n",
            "\n",
            "--- t = 2.0 ---\n",
            "Άθροισμα γραμμής P(t): [1. 1. 1.]\n",
            "\n",
            "Έλεγχος P(t+h) ≈ P(t)P(h):\n",
            "||P(t+h) - P(t)P(h)||_u,v = 2.498001805406602e-16\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# 10.4 Η ενσωματωμένη αλυσίδα\n",
        "\n",
        "Για μια CTMC με γεννήτορα $Q$, ορίζουμε:\n",
        "\n",
        "- τις παραμέτρους\n",
        "  $$\\Lambda_i = -q_{ii} \\ge 0,$$\n",
        "- τον διαγώνιο πίνακα\n",
        "  $$\\hat{\\Lambda}_{ij} = \\delta_{ij} \\Lambda_i,$$\n",
        "- τον στοχαστικό πίνακα $R$ ως:\n",
        "  $$\n",
        "  R_{ii} =\n",
        "  \\begin{cases}\n",
        "  1, & \\text{αν } \\Lambda_i = 0,\\\\\n",
        "  0, & \\text{αν } \\Lambda_i > 0,\n",
        "  \\end{cases}\n",
        "  \\qquad\n",
        "  R_{ij} =\n",
        "  \\begin{cases}\n",
        "  0, & \\text{αν } \\Lambda_i = 0,\\\\\n",
        "  \\dfrac{Q_{ij}}{\\Lambda_i}, & \\text{αν } \\Lambda_i > 0 \\text{ και } i \\ne j.\n",
        "  \\end{cases}\n",
        "  $$\n",
        "\n",
        "Ο πίνακας $R$ είναι ο πίνακας μετάβασης της *ενσωματωμένης διακριτής\n",
        "αλυσίδας* $(X_n)$, η οποία «παρακολουθεί» την CTMC στις στιγμές των\n",
        "αλμάτων. Τα άλματα συμβαίνουν σε χρόνους\n",
        "$$\n",
        "J_0 = 0, \\quad\n",
        "J_n = \\inf\\{ t > J_{n-1} : X(t) \\ne X(J_{n-1}) \\},\n",
        "$$\n",
        "και θέτουμε\n",
        "$$\n",
        "X_n = X(J_n).\n",
        "$$\n",
        "\n",
        "Στόχος μας στο κελί αυτό:\n",
        "\n",
        "1. Να κατασκευάσουμε τον πίνακα $R$ από τον $Q$.\n",
        "2. Να επαληθεύσουμε εμπειρικά ότι οι μεταβάσεις της ενσωματωμένης αλυσίδας\n",
        "   ακολουθούν (κατά προσέγγιση) τον $R$.\n",
        "3. Να ελέγξουμε ότι οι χρόνοι αναμονής στην κατάσταση $i$ είναι πράγματι\n",
        "   εκθετικοί με παράμετρο $\\Lambda_i$.\n"
      ],
      "metadata": {
        "id": "R7aLD6NFMk_W"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def build_Lambda_R_from_Q(Q):\n",
        "    \"\"\"\n",
        "    Κατασκευάζει τα Λ_i = -q_{ii} και τον στοχαστικό πίνακα R.\n",
        "    \"\"\"\n",
        "    n = Q.shape[0]\n",
        "    Lambda = -np.diag(Q)\n",
        "    R = np.zeros_like(Q, dtype=float)\n",
        "\n",
        "    for i in range(n):\n",
        "        if Lambda[i] == 0:\n",
        "            R[i, i] = 1.0\n",
        "        else:\n",
        "            for j in range(n):\n",
        "                if i != j:\n",
        "                    R[i, j] = Q[i, j] / Lambda[i]\n",
        "            # Η R[i,i] μπορεί να οριστεί ως 1 - sum_{j≠i} R[i,j]\n",
        "            R[i, i] = max(0.0, 1.0 - np.sum(R[i, :]) + R[i, i])\n",
        "\n",
        "    return Lambda, R\n",
        "\n",
        "Lambda, R_from_Q = build_Lambda_R_from_Q(Q)\n",
        "\n",
        "print(\"Λ =\", Lambda)\n",
        "print(\"\\nR (από τον Q) =\")\n",
        "print(R_from_Q)\n",
        "print(\"Άθροισμα γραμμών του R:\", R_from_Q.sum(axis=1))\n",
        "\n",
        "# Εμπειρική επαλήθευση μέσω προσομοίωσης\n",
        "rng = default_rng(321)\n",
        "T = 20.0\n",
        "n_paths = 2000\n",
        "\n",
        "counts = np.zeros_like(R_from_Q)\n",
        "waiting_times_by_state = {i: [] for i in range(Q.shape[0])}\n",
        "\n",
        "for _ in range(n_paths):\n",
        "    t = 0.0\n",
        "    x = 0  # αρχική κατάσταση\n",
        "    while t < T:\n",
        "        rate = -Q[x, x]\n",
        "        if rate <= 0:\n",
        "            break\n",
        "        wait = rng.exponential(1.0 / rate)\n",
        "        if t + wait > T:\n",
        "            break\n",
        "\n",
        "        # κρατάμε τον χρόνο αναμονής στην κατάσταση x\n",
        "        waiting_times_by_state[x].append(wait)\n",
        "\n",
        "        # επόμενο άλμα\n",
        "        rates = Q[x].copy()\n",
        "        rates[x] = 0.0\n",
        "        probs = rates / rates.sum()\n",
        "        x_next = rng.choice(Q.shape[0], p=probs)\n",
        "\n",
        "        counts[x, x_next] += 1\n",
        "        t = t + wait\n",
        "        x = x_next\n",
        "\n",
        "# Εμπειρικός πίνακας μετάβασης της ενσωματωμένης αλυσίδας\n",
        "emp_R = counts / counts.sum(axis=1, keepdims=True)\n",
        "\n",
        "print(\"\\nΕμπειρικός πίνακας μετάβασης της ενσωματωμένης αλυσίδας:\")\n",
        "print(emp_R)\n",
        "print(\"\\nΔιαφορά emp_R - R_from_Q:\")\n",
        "print(emp_R - R_from_Q)\n",
        "\n",
        "# Έλεγχος εκθετικών χρόνων αναμονής\n",
        "for i in range(Q.shape[0]):\n",
        "    ws = np.array(waiting_times_by_state[i])\n",
        "    if len(ws) > 0:\n",
        "        print(\n",
        "            f\"\\nΚατάσταση {i}: \"\n",
        "            f\"μέσος χρόνος αναμονής ≈ {ws.mean():.3f}, \"\n",
        "            f\"θεωρητικό 1/Λ_i = {1.0 / Lambda[i]:.3f}\"\n",
        "        )\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "5cLsDi0SMwUT",
        "outputId": "840605c2-2636-410e-baf8-6a89ae2f00b7"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Λ = [1.5 1.  0.5]\n",
            "\n",
            "R (από τον Q) =\n",
            "[[0.         0.66666667 0.33333333]\n",
            " [0.5        0.         0.5       ]\n",
            " [0.5        0.5        0.        ]]\n",
            "Άθροισμα γραμμών του R: [1. 1. 1.]\n",
            "\n",
            "Εμπειρικός πίνακας μετάβασης της ενσωματωμένης αλυσίδας:\n",
            "[[0.         0.66642868 0.33357132]\n",
            " [0.49463391 0.         0.50536609]\n",
            " [0.4982251  0.5017749  0.        ]]\n",
            "\n",
            "Διαφορά emp_R - R_from_Q:\n",
            "[[ 0.         -0.00023798  0.00023798]\n",
            " [-0.00536609  0.          0.00536609]\n",
            " [-0.0017749   0.0017749   0.        ]]\n",
            "\n",
            "Κατάσταση 0: μέσος χρόνος αναμονής ≈ 0.656, θεωρητικό 1/Λ_i = 0.667\n",
            "\n",
            "Κατάσταση 1: μέσος χρόνος αναμονής ≈ 0.929, θεωρητικό 1/Λ_i = 1.000\n",
            "\n",
            "Κατάσταση 2: μέσος χρόνος αναμονής ≈ 1.795, θεωρητικό 1/Λ_i = 2.000\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# 10.5 Εξισώσεις Chapman–Kolmogorov\n",
        "\n",
        "Για μια CTMC με πίνακες μετάβασης\n",
        "$$\n",
        "P(t) = (p_{ij}(t))_{i,j}, \\quad\n",
        "p_{ij}(t) = \\mathbb{P}(X(t) = j \\mid X(0) = i),\n",
        "$$\n",
        "οι εξισώσεις Chapman–Kolmogorov γράφονται:\n",
        "$$\n",
        "P(t+s) = P(t)\\,P(s), \\quad t,s \\ge 0.\n",
        "$$\n",
        "\n",
        "Επιπλέον, για τον γεννήτορα $Q$ ισχύει ότι\n",
        "$$\n",
        "P'(t) = P(t)Q = QP(t), \\quad P(0) = I,\n",
        "$$\n",
        "οπότε\n",
        "$$\n",
        "Q = P'(0) = \\lim_{h \\to 0} \\frac{P(h) - I}{h}.\n",
        "$$\n",
        "\n",
        "Στο κελί αυτό:\n",
        "\n",
        "1. Υπολογίζουμε προσεγγιστικά τα $P(t)$ με χρήση της ανάπτυξης σε σειρά\n",
        "   $$P(t) \\approx \\sum_{n=0}^N \\frac{t^n}{n!} Q^n.$$\n",
        "2. Ελέγχουμε αριθμητικά την εξίσωση\n",
        "   $$P(t+s) \\approx P(t)\\,P(s).$$\n",
        "3. Υπολογίζουμε για μικρά $h$ την ποσότητα\n",
        "   $$\\frac{P(h) - I}{h},$$\n",
        "   και συγκρίνουμε με τον $Q$ για να προσεγγίσουμε το $P'(0)$.\n"
      ],
      "metadata": {
        "id": "s3KUXSuQNxDB"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Χρησιμοποιούμε την matrix_exponential_series και τη matrix_norm_uv που ήδη ορίσαμε\n",
        "\n",
        "# Chapman–Kolmogorov: P(t+s) ≈ P(t)P(s)\n",
        "t, s = 1.25, 0.75\n",
        "P_t = matrix_exponential_series(Q, t, N=40)\n",
        "P_s = matrix_exponential_series(Q, s, N=40)\n",
        "P_t_plus_s = matrix_exponential_series(Q, t + s, N=40)\n",
        "\n",
        "diff_CK = matrix_norm_uv(P_t_plus_s - P_t @ P_s)\n",
        "print(\"Έλεγχος εξίσωσης Chapman–Kolmogorov:\")\n",
        "print(f\"||P(t+s) - P(t)P(s)||_u,v ≈ {diff_CK:.3e}\")\n",
        "\n",
        "# Προσέγγιση του Q μέσω (P(h) - I)/h για μικρά h\n",
        "h_values = [1e-1, 5e-2, 1e-2, 5e-3]\n",
        "\n",
        "print(\"\\nΠροσέγγιση Q = P'(0) μέσω (P(h)-I)/h:\")\n",
        "for h in h_values:\n",
        "    P_h = matrix_exponential_series(Q, h, N=40)\n",
        "    approx_Q = (P_h - np.eye(Q.shape[0])) / h\n",
        "    err = matrix_norm_uv(approx_Q - Q)\n",
        "    print(f\"h = {h:.3e}, ||(P(h)-I)/h - Q||_u,v ≈ {err:.3e}\")\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "ixH8g-XAN1Ij",
        "outputId": "4db7297e-7888-4742-f9e0-9ef21b6b4692"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Έλεγχος εξίσωσης Chapman–Kolmogorov:\n",
            "||P(t+s) - P(t)P(s)||_u,v ≈ 1.055e-15\n",
            "\n",
            "Προσέγγιση Q = P'(0) μέσω (P(h)-I)/h:\n",
            "h = 1.000e-01, ||(P(h)-I)/h - Q||_u,v ≈ 2.696e-01\n",
            "h = 5.000e-02, ||(P(h)-I)/h - Q||_u,v ≈ 1.392e-01\n",
            "h = 1.000e-02, ||(P(h)-I)/h - Q||_u,v ≈ 2.856e-02\n",
            "h = 5.000e-03, ||(P(h)-I)/h - Q||_u,v ≈ 1.433e-02\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# 10.6 Στάσιμη κατανομή – Οριακές πιθανότητες\n",
        "\n",
        "Για μια CTMC με γεννήτορα $Q$ μια κατανομή\n",
        "$$\\pi = (\\pi_i)_{i \\in S}$$\n",
        "λέγεται *στάσιμη* αν\n",
        "$$\n",
        "\\pi Q = 0, \\quad \\sum_i \\pi_i = 1, \\quad \\pi_i \\ge 0.\n",
        "$$\n",
        "\n",
        "Υπό κατάλληλες υποθέσεις (αδιαχωριστότητα, θετική επαναληπτικότητα κ.λπ.)\n",
        "η στάσιμη κατανομή δίνει και τις *οριακές πιθανότητες*:\n",
        "$$\n",
        "\\lim_{t \\to \\infty} \\mathbb{P}(X(t) = j) = \\pi_j.\n",
        "$$\n",
        "\n",
        "Σε αυτό το κελί:\n",
        "\n",
        "1. Υπολογίζουμε τη στάσιμη κατανομή $\\pi$ λύνοντας το γραμμικό σύστημα\n",
        "   $$\\pi Q = 0, \\quad \\sum_i \\pi_i = 1.$$\n",
        "2. Υπολογίζουμε $P(t)$ για μεγάλα $t$ και ελέγχουμε ότι οι γραμμές\n",
        "   του $P(t)$ συγκλίνουν προς $\\pi$.\n",
        "3. Προσομοιώνουμε μια μακρά τροχιά της CTMC και εκτιμούμε την κατανομή\n",
        "   χρόνου παραμονής σε κάθε κατάσταση, συγκρίνοντάς την με τη $\\pi$.\n"
      ],
      "metadata": {
        "id": "1R7lw0yVN6TW"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def stationary_distribution_ctmc(Q):\n",
        "    \"\"\"\n",
        "    Υπολογίζει στάσιμη κατανομή π για CTMC επιλύοντας:\n",
        "        π Q = 0, sum_i π_i = 1.\n",
        "    Χρησιμοποιεί την αντικατάσταση μιας εξίσωσης με sum π_i = 1.\n",
        "    \"\"\"\n",
        "    n = Q.shape[0]\n",
        "    A = Q.T.copy()\n",
        "    b = np.zeros(n)\n",
        "\n",
        "    # Αντικαθιστούμε την τελευταία γραμμή με την εξίσωση sum π_i = 1\n",
        "    A[-1, :] = 1.0\n",
        "    b[-1] = 1.0\n",
        "\n",
        "    pi = np.linalg.solve(A, b)\n",
        "    return pi\n",
        "\n",
        "pi = stationary_distribution_ctmc(Q)\n",
        "print(\"Στάσιμη κατανομή π:\")\n",
        "print(pi)\n",
        "print(\"Άθροισμα π:\", pi.sum())\n",
        "print(\"π Q (πρέπει να είναι ≈ 0):\", pi @ Q)\n",
        "\n",
        "# Έλεγχος σύγκλισης P(t) προς π (γραμμές)\n",
        "t_values = [1.0, 3.0, 5.0, 10.0, 20.0]\n",
        "print(\"\\nΣύγκριση γραμμών του P(t) με π για μεγάλα t:\\n\")\n",
        "\n",
        "for t in t_values:\n",
        "    P_t = matrix_exponential_series(Q, t, N=60)\n",
        "    # π.χ. η 0-στή γραμμή\n",
        "    row0 = P_t[0, :]\n",
        "    diff = row0 - pi\n",
        "    print(f\"t = {t:>5.1f}, row_0(P(t)) =\", row0, \"  διαφορά row_0 - π =\", diff)\n",
        "\n",
        "# Προσομοίωση μακράς τροχιάς και εμπειρική κατανομή χρόνου παραμονής\n",
        "def empirical_time_distribution(Q, x0, T, rng=None):\n",
        "    \"\"\"\n",
        "    Προσομοιώνει την CTMC μέχρι χρόνο T και επιστρέφει τις συνολικές\n",
        "    διάρκειες παραμονής σε κάθε κατάσταση.\n",
        "    \"\"\"\n",
        "    if rng is None:\n",
        "        rng = np.random.default_rng()\n",
        "    n = Q.shape[0]\n",
        "    times_in_state = np.zeros(n)\n",
        "\n",
        "    t = 0.0\n",
        "    x = x0\n",
        "\n",
        "    while t < T:\n",
        "        rate = -Q[x, x]\n",
        "        if rate <= 0:\n",
        "            times_in_state[x] += (T - t)\n",
        "            break\n",
        "\n",
        "        wait = rng.exponential(1.0 / rate)\n",
        "        if t + wait > T:\n",
        "            times_in_state[x] += (T - t)\n",
        "            break\n",
        "\n",
        "        times_in_state[x] += wait\n",
        "        t = t + wait\n",
        "\n",
        "        # επόμενο άλμα\n",
        "        rates = Q[x].copy()\n",
        "        rates[x] = 0.0\n",
        "        probs = rates / rates.sum()\n",
        "        x = rng.choice(n, p=probs)\n",
        "\n",
        "    return times_in_state\n",
        "\n",
        "rng = default_rng(999)\n",
        "T_total = 2000.0  # μακρύς συνολικός χρόνος\n",
        "times_in_state = empirical_time_distribution(Q, x0=0, T=T_total, rng=rng)\n",
        "emp_pi = times_in_state / times_in_state.sum()\n",
        "\n",
        "print(\"\\nΕμπειρική κατανομή χρόνου παραμονής (μακρά τροχιά):\")\n",
        "print(\"empirical π ≈\", emp_pi)\n",
        "print(\"θεωρητικό π =\", pi)\n",
        "print(\"διαφορά empirical - θεωρητικό =\", emp_pi - pi)\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "4swmMkbSN9Yq",
        "outputId": "26805e9c-7e7c-4f23-fb9c-52c737bc6750"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Στάσιμη κατανομή π:\n",
            "[0.1875 0.3125 0.5   ]\n",
            "Άθροισμα π: 1.0\n",
            "π Q (πρέπει να είναι ≈ 0): [0. 0. 0.]\n",
            "\n",
            "Σύγκριση γραμμών του P(t) με π για μεγάλα t:\n",
            "\n",
            "t =   1.0, row_0(P(t)) = [0.32652794 0.35741178 0.31606028]   διαφορά row_0 - π = [ 0.13902794  0.04491178 -0.18393972]\n",
            "t =   3.0, row_0(P(t)) = [0.19542753 0.32946601 0.47510647]   διαφορά row_0 - π = [ 0.00792753  0.01696601 -0.02489353]\n",
            "t =   5.0, row_0(P(t)) = [0.18837346 0.31499552 0.49663103]   διαφορά row_0 - π = [ 0.00087346  0.00249552 -0.00336897]\n",
            "t =  10.0, row_0(P(t)) = [0.18752927 0.31249343 0.4999773 ]   διαφορά row_0 - π = [ 2.92714977e-05 -6.57276907e-06 -2.26998868e-05]\n",
            "t =  20.0, row_0(P(t)) = [ 4.36655189e+13 -4.36655189e+13  1.09104383e+00]   διαφορά row_0 - π = [ 4.36655189e+13 -4.36655189e+13  5.91043828e-01]\n",
            "\n",
            "Εμπειρική κατανομή χρόνου παραμονής (μακρά τροχιά):\n",
            "empirical π ≈ [0.20474785 0.28471295 0.5105392 ]\n",
            "θεωρητικό π = [0.1875 0.3125 0.5   ]\n",
            "διαφορά empirical - θεωρητικό = [ 0.01724785 -0.02778705  0.0105392 ]\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [],
      "metadata": {
        "id": "zr-qpwMWODNv"
      }
    }
  ]
}