Py8.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### Universidad de Costa Rica\n",
    "#### IE0405 - Modelos Probabilísticos de Señales y Sistemas\n",
    "---\n",
    "\n",
    "# `Py8` - *Simulaciones aleatorias*\n",
    "\n",
    "> Conocidas como **simulaciones de Monte Carlo**, las simulaciones o experimentos aleatorios repiten una acción utilizando datos aleatorios para analizar sus resultados. El objetivo es conocer el comportamiento de un fenómeno o sistema para el cual se conoce el modelo de probabilidad de sus entradas y el modelo propio.\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introducción\n",
    "\n",
    "En contextos de sistemas de control o sistemas de comunicaciones es usual este tipo de simulaciones para probar estrategias de sintonización o esquemas de modulación digital."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## 8.1 - Simulaciones aleatorias aplicadas a transformaciones de variables aleatorias\n",
    "\n",
    "Aun sin conocer una expresión para la función de densidad $f_Y(y)$ de una variable aleatoria $Y$ producto de una transformación $Y = T(X)$, es posible conocer el efecto que tiene esta transformación sobre una muestra de datos aleatorios con una distribución de probabilidad $X$ conocida.\n",
    "\n",
    "### Ejemplo de una distribución normal con una transformación cuadrática\n",
    "\n",
    "A una distribución $X \\sim \\mathcal{N}(0, 0.7737)$ se aplica una transformación $Y = X^2$. Es útil visualizar los dos conjuntos de datos."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import stats\n",
    "import matplotlib.pyplot as plt\n",
    "import csv\n",
    "\n",
    "# Repeticiones\n",
    "N = 500\n",
    "\n",
    "# Distribución de X\n",
    "Xva = stats.truncnorm(-2,2)\n",
    "\n",
    "# Muestra de N valores aleatorios\n",
    "X = Xva.rvs(N)\n",
    "\n",
    "# Función aplicada sobre X\n",
    "Y = X**2\n",
    "\n",
    "# Crear histograma de X\n",
    "a,b,c = plt.hist(X,20)\n",
    "plt.show()\n",
    "\n",
    "# Crear histograma de Y\n",
    "a,b,c = plt.hist(Y,20)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## 8.2 - Ejemplo del canal de comunicaciones AWGN con la modulación OOK\n",
    "\n",
    "La modulación **OOK** (*on/off keying*, codificación por encendido y apagado) consiste en enviar una forma de onda sinusoidal para el bit 1, y no enviar nada para el bit 0. Esta es una modulación binaria de un bit por símbolo.\n",
    "\n",
    "#### Generación de los bits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import stats\n",
    "\n",
    "# Número de bits\n",
    "N = 10000\n",
    "\n",
    "# Variable aleatoria binaria equiprobable\n",
    "X = stats.bernoulli(0.5)\n",
    "\n",
    "# Generar bits para \"transmitir\"\n",
    "bits = X.rvs(N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Generar la señal modulada OOK\n",
    "\n",
    "Los parámetros son:\n",
    "\n",
    "* Frecuencia de la portadora: $f =$ 1000 Hz\n",
    "* Tiempo del símbolo: $T_s = T = 1/f$ (es decir, un período del símbolo es un período completo de la onda portadora)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Frecuencia de operación\n",
    "f = 1000 # Hz\n",
    "\n",
    "# Duración del período de cada símbolo (onda)\n",
    "T = 1/f # 1 ms\n",
    "\n",
    "# Número de puntos de muestreo por período\n",
    "p = 50\n",
    "\n",
    "# Puntos de muestreo para cada período\n",
    "tp = np.linspace(0, T, p)\n",
    "\n",
    "# Creación de la forma de onda de la portadora\n",
    "sinus = np.sin(2*np.pi * f * tp)\n",
    "\n",
    "# Visualización de la forma de onda de la portadora\n",
    "plt.plot(tp, sinus)\n",
    "plt.xlabel('Tiempo / s')\n",
    "plt.show()\n",
    "\n",
    "# Frecuencia de muestreo\n",
    "fs = p/T # 50 kHz\n",
    "\n",
    "# Creación de la línea temporal para toda la señal Tx\n",
    "t = np.linspace(0, N*T, N*p)\n",
    "\n",
    "# Inicializar el vector de la señal modulada Tx\n",
    "senal = np.zeros(t.shape)\n",
    "\n",
    "# Creación de la señal modulada OOK\n",
    "for k, b in enumerate(bits):\n",
    "    senal[k*p:(k+1)*p] = b * sinus\n",
    "\n",
    "# Visualización de los primeros bits modulados\n",
    "pb = 5\n",
    "plt.figure()\n",
    "plt.plot(senal[0:pb*p])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Cálculo de la potencia promedio\n",
    "\n",
    "Es conocido que:\n",
    "\n",
    "$$\n",
    "P(T) = \\frac{1}{2T}\\int_{-T}^{T}x^2(t) ~\\mathrm{d}t = A\\{x^2(t)\\}\n",
    "$$\n",
    "\n",
    "y también que:\n",
    "\n",
    "$$\n",
    "P(T) = \\frac{1}{2T}\\int_{-T}^{T}E[X^2(t)] ~\\mathrm{d}t = A\\{E[X^2(t)]\\}\n",
    "$$\n",
    "\n",
    "donde $A\\{\\cdot\\}$ es el operador de promedio temporal y $E[X^2(t)]$ es la autocorrelación del proceso aleatorio $X(t)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import integrate\n",
    "\n",
    "# Potencia instantánea\n",
    "Pinst = senal**2\n",
    "\n",
    "# Potencia promedio a partir de la potencia instantánea (W)\n",
    "Ps = integrate.trapz(Pinst, t) / (N * T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Simulación de un canal ruidoso AWGN\n",
    "\n",
    "Conociendo que la relación señal-a-ruido es\n",
    "\n",
    "$$\n",
    "SNR_{dB} = 10 \\log_{10} \\left( \\frac{P_s}{P_n} \\right)\n",
    "$$\n",
    "\n",
    "donde $P_s$ es la potencia de la señal y $P_n$ es la potencia del ruido ($n$ de *noise*)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Relación señal-a-ruido deseada\n",
    "SNR = 3\n",
    "\n",
    "# Potencia del ruido para SNR y potencia de la señal dadas\n",
    "Pn = Ps / (10**(SNR / 10))\n",
    "\n",
    "# Desviación estándar del ruido\n",
    "sigma = np.sqrt(Pn)\n",
    "\n",
    "# Crear ruido (Pn = sigma^2)\n",
    "ruido = np.random.normal(0, sigma, senal.shape)\n",
    "\n",
    "# Simular \"el canal\": señal recibida\n",
    "Rx = senal + ruido\n",
    "\n",
    "# Visualización de los primeros bits recibidos\n",
    "pb = 5\n",
    "plt.figure()\n",
    "plt.plot(Rx[0:pb*p])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Densidad espectral de potencia"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import signal\n",
    "\n",
    "# Antes del canal ruidoso\n",
    "fw, PSD = signal.welch(senal, fs, nperseg=1024)\n",
    "plt.figure()\n",
    "plt.semilogy(fw, PSD)\n",
    "plt.xlabel('Frecuencia / Hz')\n",
    "plt.ylabel('Densidad espectral de potencia / V**2/Hz')\n",
    "plt.show()\n",
    "\n",
    "# Después del canal ruidoso\n",
    "fw, PSD = signal.welch(Rx, fs, nperseg=1024)\n",
    "plt.figure()\n",
    "plt.semilogy(fw, PSD)\n",
    "plt.xlabel('Frecuencia / Hz')\n",
    "plt.ylabel('Densidad espectral de potencia / V**2/Hz')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Demodulación y decodificación de la señal\n",
    "\n",
    "Es útil aquí el conocimiento del **producto interno** de dos señales, utilizado frecuentemente para modulación *coherente* (donde la fase, es decir, la sincronización, es perfecta). Se denota como:\n",
    "\n",
    "$$\n",
    "\\langle g(t), h(t) \\rangle = \\int_0^T g(t) h(t) ~\\mathrm{d}t\n",
    "$$\n",
    "\n",
    "Y aplica que:\n",
    "\n",
    "* $\\langle g(t), g(t) \\rangle = E_g$, la energía de la señal en el intervalo de integración.\n",
    "* $\\langle g(t), h(t) \\rangle = 0$ cuando $g(t)$ y $h(t)$ son ortogonales."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pseudo-energía de la onda original (esta es suma, no integral)\n",
    "Es = np.sum(sinus**2)\n",
    "\n",
    "# Inicialización del vector de bits recibidos\n",
    "bitsRx = np.zeros(bits.shape)\n",
    "\n",
    "# Decodificación de la señal por detección de energía\n",
    "for k, b in enumerate(bits):\n",
    "    Ep = np.sum(Rx[k*p:(k+1)*p] * sinus)\n",
    "    if Ep > Es/2:\n",
    "        bitsRx[k] = 1\n",
    "    else:\n",
    "        bitsRx[k] = 0\n",
    "\n",
    "err = np.sum(np.abs(bits - bitsRx))\n",
    "BER = err/N\n",
    "\n",
    "print('Hay un total de {} errores en {} bits para una tasa de error de {}.'.format(err, N, BER))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8.3 - Ejemplo de la aproximación del número $\\pi$\n",
    "\n",
    "(Pendiente)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8.4 - Ejemplo de un proceso de Markov\n",
    "\n",
    "Con una tasa de llegada $\\lambda$ y un parámetro de tiempo de servicio $\\nu$, un proceso de Markov con un solo servidor también recibe la notación de *sistema de colas* **M/M/1** (donde la **M** viene de *Markov*). Más en general, con $s$ servidores es un sistema **M/M/s**.\n",
    "\n",
    "La simulación de un sistema del tipo **M/M/1** implica la generación de una llegada de \"clientes\" como una *corriente de Poisson*. Esto es equivalente a decir que tienen una distribución de probabilidad de *tiempo entre arribos* con distribución exponencial y parámetro $\\lambda$.\n",
    "\n",
    "Por su parte, el *tiempo de servicio* tiene también una distribución exponencial pero con parámetro $\\nu$. A la relación $\\lambda/\\nu$ usualmente se le conoce como $\\rho$.\n",
    "\n",
    "**Nota**: El tiempo de servicio se asume independiente del tiempo de llegada.\n",
    "\n",
    "#### Sobre la simulación\n",
    "\n",
    "Es posible crear una simulación de $N$ clientes con sus respectivos tiempos de servicio, distribuidos en el tiempo.\n",
    "\n",
    "Para medir el tiempo se puede utilizar una medida mínima arbitraria, pero que tenga sentido para el problema. Por ejemplo: si la tasa de llegada es de 1 persona/minuto, vale más tener una granularidad de segundos o decenas de segundos para capturar una precisión temporal adecuada.\n",
    "\n",
    "**Nota**: La selección de $N$ puede depender de la precisión deseada para el resultado. Por ejemplo: con $N = 1000$ es posible obtener una precisión de hasta el 0,1% ($1/N$).\n",
    "\n",
    "### Problema de ejemplo: un servidor web\n",
    "\n",
    "> Un servidor web es modelado como un sistema M/M/1 con una tasa de arribo de 2 solicitudes por minuto. Es deseado tener 4 o menos solicitudes en fila el 99\\% del tiempo. ¿Qué tan rápido debe ser el servicio? $\\nu$ es solicitudes atendidas por minuto.\n",
    "\n",
    "El estado $i$ es el número de clientes en el sistema. La longitud de la fila es $L_q = i - 1$ (en virtud de la solicitud que está siendo atendida en $s = 1$ servidores). Es posible encontrar que:\n",
    "\n",
    "$$\n",
    "P( \\text{5 o más clientes en el sistema} ) = \\sum_{i=5}^{\\infty} (1 - \\rho) \\rho^i  = 1 - \\sum_{i=0}^{4} (1 - \\rho) \\rho^i = \\rho^5\n",
    "$$\n",
    "\n",
    "que depende de $\\rho = \\lambda/\\nu$ y del parámetro de servicio $\\nu$ buscado. \n",
    "\n",
    "De los datos del problema: $\\lambda = 2$. Para tener una fila de 3 o menos clientes el 99\\% del tiempo se necesita:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "P( \\text{5 o más clientes en el sistema} ) = \\rho^5 & = \\left( \\frac{\\lambda}{\\nu} \\right)^5 \\leq 0.01 \\\\\n",
    "\\nu^5 & \\geq \\frac{\\lambda^5}{0.01} = \\frac{2^5}{0.01} = 3200 \\quad \\Rightarrow \\quad \\nu \\geq 5.024\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "es decir, el servidor debe atender más de 5,024 solicitudes por minuto en promedio para poder satisfacer el requisito.\n",
    "\n",
    "En la siguiente simulación, con $N = 1000$, y $\\nu = 3 < 5.024$ deberíamos obtener una probabilidad $P( \\text{5 o más clientes en el sistema} ) > 0.01$ que **no** cumple con las especificaciones.\n",
    "\n",
    "**Nota**: Observar el cambio de unidades de minutos a segundos, para lograr mayor \"granularidad\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import stats\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Número de clientes\n",
    "N = 1000\n",
    "\n",
    "# Parámetro de llegada (clientes/segundos)\n",
    "lam = 2/60\n",
    "\n",
    "# Parámetro de servicio (servicios/segundos)\n",
    "nu = 3/60\n",
    "\n",
    "# Distribución de los tiempos de llegada entre cada cliente\n",
    "X = stats.expon(scale = 1/lam)\n",
    "\n",
    "# Distribución de los tiempos de servicio a cada cliente\n",
    "Y = stats.expon(scale = 1/nu)\n",
    "\n",
    "# Intervalos entre llegadas (segundos desde último cliente)\n",
    "t_inte = np.ceil(X.rvs(N)).astype('int')\n",
    "\n",
    "# Tiempos de las llegadas (segundos desde el inicio)\n",
    "t_lleg = [t_inte[0]]\n",
    "for i in range(1, len(t_inte)):\n",
    "    siguiente = t_lleg[i-1] + t_inte[i]\n",
    "    t_lleg.append(siguiente)\n",
    "\n",
    "# Tiempos de servicio (segundos desde inicio de servicio)\n",
    "t_serv = np.ceil(Y.rvs(N)).astype('int')\n",
    "\n",
    "# Inicialización del tiempo de inicio y fin de atención\n",
    "inicio = t_lleg[0]          # primera llegada\n",
    "fin = inicio + t_serv[0]    # primera salida\n",
    "\n",
    "# Tiempos en que recibe atención cada i-ésimo cliente (!= que llega)\n",
    "t_aten = [inicio]\n",
    "for i in range(1, N):\n",
    "    inicio = np.max((t_lleg[i], fin))\n",
    "    fin = inicio + t_serv[i]\n",
    "    t_aten.append(inicio)\n",
    "\n",
    "# Inicialización del vector temporal para registrar eventos\n",
    "t = np.zeros(t_aten[-1] + t_serv[-1] + 1)\n",
    "\n",
    "# Asignación de eventos de llegada (+1) y salida (-1) de clientes\n",
    "for c in range(N):\n",
    "    i = t_lleg[c]\n",
    "    t[i] += 1\n",
    "    j = t_aten[c] + t_serv[c]\n",
    "    t[j] -= 1\n",
    "\n",
    "# Umbral de P o más personas en sistema (hay P - 1 en fila)\n",
    "P = 5\n",
    "\n",
    "# Instantes (segundos) de tiempo con P o más solicitudes en sistema\n",
    "frecuencia = 0\n",
    "\n",
    "# Proceso aleatorio (estados n = {0, 1, 2...})\n",
    "Xt = np.zeros(t.shape)\n",
    "\n",
    "# Inicialización de estado n\n",
    "n = 0\n",
    "\n",
    "# Recorrido del vector temporal y conteo de clientes (estado n)\n",
    "for i, c in enumerate(t):\n",
    "    n += c # sumar (+1) o restar (-1) al estado\n",
    "    Xt[i] = n\n",
    "    if Xt[i] >= P: \n",
    "        frecuencia += 1\n",
    "\n",
    "# Fracción de tiempo con P o más solicitudes en sistema\n",
    "fraccion = frecuencia / len(t)\n",
    "\n",
    "# Resultados\n",
    "print('Parámetro lambda = ', str(lam*60))\n",
    "print('Parámetro nu = ', str(nu*60))\n",
    "print('Tiempo con más de {} solicitudes en fila:'.format(P-2))\n",
    "print('\\t {:0.2f}%'.format(100*fraccion))\n",
    "if fraccion <= 0.01:\n",
    "    print('\\t Sí cumple con la especificación.')\n",
    "else:\n",
    "    print('\\t No cumple con la especificación.') \n",
    "print('Simulación es equivalente a {:0.2f} horas.'.format(len(t)/3600))\n",
    "\n",
    "# Gráfica de X(t) (estados del sistema)\n",
    "plt.figure()\n",
    "plt.plot(Xt)\n",
    "plt.plot(range(len(t)), (P-1)*np.ones(t.shape))\n",
    "plt.legend(('$X(t) = n$', '$L_q = $' + str(P-2)))\n",
    "plt.ylabel('Clientes en el sistema, $n$')\n",
    "plt.xlabel('Tiempo, $t$ / segundos')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### Más información\n",
    "\n",
    "* Función [scipy.integrate.trapz](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.trapz.html)\n",
    "* Modulación [OOK](https://en.wikipedia.org/wiki/On%E2%80%93off_keying)\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "**Universidad de Costa Rica**\n",
    "\n",
    "Facultad de Ingeniería\n",
    "\n",
    "Escuela de Ingeniería Eléctrica\n",
    "\n",
    "---"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}