Un proceso de decisión markoviano (MDP) es un modelo secuencial de decisiones hechas bajo incertidumbre, el cual toma en cuenta los resultados a corto tiempo dadas las deciciones tomadas y las oportunidad de hacer decisiones en el futuro. En concreto el término MDP fue introducido por Bellman para hacer una descipción de un proceso estocástico controlado por una secuencia de acciones bajo condiciones de incertidumbre (Udom, 2013). Además, un MDP debe de tener 7 caracterísitcas (Kallenberg, 2016):
El espacio de estados (N): conjunto de estados posibles.
Conjunto de acciones ( \(A_i\) ): Dado que el agente observa que el sistema está en un estado i (\(S_{i}\)), entonces este debe elegir una acción entre el conunto de acciones que dependen del estado observado.
Tiempos de decisión: El intervalo temporal entre las decisiones debe ser constante (MPD) o al azar (semi-MPD). -Recompenzas/costos inmediatos: Dado el estado en el que se encuentra el sistema y la acción elegida por el agente, después una recompensa (o costo) inmedianta es obtenida. Estas recomenpensas dependen del estado observado, la acción tomada y del punto decisional en el tiempo, pero no de la historia del proceso. -Transiciones de probabilidades: Dado el estado del sistema y la acción elegida, el siguiente estado en el próximo punto decisional en el tiempo es determinado por una ley de transición (propiedad markoviana).
Horizonte de planeación: Resultado de los puntos en el tiempo en que el sistema es controlado. -Criterio óptimo: El objetivo de un problema de decisión Markoviano es determinar una política (regla decisional) para cada punto decisional en el tiempo que optimice el redimiento del sistema. Esto es medido por una función de utilidad, la cual asigna a cada política un valor tomando en cuenta el estado inicial.
Un proceso markoviano parcialmente observable (POMDP) es una combinación de un MDP regular y las dinámicas del sistema con un modelo Markoviano oculto, ya que se conectan los estados probabilísticos no observables del sistema con las observaciones (Hahsler y Kamalzadeh, 2021). Es decir, el agente no puede directamente observar el estado completo del sistema (por lo que se usan probabilidades de observación), pero puede hacer observaciones que dependen del estado y el agente usa estas observaciones para crearse una creencia sobre en qué estado está el sistema en ese momento. Esta creencia se llama estado de creencia y es expresado como una distribución de probabilidad sobre todos los posibles estados (Hahsler et al., 2021). La solución de un POMDP es una política que prescriba qué acción se debe tomar en cada estado de creencia.
Cuando se trabaja con distribuciones de probabilidad y actualizaciones de estas probabilidades la estadísitica bayesiana es la más últil para abordar el problema, ya que la idea de esta es usar las distrubuciones de probabilidad no sólo en las variables observables, como en la estadística frecuentista, sino también en los parámetros que generan esas variables observadas (Verguts, 2022). De hecho la inferencia bayesiana es una técnica general para el aprendizaje de parámetros desconocidos de un modelo de probabilidad, dadas las observaciones de ese modelo. En concreto la probabilidad bayesiana está basada en el Teorema de Bayes que conecta el grado de creencia en una crencia previa (prior) y una después (posterior) de que se ha recopilado evidencia, tomando en cuenta las probabilidades condicionales implicadas (Lin, 2013):
\[P(p|X) = \frac{P(X|p) P(p)}{P(X)}
\] Donde: - P ( p | X ) = Posterior - P ( p ) = Prior - P ( X | P ) = Verosimilitud (probabilidad de X dado p) - P ( X ) = Evidencia (toma en cuenta toda la evidencia de X)
El control óptimo es un aspecto de la optimización en la que los parámetros de entrada (control) de un sistem dinámico son manipulados para alcanzar los resultados deseados, ya sea minimizando la función de costo o maximizando la función de ganancia asociada con la trjectoria de control del sistema (Udom, 2021). Los MDP son un puente entre el control óptimo estocástico y el control óptimo determinístico, tanto así que los criterios óptimos más usandos en los MDP son la minimización del horizonte finito de costo esperado o criterios relacionados a este (Udom, 2021).
En la teoría de control óptimo la creencia es es una probabilidad posterior sobre la variable latente dadas todas las medidas anteriores, por lo que se considera que maximizar la utilidad esperada con respecto a la posterior es óptimo. Además, se puede parametrizar esta función con un umbral, es así que cuando tu creencia cae debajo de un umbral entonces se toma una acción diferente (Wu, Saxena, y Pitkow, 2021).
El sigueinte simulador se hizo con código adaptado de un tutorial de Neuromacth Academy y sólo se modificó para que en una sola gráfica estuvieran todos los sliders y así se entendiera mejor cómo afecta la modificación de cada parámetro (Wu, Saxena, y Pitkow, 2021).
Este simulador sirve para estudiar en qué consiste un POMDP e ideas básicas del control óptimo, donde se pueden modificar parámetros como: - La probabilidad de cambio de estado (Prob_quedarse), cuando se tiene un espacio de estados binario. - El umbral para decidir qué probabilidad se necesita tener del estado de creencia para seleccionar la acción de cambiar o quedarse (Umbral). - Probabilidad de obtener una recompensa dado que no se tomó la acción adecuada (Prob_pesca_LO). - Probabilidad de obtener una recompensa dado que se tomó la acción adecuada (Prob_pesca_ML).
El problema que se platea es el siguiente:
Un hombre fue a pescar a un lago que se divide por un puente de madera (izquierda-derecha), es decir tiene sólo dos sitios dónde pescar y obtener una recompensa (peces). Además, él no sabe dónde está el cardúmen (banco de peces) y no puede ver a través del agua (representa un modelo markoviano oculto porque hay dinámicas ocultas). Basándose en cuándo y dónde capturó peces, el lector puede actualizar su creencia sobre la localización de los peces (están a la derecha o a la izquierda): se actualiza la posterior de dónde están los peces dadas las observaciones anteriores. El lector deberá controlar su posición para obtener la mayor cantidad de peces mientras minimiza el costo que implica cambiar de qué lado pescar. Además, evaluará la calidad de diferentes políticas de control (umbrales) para elegir acciones y cómo este puede ser diferente dependiendo de las probabilidades de del sistema.
Dinámica de estados: Dos posibles ubicaciones de los peces - derecha o izquierda. Sin embargo, aunque el cardúmen esté en un lado, puede que se hayan quedado unos cuantos peces en el lado opuesto (súcede cuando Prob_pescar_LO > 0).
Acciones: El agente puede decidir entre quedarse pescando en la ubicación actual o cambiar de ubicación.
Recompenzas: El agente obtiene una recompensa por cada pez capturado (si se está del mismo lado que el cardúmen entonces el agente capturará más peces con Prob_pesca_ML). Si está en el lado opuesto al cardúmen capturará peces con una probabilidad Prob_pesca_LO.
Costos: El agente paga un precio de por cada vez que cambia de lado- aunque para esta verisón del simulador todavía no se toma en cuenta el costo.
Maximización de la utilidad: Se seguirá una política que prescribe qué hacer en cada situación, la cual se determina por la ubicación del agente y su creencia posterior de la localización del cardúmen. La política es el umbral que indica cuándo se debe cambiar de lado para pescar.
En este simulador no se buscará la política óptima, pero al configurar ciertos parámetros con los sliders y sólo cambiando el umbral el lector puede tratar de averiguar cuál es el umbral adecuado dada esa configuración de parámetros.
Funciones.
Haz Clic para ver el Código
```{python}# Importsimport numpy as npfrom math import iscloseimport matplotlib.pyplot as plt# Figure Settingsimport ipywidgets as widgetsfrom IPython.display import HTML%config InlineBackend.figure_format ='retina'plt.style.use("https://raw.githubusercontent.com/NeuromatchAcademy/course-content/master/nma.mplstyle")# Plotting Functionsdef plot_fish(fish_state, ax=None):""" Plot the fish dynamics (states across time) """ T =len(fish_state) offset =3ifnot ax: fig, ax = plt.subplots(1, 1, figsize=(12, 3.5)) x = np.arange(0, T, 1) y = offset * (fish_state*2-1) ax.plot(y, color='cornflowerblue', markersize=10, linewidth=3.0, zorder=0) ax.fill_between(x, y, color='cornflowerblue', alpha=.3) ax.set_xlabel('Tiempo') ax.set_ylabel('Ubicación del cardúmen') ax.set_xlim([0, T]) ax.set_xticks([]) ax.xaxis.set_label_coords(1.05, .54) ax.set_ylim([-(offset+.5), offset+.5]) ax.set_yticks([-offset, offset]) ax.set_yticklabels(['Izda', 'Dcha']) ax.spines['bottom'].set_position('center')# plt.savefig('fish_loc.pdf')def plot_measurement(measurement, ax=None):""" Plot the measurements """ T =len(measurement) rel_pos =3 red_y = [] blue_y = []for idx, value inenumerate(measurement):if value ==0: blue_y.append([idx, -rel_pos])else: red_y.append([idx, rel_pos]) red_y = np.asarray(red_y) blue_y = np.asarray(blue_y)ifnot ax: fig, ax = plt.subplots(1, 1, figsize=(12, 3.5))iflen(red_y) >0: ax.plot(red_y[:, 0], red_y[:, 1], '*', markersize=8, color='crimson')iflen(blue_y) >0: ax.plot(blue_y[:, 0], blue_y[:, 1], '*', markersize=8, color='royalblue') ax.set_xlabel('Tiempo', fontsize=18) ax.set_ylabel('¿Pescó?') ax.set_xlim([0, T]) ax.set_xticks([]) ax.xaxis.set_label_coords(1.05, .54) ax.set_ylim([-rel_pos -.5, rel_pos +.5]) ax.set_yticks([-rel_pos, rel_pos]) ax.set_yticklabels(['No', '¡Sí!']) ax.spines['bottom'].set_position('center')# plt.savefig('measurement.pdf')def plot_act_loc(loc, act, ax_loc=None):""" Plot the action and location of T time points """ T =len(act)ifnot ax_loc: fig, ax_loc = plt.subplots(1, 1, figsize=(12, 2.5)) loc = loc*2-1 act_down = [] act_up = []for t inrange(1, T):if loc[t-1] ==-1and loc[t] ==1: act_up.append([t -0.5, 0])if loc[t-1] ==1and loc[t] ==-1: act_down.append([t -0.5, 0]) act_down = np.array(act_down) act_up = np.array(act_up) ax_loc.plot(loc, 'g.-', markersize=8, linewidth=5)iflen(act_down) >0: ax_loc.plot(act_down[:, 0], act_down[:, 1], 'rv', markersize=18, zorder=10, label='switch')iflen(act_up) >0: ax_loc.plot(act_up[:, 0], act_up[:, 1], 'r^', markersize=18, zorder=10) ax_loc.set_xlabel('Tiempo') ax_loc.set_ylabel('Tu estado') ax_loc.set_xlim([0, T]) ax_loc.set_xticks([]) ax_loc.xaxis.set_label_coords(1.05, .54)iflen(act_down) >0: ax_loc.legend(loc="upper right")eliflen(act_down) ==0andlen(act_up) >0: ax_loc.plot(act_up[:, 0], act_up[:, 1], 'r^', markersize=18, zorder=10, label='switch') ax_loc.legend(loc="upper right") ax_loc.set_ylim([-1.1, 1.1]) ax_loc.set_yticks([-1, 1]) ax_loc.tick_params(axis='both', which='major') ax_loc.set_yticklabels(['Izda', 'Dcha']) ax_loc.spines['bottom'].set_position('center')def plot_belief(belief, ax1=None, choose_policy=None):""" Plot the belief dynamics of T time points """ T = belief.shape[1]ifnot ax1: fig, ax1 = plt.subplots(1, 1, figsize=(12, 2.5)) ax1.plot(belief[1, :], color='midnightblue', markersize=10, linewidth=3.0) ax1.set_xlabel('Tiempo') ax1.set_ylabel('Creencia (dcha)') ax1.set_xlim([0, T]) ax1.set_xticks([]) ax1.xaxis.set_label_coords(1.05, 0.05) ax1.set_yticks([0, 1]) ax1.set_ylim([0, 1.1]) labels = [item.get_text() for item in ax1.get_yticklabels()] ax1.set_yticklabels([' 0', ' 1'])""" if choose_policy == "threshold": ax2 = ax1.twinx() ax2.plot(time_range, threshold * np.ones(time_range.shape), 'r--') ax2.plot(time_range, (1 - threshold) * np.ones(time_range.shape), 'c--') ax2.set_yticks([threshold, 1 - threshold]) ax2.set_ylim([0, 1.1]) ax2.tick_params(axis='both', which='major', labelsize=18) labels = [item.get_text() for item in ax2.get_yticklabels()] labels[0] = 'threshold to switch \n from left to right' labels[-1] = 'threshold to switch \n from right to left' ax2.set_yticklabels(labels) """def plot_dynamics(belief, loc, act, meas, fish_state, choose_policy):""" Plot the dynamics of T time points """if choose_policy =='threshold': fig, [ax0, ax_bel, ax_loc, ax1] = plt.subplots(4, 1, figsize=(12, 9)) plot_fish(fish_state, ax=ax0) plot_belief(belief, ax1=ax_bel) plot_measurement(meas, ax=ax1) plot_act_loc(loc, act, ax_loc=ax_loc)else: fig, [ax0, ax_bel, ax1] = plt.subplots(3, 1, figsize=(12, 7)) plot_fish(fish_state, ax=ax0) plot_belief(belief, ax1=ax_bel) plot_measurement(meas, ax=ax1) plt.tight_layout()#plt.savefig('active dynamics.pdf') plt.show()def belief_histogram(belief, bins=100):""" Plot the histogram of belief states """ fig, ax = plt.subplots(1, 1, figsize=(8, 6)) ax.hist(belief, bins) ax.set_xlabel('belief', fontsize=18) ax.set_ylabel('count', fontsize=18) plt.show()# Helper Functions# To generate a binomial with fixed "noise",# we generate a sequence of T numbers uniformly at randomT =100rnd_tele = np.random.uniform(0, 1, T)rnd_high_rwd = np.random.uniform(0, 1, T)rnd_low_rwd = np.random.uniform(0, 1, T)def get_randomness(T):global rnd_teleglobal rnd_high_rwdglobal rnd_low_rwd rnd_tele = np.random.uniform(0, 1, T) rnd_high_rwd = np.random.uniform(0, 1, T) rnd_low_rwd = np.random.uniform(0, 1, T)def binomial_tele(p):return np.array([1if p > rnd_tele[i] else0for i inrange(T)])def getRandomness(p, largeT):global rnd_teleglobal rnd_high_rwdglobal rnd_low_rwd rnd_tele = np.random.uniform(0, 1, largeT) rnd_high_rwd = np.random.uniform(0, 1, largeT) rnd_low_rwd = np.random.uniform(0, 1, largeT)return [np.array([1if p > rnd_tele[i] else0for i inrange(T)]), rnd_high_rwd, rnd_low_rwd]# def binomial_high_rwd(p):# return np.array([1 if p > rnd_high_rwd[i] else 0 for i in range(T)])# def binomial_low_rwd(p):# return np.array([1 if p > rnd_low_rwd[i] else 0 for i in range(T)])class ExcerciseError(AssertionError):passclass binaryHMM():def__init__(self, params, fish_initial=0, loc_initial=0):self.params = paramsself.fish_initial = fish_initialself.loc_initial = loc_initialdef fish_dynamics(self):""" fish state dynamics according to telegraph process Returns: fish_state (numpy array of int) """ p_stay, _, _, _ =self.params fish_state = np.zeros(T, int) # 0: left side and 1: right side# initialization fish_state[0] =self.fish_initial tele_operations = binomial_tele(p_stay) # 0: switch and 1: stayfor t inrange(1, T):# we use logical operation NOT XOR to determine the next state fish_state[t] =int(not(fish_state[t-1] ^ tele_operations[t]))return fish_statedef generate_process_lazy(self):""" fish dynamics and rewards if you always stay in the initial location without changing sides Returns: fish_state (numpy array of int): locations of the fish loc (numpy array of int): left or right site, 0 for left, and 1 for right rwd (numpy array of binary): whether a fish was catched or not """ _, p_low_rwd, p_high_rwd, _ =self.params fish_state =self.fish_dynamics() rwd = np.zeros(T, int) # 0: no food, 1: get foodfor t inrange(0, T):# new measurementif fish_state[t] !=self.loc_initial: rwd[t] =1if p_low_rwd > rnd_low_rwd[t] else0else: rwd[t] =1if p_high_rwd > rnd_high_rwd[t] else0# rwd[t] = binomial(1, p_rwd_vector[(fish_state[t] == loc[t]) * 1])return fish_state, self.loc_initial*np.ones(T), rwdclass binaryHMM_belief(binaryHMM):def__init__(self, params, fish_initial=0, loc_initial=1, choose_policy='threshold'): binaryHMM.__init__(self, params, fish_initial, loc_initial)self.choose_policy = choose_policydef generate_process(self):""" fish dynamics and measurements based on the chosen policy Returns: belief (numpy array of float): belief on the states of the two sites act (numpy array of string): actions over time loc (numpy array of int): left or right site measurement (numpy array of binary): whether a reward is obtained fish_state (numpy array of int): fish locations """ p_stay, low_rew_p, high_rew_p, threshold =self.params fish_state =self.fish_dynamics() # 0: left side; 1: right side loc = np.zeros(T, int) # 0: left side, 1: right side measurement = np.zeros(T, int) # 0: no food, 1: get food act = np.empty(T, dtype='object') # "stay", or "switch" belief = np.zeros((2, T), float) # the probability that the fish is on the left (1st element)# or on the right (2nd element),# the beliefs on the two boxes sum up to be 1 rew_prob = np.array([low_rew_p, high_rew_p])# initialization loc[0] =self.loc_initial measurement[0] =0 belief_0 = np.random.random(1)[0] belief[:, 0] = np.array([belief_0, 1- belief_0]) act[0] =self.policy(threshold, belief[:, 0], loc[0])for t inrange(1, T):if act[t -1] =="stay": loc[t] = loc[t -1]else: loc[t] =int(not(loc[t -1] ^0))# new measurement# measurement[t] = binomial(1, rew_prob[(fish_state[t] == loc[t]) * 1])if fish_state[t] != loc[t]: measurement[t] =1if low_rew_p > rnd_low_rwd[t] else0else: measurement[t] =1if high_rew_p > rnd_high_rwd[t] else0 belief[0, t] =self.belief_update(belief[0, t -1] , loc[t], measurement[t], p_stay, high_rew_p, low_rew_p) belief[1, t] =1- belief[0, t] act[t] =self.policy(threshold, belief[:, t], loc[t])return belief, loc, act, measurement, fish_statedef policy(self, threshold, belief, loc):""" chooses policy based on whether it is lazy policy or a threshold-based policy Args: threshold (float): the threshold of belief on the current site, when the belief is lower than the threshold, switch side belief (numpy array of float): the belief on the two sites loc (int) : the location of the agent Returns: act (string): "stay" or "switch" """ifself.choose_policy =="threshold": act = policy_threshold(threshold, belief, loc)ifself.choose_policy =="lazy": act = policy_lazy(belief, loc)return actdef belief_update(self, belief_past, loc, measurement, p_stay, high_rew_p, low_rew_p):""" using PAST belief on the LEFT box, CURRENT location and and measurement to update belief """ rew_prob_matrix = np.array([[1- high_rew_p, high_rew_p], [1- low_rew_p, low_rew_p]])# update belief posterior, p(s[t] | measurement(0-t), act(0-t-1)) belief_0 = (belief_past * p_stay + (1- belief_past) * (1- p_stay)) *\ rew_prob_matrix[(loc +1) //2, measurement] belief_1 = ((1- belief_past) * p_stay + belief_past * (1- p_stay)) *\ rew_prob_matrix[1-(loc +1) //2, measurement] belief_0 = belief_0 / (belief_0 + belief_1)return belief_0```
Widget.
Haz Clic para ver el Código
```{python}display(HTML('''<style>.widget-label { min-width: 15ex !important; }</style>'''))@widgets.interact(Umbral=widgets.FloatSlider(.2, description="Umbral", min=0., max=1., step=.01), Prob_quedarse=widgets.FloatSlider(.96, description="Prob_quedarse", min=0, max=1., step=.01), Prob_pesca_lado_opuesto=widgets.FloatSlider(.1, description="Prob_pesca_LO", min=0., max=1., step=.01), Prob_pesca_mismo_lado=widgets.FloatSlider(.3, description="Prob_pesca_ML", min=0., max=1., step=.01), Nueva_semilla=widgets.ToggleButtons(options=['Reusing', 'Refreshing'], description='Nueva semilla:', disabled=False, button_style='', # 'success', 'info', 'warning', 'danger' or '', icons=['check'] *2 ))def update_ex_4(Umbral, Nueva_semilla, Prob_quedarse, Prob_pesca_lado_opuesto, Prob_pesca_mismo_lado):""" p_stay: probability fish stay high_rew_p: p(catch fish) when you're on their side low_rew_p : p(catch fish) when you're on other side threshold: threshold of belief below which switching is taken """if Nueva_semilla =="Refreshing": get_randomness(T) params = [Prob_quedarse, Prob_pesca_mismo_lado, Prob_pesca_lado_opuesto, Umbral]#### initial condition for fish [fish_initial] and you [loc_initial] #### binaryHMM_test = binaryHMM_belief(params, fish_initial=0, loc_initial=0, choose_policy="threshold") belief, loc, act, measurement, fish_state = binaryHMM_test.generate_process() plot_dynamics(belief, loc, act, measurement, fish_state, binaryHMM_test.choose_policy)plt.show()```
El umbral en el anterior simulador funciona de la siguiente manera: Si la creencia (el cardúmen está a la derecha) es menor al umbral cambia de lado.
La probabilidades se midel de 0 a 1, siendo 0 nada probable y 1 100% probable.
Referencias
Farrell, S. y Lewandowsky, S. (2018). Bayesian Parameter Estimation. En Cambridge University Press (Ed.) Computational Modeling of Cognition Behavior (pp. 126-143). Cambridge University Press.