Cadenas de Markov absorbentes en Sage

1618 days ago by pang

Cadenas de Markov absorbentes en Sage

Presentado por Javier Sanz-Cruzado y Pablo Angulo en las III Jornadas Sage/Python celebradas en la Universidade de Vigo los días 21 y 22 de Junio de 2012.

 

4 Matemáticos quieren colaborar

  • 4 matemáticos en un departamento.
  • Cada uno trabaja en un área distinta (Geometría, Análisis, Topología, Estadística).
  • Deciden seguir el siguiente esquema:
    • Cada noche juegan a las cartas, y como resultado hay un ganador y un perdedor.
    • El perdedor se pasa al área del ganador.
    • Cuando los cuatro trabajen en el mismo área, escriben un artículo, y vuelven a empezar.

Ref: Markov Chains for Collaboration (Robert Mena y Will Murray), Mathematics Magazine

 

 

Ejemplo

(Fulano) Geometría, (Mengano) Análisis, (Zutano) Topología, (Renano) Estadística

Noche 1: Fulano gana, Renano pierde

(Fulano y Renano) Geometría, (Mengano) Análisis, (Zutano) Topología, Estadística

Noche 2: Mengano gana, Fulano pierde

(Fulano y Mengano) Análisis, (Renano) Geometría, (Zutano) Topología, Estadística

Noche 3: Mengano gana, Zutano pierde

(Fulano, Mengano y Zutano) Análisis, (Renano) Geometría, Topología, Estadística

Noche 4: Mengano gana, Fulano pierde

(Fulano, Mengano y Zutano) Análisis, (Renano) Geometría, Topología, Estadística

Noche 5: Zutano gana, Renano pierde

(Fulano, Mengano, Zutano y Renano) Análisis, Geometría, Topología, Estadística

 

 

Modelizando

Modelizamos el proceso con una cadena de Markov.

Nos fijamos sólo en cuántas áreas distintas quedan, y cuantas personas trabajan en cada área: Particiones de 4.

Ejemplo: En el estado (2,1,1) tenemos cuatro matemáticos trabajando en tres áreas: ((1,2), (3), (4))

  • Escogemos un ganador y un perdedor: $\frac{1}{4\cdot(4 -1)}$ posibilidades. Por ejemplo, gana 4 y pierde 2.
  • Sumamos 1 miembro al equipo del ganador y restamos 1 al del perdedor. Por ejemplo, pasamos a (1,1,2)
  • Ordenamos y eliminamos los equipos vacíos. En el ejemplo, pasamos a (2,1,1)

Noche 2: Mengano gana, Fulano pierde

def grupo(ls, j): 'Indice del grupo al que pertenece el matematico j-esimo' k = 0 for l,a in enumerate(ls): k+=a if k>=j: return l def pp(tupla): 'pretty print de una tupla' return ','.join(str(d) for d in tupla) def opciones(estadoi): N = sum(estadoi) p = 1/(N*(N-1)) opciones = {} #Si ya estan trabajando juntos, volvemos a empezar if len(estadoi) == 1: return {pp((1,)*N):1} #Elegimos un ganador y un perdedor for ganador, perdedor in ((a,b) for a in [1..N] for b in [1..N] if a != b): g1 = grupo(estadoi, ganador) g2 = grupo(estadoi, perdedor) #Sumamos 1 miembro al equipo del ganador y restamos 1 al del perdedor estadof = list(estadoi) estadof[g1] += 1 estadof[g2] -= 1 #Ordenamos y eliminamos los equipos vacíos estadof.sort(reverse=True) if estadof[-1] == 0: estadof = estadof[:-1] if pp(estadof) in opciones: opciones[pp(estadof)] += p else: opciones[pp(estadof)] = p return opciones def elegir_proyecto(N): estados = [tuple(ls) for ls in Partitions(N)] transiciones = dict((pp(estado), opciones(estado) ) for estado in estados) g = DiGraph(transiciones) return g N=4 g = elegir_proyecto(N) g.plot(edge_labels=True, layout='circular', vertex_size = 130*N).show(figsize = 8, ymin=-1.1) 
       

Distribución de probabilidad estable

En el campo de cadenas de Markov, suele ser interesante la distribución de probabilidad estable.

Pero, ¿qué significa en este caso?

Mejor aún, ¿qué preguntas queremos responder?

Quizá la pregunta más obvia es:

¿Cuánto tiempo tardan en ponerse de acuerdo sobre el área en la que van a trabajar?

Si hacemos el cálculo...

def matriz_probs(g): '''Matriz de probabilidades de transicion''' vs = g.vertices() M = matrix(QQ, len(vs)) vs = g.vertices() for v1,v2,p in g.edges(): j1 = vs.index(v1) j2 = vs.index(v2) M[j1,j2] = p return M 
       
N=4 #Grafo que representa la cadena g = elegir_proyecto(N) vs = g.vertices() M = matriz_probs(g) #Distribucion de probabilidad estable probs = (M - 1).left_kernel().basis()[0] probs /= sum(probs) #Dibujamos p = g.plot(edge_labels=True, layout='circular', vertex_size = 130*N, vertex_colors=dict(((1-p,1-p,1-0.001*random()),[v]) for v,p in zip(vs,probs))) print dict(zip(g.vertices(), probs)) p.show(axes = False, figsize = 8, ymin=-1.1) 
       
{'2,2': 1/5, '4': 1/10, '1,1,1,1': 1/10, '3,1': 2/5, '2,1,1': 1/5}
{'2,2': 1/5, '4': 1/10, '1,1,1,1': 1/10, '3,1': 2/5, '2,1,1': 1/5}

Tiempo medio que dura el proceso

  • Los caminos aleatorios que comienzan en (1,1,1,1) eventualmente llegan a (4).
  • Desde (4) pasan a (1,1,1,1) y se repite el proceso.
  • La probabilidad estable del vértice (4) es 1/10.
  • Esto siginifica que los caminos aleatorios pasan la décima parte del tiempo pasando de (4) a (1,1,1,1) (una sóla noche)
  • El resto del tiempo, lo pasan viajando del (1,1,1,1) al (4).
  • Luego el viaje del (1,1,1,1) al (4) dura en promedio 9 noches.

Pero con este truco sólo sacamos el tiempo medio. Ésta no es la forma adecuada de estudiar el problema.



Cadenas de Markov absorbentes

  • Uno (o varios) estados terminales.
  • Si sólo hay un estado terminal, la probabilidad estable es trivial!

Modificamos el ejemplo anterior: al llegar al estado (4), nos quedamos allí indefinidamente.

def opcionesA(estadoi): N = sum(estadoi) p = 1/(N*(N-1)) opciones = {} #Si ya estan trabajando juntos, no hacemos nada especial # if len(estadoi) == 1: # return {pp((1,)*N):1} #Elegimos un ganador y un perdedor for ganador, perdedor in ((a,b) for a in [1..N] for b in [1..N] if a != b): g1 = grupo(estadoi, ganador) g2 = grupo(estadoi, perdedor) #Sumamos 1 miembro al equipo del ganador y restamos 1 al del perdedor estadof = list(estadoi) estadof[g1] += 1 estadof[g2] -= 1 #Ordenamos y eliminamos los equipos vacíos estadof.sort(reverse=True) if estadof[-1] == 0: estadof = estadof[:-1] if pp(estadof) in opciones: opciones[pp(estadof)] += p else: opciones[pp(estadof)] = p return opciones def elegir_proyectoA(N): estados = [tuple(ls) for ls in Partitions(N)] transiciones = dict((pp(estado), opcionesA(estado) ) for estado in estados) g = DiGraph(transiciones) return g 
       
N=4 #Grafo que representa la cadena g = elegir_proyectoA(N) vs = g.vertices() M = matriz_probs(g) #Distribucion de probabilidad estable probs = (M - 1).left_kernel().basis()[0] probs /= sum(probs) #Dibujamos p = g.plot(edge_labels=True, layout='circular', vertex_size = 130*N, vertex_colors=dict(((1-p,1-p,1-0.001*random()),[v]) for v,p in zip(vs,probs))) print dict(zip(g.vertices(), probs)) p.show(axes = False, figsize = 8, ymin=-1.1) 
       
{'2,2': 0, '4': 1, '1,1,1,1': 0, '3,1': 0, '2,1,1': 0}
{'2,2': 0, '4': 1, '1,1,1,1': 0, '3,1': 0, '2,1,1': 0}

Ahora podemos calcular la probabilidad de que el proceso termine antes de tiempo k:

Es la entrada en la primera fila y última columna de la matriz M^k!

N=4 g = elegir_proyectoA(N) vs = g.vertices() M = matriz_probs(g) n3 = lambda x: x.n(digits = 3) for k in [0..8]: probs = vector([1,0,0,0,0])*M^k probs /= sum(probs) p = g.plot(edge_labels=True, layout='circular', vertex_size = 400, vertex_colors=dict(((1-p,1-p,1-0.001*random()),[v]) for v,p in zip(vs,probs))) print k, dict(zip(g.vertices(), map(n3,probs))) p.show(axes = False, figsize = 8, ymin=-1.1) 
       
0 {'2,2': 0.000, '4': 0.000, '1,1,1,1': 1.00, '3,1': 0.000, '2,1,1':
0.000}

1 {'2,2': 0.000, '4': 0.000, '1,1,1,1': 0.000, '3,1': 0.000, '2,1,1':
1.00}

2 {'2,2': 0.167, '4': 0.000, '1,1,1,1': 0.000, '3,1': 0.333, '2,1,1':
0.500}

3 {'2,2': 0.222, '4': 0.0833, '1,1,1,1': 0.000, '3,1': 0.444, '2,1,1':
0.250}

4 {'2,2': 0.227, '4': 0.194, '1,1,1,1': 0.000, '3,1': 0.454, '2,1,1':
0.125}

5 {'2,2': 0.210, '4': 0.308, '1,1,1,1': 0.000, '3,1': 0.420, '2,1,1':
0.0625}

6 {'2,2': 0.185, '4': 0.413, '1,1,1,1': 0.000, '3,1': 0.371, '2,1,1':
0.0312}

7 {'2,2': 0.160, '4': 0.505, '1,1,1,1': 0.000, '3,1': 0.319, '2,1,1':
0.0156}

8 {'2,2': 0.136, '4': 0.585, '1,1,1,1': 0.000, '3,1': 0.271, '2,1,1':
0.00781}
0 {'2,2': 0.000, '4': 0.000, '1,1,1,1': 1.00, '3,1': 0.000, '2,1,1': 0.000}

1 {'2,2': 0.000, '4': 0.000, '1,1,1,1': 0.000, '3,1': 0.000, '2,1,1': 1.00}

2 {'2,2': 0.167, '4': 0.000, '1,1,1,1': 0.000, '3,1': 0.333, '2,1,1': 0.500}

3 {'2,2': 0.222, '4': 0.0833, '1,1,1,1': 0.000, '3,1': 0.444, '2,1,1': 0.250}

4 {'2,2': 0.227, '4': 0.194, '1,1,1,1': 0.000, '3,1': 0.454, '2,1,1': 0.125}

5 {'2,2': 0.210, '4': 0.308, '1,1,1,1': 0.000, '3,1': 0.420, '2,1,1': 0.0625}

6 {'2,2': 0.185, '4': 0.413, '1,1,1,1': 0.000, '3,1': 0.371, '2,1,1': 0.0312}

7 {'2,2': 0.160, '4': 0.505, '1,1,1,1': 0.000, '3,1': 0.319, '2,1,1': 0.0156}

8 {'2,2': 0.136, '4': 0.585, '1,1,1,1': 0.000, '3,1': 0.271, '2,1,1': 0.00781}
 
       

Representamos la convergencia en forma de animaciones:

  • Representación en forma de grafo
def animateMarkovChain(g, prob_inicial, longitud = 10): vs = g.vertices() probs = vector(prob_inicial[v] for v in range(len(vs))) def plotg(probs, j): return (g.plot(edge_labels=True, layout='circular', vertex_size = 400, vertex_colors=dict(((1-p,1-p,1-0.001*random()),[v]) for v,p in zip(vs,probs)))+ text('%d: %s'%(j,map(n3,probs)), (0,-1.2), fontsize=12)) ps = [plotg(probs,0), plotg(probs,0)] M = matriz_probs(g) for j in xrange(longitud): probs = probs*M ps.append(plotg(probs, j)) return animate(ps, axes = False, figsize = 10) 
       
show(animateMarkovChain(g, [1,0,0,0,0]), delay = 50) 
       

Otra representación interesante usando Matrix Plot

animate([(M^i).plot() for i in [1..20]]).show(delay=30) 
       

Distribución del tiempo de llegada

... usamos una variable simbólica k

Problema: No se puede elevar una matriz a un exponente directamente

var('k') M^k 
       
Traceback (click to the left of this block for traceback)
...
NotImplementedError: non-integral exponents not supported
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_13.py", line 10, in <module>
    exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("dmFyKCdrJykKTV5r"),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))
  File "", line 1, in <module>
    
  File "/tmp/tmpZlr2xh/___code___.py", line 3, in <module>
    exec compile(u'M**k
  File "", line 1, in <module>
    
  File "matrix0.pyx", line 4787, in sage.matrix.matrix0.Matrix.__pow__ (sage/matrix/matrix0.c:26742)
  File "element.pyx", line 1790, in sage.structure.element.RingElement.__pow__ (sage/structure/element.c:14901)
  File "element.pyx", line 3363, in sage.structure.element.generic_power_c (sage/structure/element.c:26025)
NotImplementedError: non-integral exponents not supported

Así que usamos la ténica aprendida en Álgebra Lineal:

$$M^k=P\cdot D^k\cdot P^{-1}$$

donde $M=P\cdot D\cdot P^{-1}$ es la descomposición de Jordan.

Por el momento asumimos que la matriz es diagonalizable.

def pot(D,k): return diagonal_matrix([l^k for l in D.diagonal()]) D,P=M.jordan_form(transformation=True) show(P*pot(D,k)*~P) 
       

La prob que nos interesa es ...

pcum = (P*pot(D,k)*~P)[0,-1] p = pcum(k=k) - pcum(k=k-1) p 
       
(1/2)^k - (1/2)^(k - 1) - 9/5*(5/6)^k + 9/5*(5/6)^(k - 1)
(1/2)^k - (1/2)^(k - 1) - 9/5*(5/6)^k + 9/5*(5/6)^(k - 1)

Ojo! pcum(k=0) es incorrecto
luego p(k=1)   es incorrecto

Sabemos que p(k=1)=0!

point2d([(k0, p(k=k0)) for k0 in [1..20]]) 
       
  • La fórmula que hemos usado no funciona cuando k<=0
  • La matriz M es singular!
  • La técnica anterior funciona para k<=0 si y sólo si la matriz es no singular.
for k0 in [1,0,-1]: print k0,':' show(M^k0) show( (P*pot(D,k)*~P)(k=k0) ) 
       
1 :


0 :


-1 :
Traceback (click to the left of this block for traceback)
...
ZeroDivisionError: input matrix must be nonsingular
1 :


0 :


-1 :
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_17.py", line 10, in <module>
    exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("Zm9yIGswIGluIFsxLDAsLTFdOgogICAgcHJpbnQgazAsJzonCiAgICBzaG93KE1eazApCiAgICBzaG93KCAoUCpwb3QoRCxrKSp+UCkoaz1rMCkgKQ=="),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))
  File "", line 1, in <module>
    
  File "/tmp/tmpYu2_Pr/___code___.py", line 3, in <module>
    exec compile(u"for k0 in [_sage_const_1 ,_sage_const_0 ,-_sage_const_1 ]:\n    print k0,':'\n    show(M**k0)\n    show( (P*pot(D,k)*~P)(k=k0) )" + '\n', '', 'single')
  File "", line 3, in <module>
    
  File "matrix0.pyx", line 4787, in sage.matrix.matrix0.Matrix.__pow__ (sage/matrix/matrix0.c:26742)
  File "element.pyx", line 1790, in sage.structure.element.RingElement.__pow__ (sage/structure/element.c:14901)
  File "element.pyx", line 3385, in sage.structure.element.generic_power_c (sage/structure/element.c:26349)
  File "matrix_rational_dense.pyx", line 650, in sage.matrix.matrix_rational_dense.Matrix_rational_dense.__invert__ (sage/matrix/matrix_rational_dense.c:9350)
  File "matrix_rational_dense.pyx", line 756, in sage.matrix.matrix_rational_dense.Matrix_rational_dense.__invert__main (sage/matrix/matrix_rational_dense.c:9945)
  File "matrix_rational_dense.pyx", line 2658, in sage.matrix.matrix_rational_dense.Matrix_rational_dense._invert_pari (sage/matrix/matrix_rational_dense.c:24393)
ZeroDivisionError: input matrix must be nonsingular

Afortunadamente, la solución es sencilla: no necesitamos trabajar con k0=0.

point2d([(k0, p(k=k0)) for k0 in [2..20]]) 
       

Ahora podríamos calcular la media, varianza...

pero Sage dice que no sabe!

sum(p,k,2,oo) 
       
-1/25*sum(-(9*2^k*5^k - 25*6^k)*2^(-k)*6^(-k), k, 2, +Infinity)
-1/25*sum(-(9*2^k*5^k - 25*6^k)*2^(-k)*6^(-k), k, 2, +Infinity)

Sin embargo, sabe sumar las partes de p:

print -1/25*sum(-(9*2^k*5^k )*2^(-k)*6^(-k), k, 2, +Infinity) print -1/25*sum(-( - 25*6^k)*2^(-k)*6^(-k), k, 2, +Infinity) 
       
3/2
-1/2
3/2
-1/2

Así que le ayudamos, expandiendo la expresión simbólica, y luego sumando cada sumando por separado.

def mi_suma(p, v, li, ls): return sum(sum(o, v, li, ls ) for o in p.expand().operands() ) print 'Probabilidad total: ', mi_suma(p,k,2,oo) m = mi_suma(p*k,k,2,oo) print 'Media: ', m print 'Varianza: ', mi_suma(p*(k-m)^2,k,2,oo) 
       
Probabilidad total:  1
Media:  9
Varianza:  32
Probabilidad total:  1
Media:  9
Varianza:  32

Cálculo directo del tiempo medio

Buscamos ahora otro método para calcular el tiempo medio.

       $EX_i$ $\rightarrow$ Tiempo medio en alcanzar el estado absorbente desde el estado $i$

Queremos resolver un sistema de ecuaciones lineales, con ecuaciones del tipo 

$$EX_i = 1 + p(\text{Ir al nodo a}) · EX_a + p(\text{Ir al nodo b}) · EX_b ....$$

para un estado $i$ que no es absorbente, y la ecuación:

$$EX_i = 0$$

para un estado $i$ es absorbente.

 

El sistema en forma matricial

Sea $E$ el vector:

$$(EX_1, \dots, EX_N)$$

$E$ resuelve el sistema:

$$E\cdot A = v$$

donde

$A \rightarrow$ Matriz cuyas filas son:

  • Si el vértice i-ésimo no es absorbente
    • $1 - M[i][i]$ en la posición i-ésima ($M \rightarrow $Matriz probabilidades de transición)
    • $- M[i][j]$ en el resto
  • Si el vértice i-ésimo es absorbente
    • 1 en la posición i
    • 0 en el resto

$v \rightarrow$ Vector de términos independientes cuyos elementos son:

  • 1 si el vértice i-ésimo no es absorbente
  • 0 si el vértice i-ésimo es absorbente
def matriz_fundamental(g): M = matriz_probs(g) rows = [] for k,row in enumerate(M.rows()): if M[k,k]==1: #estado absorbente rows.append(row) else: row = -row row[k] +=1 rows.append(row) return matrix(QQ,rows) def vector_fundamental(g): M = matriz_probs(g) L = M.ncols() row = [0 if M[k,k]==1 else 1 for k in range(L)] return vector(QQ,row) show(matriz_fundamental(g)) show(vector_fundamental(g)) print 'Tiempo medio en llegar al estado absorbente' g = elegir_proyectoA(N) for v,t in zip(g.vertices() , matriz_fundamental(g)\vector_fundamental(g)): print v,':',t 
       


Tiempo medio en llegar al estado absorbente
1,1,1,1 : 9
2,1,1 : 8
2,2 : 7
3,1 : 11/2
4 : 0
Tiempo medio en llegar al estado absorbente
1,1,1,1 : 9
2,1,1 : 8
2,2 : 7
3,1 : 11/2
4 : 0

El tiempo medio partiendo del estado (1,1,1,1) se corresponde con el calculado anteriormente.




Araña en un cubo

  • Tenemos una araña en un vértice de un cubo
  • Tenemos veneno en dos vértices del cubo
  • La araña se mueve siguiendo las aristas del cubo, de un vértice a otro
  • Si llega a un vértice con veneno, muere.

Buscamos saber cual es la probabilidad de que muera en un vértice u en otro según el vértice del que parta la araña

Ejemplo

Tenemos que los venenos se hallan en los vértices 2 y 5

Tenemos una cadena de Markov absorbente con varios estados terminales

def grafoSpider(posVen1, posVen2): '''Convierte el grafo del cubo en el grafo que representa la cadena de Markov en funcion de la posicion de los venenos''' g = graphs.CubeGraph(3) g.allow_loops(True) grafoDir = g.adjacency_matrix() g1 = DiGraph(grafoDir, loops=True) for ed in g1.edges(): g1.set_edge_label(ed[0],ed[1], 1/3) for i in range(len(g1.vertices())): if g1.has_edge(posVen1, i): g1.delete_edge(posVen1, i) if g1.has_edge(posVen2, i): g1.delete_edge(posVen2, i) g1.add_edge(posVen1,posVen1, 1) g1.add_edge(posVen2, posVen2, 1) return g1 grafoSpider(2,5).plot(edge_labels=True, layout='spring', vertex_size = 300).show(figsize = 6) 
       

Ejemplo partiendo la araña del vértice 0

animateMarkovChain(grafoSpider(2,5), [1,0,0,0,0,0,0,0]).show(delay = 50) 
       
M1 = matrix([row/sum(row) for row in grafoSpider(2,5).adjacency_matrix().rows()]) animate([(M1^i).plot() for i in [1..20]]).show(delay=30) 
       

En este caso, la distribución de probabilidad estable no es trivial: Se reparte entre los diferentes estados absorbentes.

Vamos a calcular estas probabilidades para cada vértice inicial:

  • $p_X(v_Y) \rightarrow$ Probabilidad de llegar al veneno X desde el vértice Y

Las ecuaciones ahora serán del tipo: $$p_X(v_0) = p_X(v_0) · p(\text{Llegar  a }v_0) + p_X(v_1) · p(\text{Llegar  a }v_1) + ...$$

para un estado no absorbente y para los estados absorbentes sabemos que:

  • $p_{Ven1}(v_{Ven1}) = p_{Ven2}(v_{Ven2}) = 1$
  • $p_{Ven1}(v_{Ven2}) = p_{Ven2}(v_{Ven1}) = 0$

El sistema en forma matricial

Sea $P$ el vector:

$$(p_X(v_1), \dots, p_X(v_N))$$

$P$ resuelve el sistema:

$$P\cdot A = v$$

donde

$A \rightarrow$ Matriz cuyas filas son:

  • Si el vértice i-ésimo no es absorbente
    • $1 - M[i][i]$ en la posición i-ésima ($M \rightarrow $Matriz de transición)
    • $- M[i][j]$ en el resto
  • Si el vértice i-ésimo es absorbente
    • 1 en la posición i
    • 0 en el resto

$v \rightarrow$ Vector de términos independientes cuyos elementos son:

  • 1 si el vértice i-ésimo es el vértice absorbente al que queremos llegar
  • 0 para el resto

Por ejemplo, calculamos la probabilidad de morir en el veneno 1:

print 'PROBABILIDAD DE CAER EN VENENO 1' posVen1 = 2 posVen2 = 5 g = grafoSpider(posVen1, posVen2) print 'Matriz del sistema' show(matriz_fundamental(g)) print 'Terminos independientes del sistema' probSpider = [0]*g.num_verts() probSpider[posVen1] = 1 vectorProbSpider = vector(probSpider) show(vectorProbSpider) print 'Solucion del sistema' for v,t in zip(g.vertices() , matriz_fundamental(g)\vectorProbSpider): print v,':',t 
       
PROBABILIDAD DE CAER EN VENENO 1
Matriz del sistema

Terminos independientes del sistema

Solucion del sistema
0 : 3/5
1 : 2/5
2 : 1
3 : 3/5
4 : 2/5
5 : 0
6 : 3/5
7 : 2/5
PROBABILIDAD DE CAER EN VENENO 1
Matriz del sistema

Terminos independientes del sistema

Solucion del sistema
0 : 3/5
1 : 2/5
2 : 1
3 : 3/5
4 : 2/5
5 : 0
6 : 3/5
7 : 2/5

Para el veneno 2, tenemos que:

$$p_{Ven2}(v_i) = 1 - p_{Ven1}(v_i)$$

Tiempo medio en morir

Calculamos con la técnica anterior el tiempo medio en morir (es decir, en llegar a uno de los dos estados absorbentes).

print 'TIEMPO MEDIO EN MORIR' print 'Matriz del sistema' show(matriz_fundamental(g)) print 'Terminos independientes del sistema' show(vector_fundamental(g)) print 'Solucion del sistema' for v,t in zip(g.vertices() , matriz_fundamental(g)\vector_fundamental(g)): print v,':',t 
       
TIEMPO MEDIO EN MORIR
Matriz del sistema

Terminos independientes del sistema

Solucion del sistema
0 : 3
1 : 3
2 : 0
3 : 3
4 : 3
5 : 0
6 : 3
7 : 3
TIEMPO MEDIO EN MORIR
Matriz del sistema

Terminos independientes del sistema

Solucion del sistema
0 : 3
1 : 3
2 : 0
3 : 3
4 : 3
5 : 0
6 : 3
7 : 3

La misma técnica funciona para cualquier grafo absorbente.

#Descomenta una de las tres lineas siguientes #g = DiGraph(graphs.CirculantGraph(7,4), loops=True) #g = DiGraph(graphs.DesarguesGraph(), loops=True) g = DiGraph(graphs.PetersenGraph(), loops=True) #Añadimos dos venenos al azar pos1, pos2 = Combinations(g.vertices(), 2).random_element() g.add_edges([(pos1,'V1'),('V1','V1'),(pos2,'V2'), ('V2','V2')]) d = g.get_pos() d['V1'] = (-1.5,0) d['V2'] = (1.5,0) g.set_pos(d) for v1 in g.vertices(): d = g.out_degree(v1) for v2 in g.neighbors_out(v1): g.set_edge_label(v1, v2, 1/d) g.plot(edge_labels=True, vertex_size = 300).show(figsize = 8) 
       
print 'PROBABILIDAD DE CAER EN VENENO 2' posVen2 = g.num_verts() - 1 probPetersen = [0]*g.num_verts() probPetersen[posVen2] = 1 vectorprobPetersen = vector(probPetersen) print zip(g.vertices() , matriz_fundamental(g)\vectorprobPetersen) probs = matriz_fundamental(g)\vectorprobPetersen g.plot(vertex_colors=dict( ( (1-p,1-p,1-0.001*random()),[v]) for v,p in zip(g.vertices(),probs) ), vertex_size = 300 ).show(figsize = 8) 
       
PROBABILIDAD DE CAER EN VENENO 2
[(0, 1/2), (1, 1/2), (2, 1/2), (3, 13/28), (4, 5/14), (5, 9/14), (6,
1/2), (7, 15/28), (8, 15/28), (9, 13/28), ('V1', 0), ('V2', 1)]
PROBABILIDAD DE CAER EN VENENO 2
[(0, 1/2), (1, 1/2), (2, 1/2), (3, 13/28), (4, 5/14), (5, 9/14), (6, 1/2), (7, 15/28), (8, 15/28), (9, 13/28), ('V1', 0), ('V2', 1)]
 
       

Tiempo de llegada infinito

¿Qué ocurriría si la probabilidad de no llegar a alguno de los venenos fuese >0?

En este caso, el grafo no es absorbente.

Hemos de estudiar:

  • Probabilidad de llegar a cada veneno
  • Tiempo medio

PROBABILIDAD DE LLEGAR A CADA VENENO

grafo = grafoSpider(2,5) grafo.add_vertices([8,9,10,11]) grafo.add_edges([(8,9,1), (9,10,1), (10,8,1), (11,8,1/2), (11,6,1/2)]) show(grafo, layout = 'circular', figsize = 8, edge_labels = True) 
       
M2 = matrix([row/sum(row) for row in grafo.adjacency_matrix().rows()]) animate([(M2^i).plot() for i in [1..20]]).show(delay=30) 
       

 

Observamos que:

  • Si partimos de cualquier nodo del 1 al 7, la probabilidad de caer en el veneno se estabiliza 
  • Si partimos de los nodos 8, 9 o 10, la probabilidad no se estabiliza
    • p(Llegar a algun veneno) = 0
  • Si partimos del vértice 11, la probabilidad en los nodos con veneno se estabiliza, mientras que no ocurre lo mismo en el bucle.
    • $p_1(v_{11}) + p_2(v_{11}) < 1$
posVen1 = 2 posVen2 = 5 
       
print 'PROBABILIDAD DE CAER EN EL VENENO 1' print 'Matriz del sistema' show(matriz_fundamental(grafo)) probSpider = [0]*grafo.num_verts() probSpider[posVen1] = 1 vectorProbSpider = vector(probSpider) print 'Terminos independientes del sistema' show(vectorProbSpider) print 'Solucion del sistema' for v,t in zip(grafo.vertices() , matriz_fundamental(grafo)\vectorProbSpider): print v,':',t 
       
PROBABILIDAD DE CAER EN EL VENENO 1
Matriz del sistema

Terminos independientes del sistema

Solucion del sistema
0 : 3/5
1 : 2/5
2 : 1
3 : 3/5
4 : 2/5
5 : 0
6 : 3/5
7 : 2/5
8 : -3/5
9 : -3/5
10 : -3/5
11 : 0
PROBABILIDAD DE CAER EN EL VENENO 1
Matriz del sistema

Terminos independientes del sistema

Solucion del sistema
0 : 3/5
1 : 2/5
2 : 1
3 : 3/5
4 : 2/5
5 : 0
6 : 3/5
7 : 2/5
8 : -3/5
9 : -3/5
10 : -3/5
11 : 0

¡Nos encontramos un problema!

Las ecuaciones

$$p_{Ven1}(8) = p_{Ven1}(9)  $$

$$p_{Ven1}(9) = p_{Ven1}(10)  $$

$$p_{Ven1}(10) = p_{Ven1}(8)  $$

 no tienen solución

print 'Determinante: ',det(matriz_fundamental(grafo)) matriz_fundamental(grafo).right_kernel() 
       
Determinante:  0
Vector space of degree 12 and dimension 1 over Rational Field
Basis matrix:
[  0   0   0   0   0   0   0   0   1   1   1 1/2]
Determinante:  0
Vector space of degree 12 and dimension 1 over Rational Field
Basis matrix:
[  0   0   0   0   0   0   0   0   1   1   1 1/2]

Sin embargo, sabemos que en estos vértices la probabilidad de llegar al veneno es 0.

Hemos de modificar el sistema. ¿Cómo detectamos estos vértices?

Idea: Distancias entre nodos del grafo

print grafo.distance(1,2) print grafo.distance(8,2) 
       
2
+Infinity
2
+Infinity

Tenemos que:

  • $d(v_i, v_{VenenoX}) = \infty \Leftrightarrow p_X(v_i) = 0$

Por  tanto, modificamos el sistema para tratar con estos problemas:

Para el estado absorbente objetivo:

$$p_X(X)=1$$

Para otro estado absorbente $v_i$:

$$p_X(v_i)=0$$

Si $d(v_i, v_{VenenoX})=\infty$:

$$p_X(v_i)=0$$

Si $d(v_i, v_{VenenoX})<\infty$:

$$p_X(v_i) = p_X(v_0) · p(\text{Llegar  a }v_0) + p_X(v_1) · p(\text{Llegar  a }v_1) + ...$$

def generaVectoresSpiderGeneral(pos, posVen1, posVen2, h): ''' genera vectores para situarlos en una matriz para resolver un sistema de ecuaciones lineales, dado un grafo dirigido''' counter = 0 ls=[] mat = h.adjacency_matrix() #Calculamos las posibles elecciones que tiene la arana desde el vertice pos e inicializamos la lista que se #convertira en el vector para la matriz A, tal que Ax=b, donde b es el vector solucion del sistema for i in range(len(h.vertices())): if h.has_edge(pos, i) == true: #Contamos el numero de aristas entre el vertice inicial y el destino counter += mat[pos][i] ls.append(0) ls[pos] = 1 if(counter != 0 and pos != posVen1 and pos != posVen2): #Si podemos llegar hasta el veneno if h.distance(pos, posVen1) != +oo or h.distance(pos, posVen2) != +oo: for j in range(len(h.vertices())): if h.has_edge(pos, j) == true: ls[j] = ls[j] - mat[pos][j]/counter v = vector(ls) return v def generaMatrizSpiderGeneral(posVen1, posVen2, h): '''Genera la matriz para resolver las probabilidades de que la arana caiga en el veneno''' mat = [] for i in range(len(h.vertices())): mat.append(generaVectoresSpiderGeneral(i, posVen1, posVen2, h)) m = matrix(mat) return m 
       
def printProbabilidadVeneno(posVen1, posVen2, h): '''Imprime las probabilidades''' print 'VENENO 1' print 'Matriz del sistema' show(generaMatrizSpiderGeneral(posVen1, posVen2, h)) print 'Terminos independientes' probSpider = [0]*grafo.num_verts() probSpider[posVen1] = 1 vectorProbSpider = vector(probSpider) show(vectorProbSpider) value = generaMatrizSpiderGeneral(posVen1, posVen2, h)\vectorProbSpider print 'Solucion del sistema' for v,t in zip(h.vertices() , value): print v,':',t print '' print 'VENENO 2' print 'Matriz del sistema' show(generaMatrizSpiderGeneral(posVen1, posVen2, h)) print 'Terminos independientes' probSpider = [0]*grafo.num_verts() probSpider[posVen2] = 1 vectorProbSpider = vector(probSpider) show(vectorProbSpider) value2 = generaMatrizSpiderGeneral(posVen1, posVen2, h)\vectorProbSpider print 'Solucion del sistema' for v,t in zip(h.vertices() , value2): print v,':',t print '' print 'SUMA' for i in range(len(value)): print 'p1(v%d) + p2(v%d) = %f'%(i,i,value[i] + value2[i]) 
       
printProbabilidadVeneno(2,5,grafo) 
       
VENENO 1
Matriz del sistema

Terminos independientes

Solucion del sistema
0 : 3/5
1 : 2/5
2 : 1
3 : 3/5
4 : 2/5
5 : 0
6 : 3/5
7 : 2/5
8 : 0
9 : 0
10 : 0
11 : 3/10

VENENO 2
Matriz del sistema

Terminos independientes

Solucion del sistema
0 : 2/5
1 : 3/5
2 : 0
3 : 2/5
4 : 3/5
5 : 1
6 : 2/5
7 : 3/5
8 : 0
9 : 0
10 : 0
11 : 1/5

SUMA
p1(v0) + p2(v0) = 1.000000
p1(v1) + p2(v1) = 1.000000
p1(v2) + p2(v2) = 1.000000
p1(v3) + p2(v3) = 1.000000
p1(v4) + p2(v4) = 1.000000
p1(v5) + p2(v5) = 1.000000
p1(v6) + p2(v6) = 1.000000
p1(v7) + p2(v7) = 1.000000
p1(v8) + p2(v8) = 0.000000
p1(v9) + p2(v9) = 0.000000
p1(v10) + p2(v10) = 0.000000
p1(v11) + p2(v11) = 0.500000
VENENO 1
Matriz del sistema

Terminos independientes

Solucion del sistema
0 : 3/5
1 : 2/5
2 : 1
3 : 3/5
4 : 2/5
5 : 0
6 : 3/5
7 : 2/5
8 : 0
9 : 0
10 : 0
11 : 3/10

VENENO 2
Matriz del sistema

Terminos independientes

Solucion del sistema
0 : 2/5
1 : 3/5
2 : 0
3 : 2/5
4 : 3/5
5 : 1
6 : 2/5
7 : 3/5
8 : 0
9 : 0
10 : 0
11 : 1/5

SUMA
p1(v0) + p2(v0) = 1.000000
p1(v1) + p2(v1) = 1.000000
p1(v2) + p2(v2) = 1.000000
p1(v3) + p2(v3) = 1.000000
p1(v4) + p2(v4) = 1.000000
p1(v5) + p2(v5) = 1.000000
p1(v6) + p2(v6) = 1.000000
p1(v7) + p2(v7) = 1.000000
p1(v8) + p2(v8) = 0.000000
p1(v9) + p2(v9) = 0.000000
p1(v10) + p2(v10) = 0.000000
p1(v11) + p2(v11) = 0.500000

TIEMPO MEDIO

Empleando lo anterior, tenemos que:

  • $d(v_i, v_{VenenoX}) = \infty \Rightarrow p_X(v_i) = 0 \Rightarrow EX_i = \infty$

Además:

  • Si $d(v_i, v_j) \ne \infty $ donde $d(v_j, v_{VenenoX}) = \infty$, entonces $EX_i  = \infty$

Calculamos estos vértices:

def safeVertex(posVen1, posVen2, h): '''Halla los vertices que no estan conectados a los vertices con el veneno y devuelve una lista con ellos''' ls = [] for i in h.vertices(): if(i != posVen1 and i!= posVen2): if(h.distance(i, posVen1) == +oo and h.distance(i, posVen2) == +oo): ls.append(i) return ls def infinityAverage(posVen1, posVen2, g): '''Calcula los vertices que tienen media infinita en un grafo''' seguros = safeVertex(posVen1, posVen2, g) infinitos = [] vert = g.vertices() for i in range (len(vert)): for j in range (len(seguros)): if(vert[i] != posVen1 and vert[i] != posVen2): if g.distance(vert[i],seguros[j]) != +oo: infinitos.append(i) break return infinitos 
       
print 'Vertices seguros:',safeVertex(2,5,grafo) print 'Vertices con media infinita:',infinityAverage(2,5,grafo) 
       
Vertices seguros: [8, 9, 10]
Vertices con media infinita: [8, 9, 10, 11]
Vertices seguros: [8, 9, 10]
Vertices con media infinita: [8, 9, 10, 11]

Con esta función, y el sistema de ecuaciones empleado anteriormente, calcularemos de nuevo el tiempo medio

def mediaVenenoGeneral(posVen1, posVen2, h): '''Calcula la media de tiempo que tarda en llegar para cualquier grafo dirigido''' ls = [] res = [] resultado = generaMatrizSpiderGeneral(posVen1, posVen2, h)\vector_fundamental(h) listaInf = infinityAverage(posVen1, posVen2, h) for i in range(len(h.vertices())): if i in listaInf: res.append(+oo) else: res.append(resultado[i]) return res 
       
print 'TIEMPO MEDIO' print 'Matriz del sistema' show(generaMatrizSpiderGeneral(2, 5, grafo)) print 'Terminos independientes' show(vector_fundamental(grafo)) value = mediaVenenoGeneral(2,5,grafo) print 'Solucion del sistema' for v,t in zip(grafo.vertices() , value): print v,':',t 
       
TIEMPO MEDIO
Matriz del sistema

Terminos independientes

Solucion del sistema
0 : 3
1 : 3
2 : 0
3 : 3
4 : 3
5 : 0
6 : 3
7 : 3
8 : +Infinity
9 : +Infinity
10 : +Infinity
11 : +Infinity
TIEMPO MEDIO
Matriz del sistema

Terminos independientes

Solucion del sistema
0 : 3
1 : 3
2 : 0
3 : 3
4 : 3
5 : 0
6 : 3
7 : 3
8 : +Infinity
9 : +Infinity
10 : +Infinity
11 : +Infinity
 
       
 
       

Lanzamiento de monedas

  • Tenemos una moneda no lastrada
  • Lanzamos la moneda hasta obtener cuatro caras

¿Cuántas veces hay que lanzar la moneda para obtener cuatro caras?

Dos formas de solucionar el problema:

  • Distribución binomial negativa
  • Usando una Cadena de Markov

Usando cadenas de Markov

def T(altura): d = {} for n in [0..altura-1]: d[n] = {(n+1):1/2, n:1/2} d[altura] = {altura: 1} g = DiGraph(d) g.set_pos(dict(( n, (n,0)) for n in [0..altura])) return g g = T(4) g.plot(edge_labels=True, vertex_size = 1000, loop_size=0.15 ).show(figsize=10) 
       

Observamos la evolución de la distribución:

vs = g.vertices() #Matriz de probabilidades de transicion M = matrix(QQ, len(vs)) vs = g.vertices() for v1,v2,p in g.edges(): j1 = vs.index(v1) j2 = vs.index(v2) M[j1,j2] = p pl = [] for k in [0..20]: probs = vector([1] + [0]*(len(vs)-1))*M^k probs /= sum(probs) g.allow_loops(True) pl.append(g.plot(edge_labels=True,loops = True, vertex_size = 500, vertex_colors=dict(((1-p,1-p,1-0.001*random()),[v]) for v,p in zip(vs,probs))) + text('%d : %s'%(k,map(n3,probs)), (2,-0.5), fontsize=12)) animate(pl, axes = False, figsize = 10).show(delay = 50) 
       
/usr/local/sage-5.9/local/lib/python2.7/site-packages/sage/graphs/generi\
c_graph.py:14436: DeprecationWarning: You provided loops as an argument
to a function which has always silently ignored its inputs. This method
may soon be updated so that the method raises an exception instead of
this warning, which will break your code : to be on the safe side,
update it !
See http://trac.sagemath.org/13891 for details.
  return GraphPlot(graph=self, options=options)
/usr/local/sage-5.9/local/lib/python2.7/site-packages/sage/graphs/generic_graph.py:14436: DeprecationWarning: You provided loops as an argument to a function which has always silently ignored its inputs. This method may soon be updated so that the method raises an exception instead of this warning, which will break your code : to be on the safe side, update it !
See http://trac.sagemath.org/13891 for details.
  return GraphPlot(graph=self, options=options)
show(M) print(M.jordan_form()) 
       

[  1|  0   0   0   0]
[---+---------------]
[  0|1/2   1   0   0]
[  0|  0 1/2   1   0]
[  0|  0   0 1/2   1]
[  0|  0   0   0 1/2]
[  1|  0   0   0   0]
[---+---------------]
[  0|1/2   1   0   0]
[  0|  0 1/2   1   0]
[  0|  0   0 1/2   1]
[  0|  0   0   0 1/2]

Observamos que la matriz de transición no es diagonalizable.

Tenemos que construir una nueva función para elevar cualquier matriz de Jordan a una variable simbólica

def potJordan(D, k): '''Eleva una matriz de Jordan (no necesariamente diagonal) a una variable simbolica''' dim = D.ncols() M = matrix(SR, dim, dim) i = 0 while i<=dim-1: #Autovalores simples, o con multiplicidad algebraica == multiplicidad geometrica if (i== dim-1) or D[i,i+1] == 0: M[i,i] = SR(D[i,i])^k i = i + 1 else: #buscamos la dimension del subespacio invariante counter = 1 for j in range(i, dim-1): if D[j,j+1] == 0: break counter += 1 for j in srange(i, i+counter): for t in srange(j,(i+counter)): M[j,t] = binomial(k, t-j)*SR(D[i,i])^(k-(t-j)) i += counter return M 
       
var('k') D,P = M.jordan_form(transformation = True) show(potJordan(D,k)) 
       
matriz = (P*potJordan(D,k)*~P) show(matriz) 
       

Hallamos la media y la varianza

pcum = (P*potJordan(D,k)*~P)[0,-1] p = pcum(k=k) - pcum(k=k-1) p 
       
(1/2)^(k - 1) + 1/48*(k - 3)*(k - 2)*(k - 1)*(1/2)^(k - 4) - 1/48*(k -
2)*(k - 1)*k*(1/2)^(k - 3) + 1/8*(k - 2)*(k - 1)*(1/2)^(k - 3) - 1/8*(k
- 1)*k*(1/2)^(k - 2) + 1/2*(k - 1)*(1/2)^(k - 2) - 1/2*k*(1/2)^(k - 1) -
(1/2)^k
(1/2)^(k - 1) + 1/48*(k - 3)*(k - 2)*(k - 1)*(1/2)^(k - 4) - 1/48*(k - 2)*(k - 1)*k*(1/2)^(k - 3) + 1/8*(k - 2)*(k - 1)*(1/2)^(k - 3) - 1/8*(k - 1)*k*(1/2)^(k - 2) + 1/2*(k - 1)*(1/2)^(k - 2) - 1/2*k*(1/2)^(k - 1) - (1/2)^k
print 'Probabilidad total: ', sum(p,k,4,oo) m = sum(p*k,k,4,oo) print 'Media: ', m print 'Varianza: ', sum(p*(k-m)^2,k,4,oo) 
       
Probabilidad total:  1
Media:  8
Varianza:  8
Probabilidad total:  1
Media:  8
Varianza:  8
 
       
 
       

Distribución Binomial Negativa

Comparamos con la distribución del tiempo de llegada con la distribución de probabilidad binomial negativa

$$f(k) \equiv \Pr(X = k) = {k-1 \choose k-r} (1-p)^r p^{k-r} \quad\text{for }k = r, r+1, r+2, \dots, $$

p.simplify_full() 
       
1/3*(k^3 - 6*k^2 + 11*k - 6)*2^(-k - 1)
1/3*(k^3 - 6*k^2 + 11*k - 6)*2^(-k - 1)
r = 4 p0 = 1/2 p_bn = binomial(k-1,k-r)*p0^(k-r)*(1-p0)^r p_bn.simplify_full() 
       
1/3*(k^3 - 6*k^2 + 11*k - 6)*2^(-k - 1)
1/3*(k^3 - 6*k^2 + 11*k - 6)*2^(-k - 1)
 
       

Paquetes de cereales

En una determinada marca de cereales, con cada paquete viene un regalo con las siguientes condiciones:

  • Hay 2 regalos distintos
  • Todos los regalos tienen la misma posibilidad de aparecer

Queremos saber cuántos paquetes tendremos que comprar para obtener dos veces cada regalo

def posibilidadesRegalosCereales(numRegalos, numVeces): '''Devuelve una lista de los nodos del grafo''' ls = [] ls.append(tuple([0]*numRegalos)) max_suma = numRegalos * numVeces for i in range(1, max_suma + 1): #Hallamos todas las formas de sumar cada numero desde 1 hasta el numero total de regalos, y las permutaciones de cada #suma suma = Partitions(i, max_part = numVeces, max_length = numRegalos).list() for addit in suma: lista = list(addit) if len(lista) < numRegalos: for k in range(len(lista), numRegalos): lista.append(0) perm = Permutations(lista).list() for j in range(len(perm)): ls.append(tuple(perm[j])) return ls def posibilidadesRegalosCereales(numRegalos, numVeces): return [tuple(l) for l in CartesianProduct(*[[0..numVeces]]*numRegalos)] from collections import defaultdict def generaGrafoCereales(numRegalos, numVeces): '''Genera el grafo que modeliza la cadena de Markov''' lista_nodos = posibilidadesRegalosCereales(numRegalos, numVeces) dg = {} for v in lista_nodos: d = defaultdict(QQ) for k in range(numRegalos): vv = list(v) vv[k] += 1 vv = tuple(min(k, numVeces) for k in vv) d[vv] += 1/numRegalos dg[v] = d g = DiGraph(dg) xs = {} hs = {} for v in lista_nodos: h = sum(v) if h in hs: hs[h] += 1 else: hs[h] = 1 xs[v] = hs[h] g.set_pos(dict(( v, (xs[v], -sum(v) )) for v in lista_nodos )) return g 
       
posibilidadesRegalosCereales(3,2) 
       
[(0, 0, 0), (0, 0, 1), (0, 0, 2), (0, 1, 0), (0, 1, 1), (0, 1, 2), (0,
2, 0), (0, 2, 1), (0, 2, 2), (1, 0, 0), (1, 0, 1), (1, 0, 2), (1, 1, 0),
(1, 1, 1), (1, 1, 2), (1, 2, 0), (1, 2, 1), (1, 2, 2), (2, 0, 0), (2, 0,
1), (2, 0, 2), (2, 1, 0), (2, 1, 1), (2, 1, 2), (2, 2, 0), (2, 2, 1),
(2, 2, 2)]
[(0, 0, 0), (0, 0, 1), (0, 0, 2), (0, 1, 0), (0, 1, 1), (0, 1, 2), (0, 2, 0), (0, 2, 1), (0, 2, 2), (1, 0, 0), (1, 0, 1), (1, 0, 2), (1, 1, 0), (1, 1, 1), (1, 1, 2), (1, 2, 0), (1, 2, 1), (1, 2, 2), (2, 0, 0), (2, 0, 1), (2, 0, 2), (2, 1, 0), (2, 1, 1), (2, 1, 2), (2, 2, 0), (2, 2, 1), (2, 2, 2)]
numRegalos = 3 numVeces = 2 g = generaGrafoCereales(numRegalos, numVeces) pl = g.plot(edge_labels=True, vertex_size = numRegalos*numVeces*600, loop_size=0.2) show(pl, axes = False, figsize = 12, xmin = -1, ymin=-numRegalos*numVeces-0.3) 
       
M = matriz_probs(g) 
       
pl = [] vs = g.vertices() for k in [0..20]: probs = vector([1] + [0]*(len(vs)-1))*M^k probs /= sum(probs) g.allow_loops(True) pl.append(g.plot(edge_labels=True,loop_size = 0.2, vertex_size = numRegalos*numVeces*600, vertex_colors=dict(((1-p,1-p,1-0.001*random()),[v]) for v,p in zip(vs,probs)))) animate(pl, axes = False, figsize = 12).show(delay = 50) 
       

Tenemos de nuevo una cadena de Markov absorbente

Calculamos el tiempo medio en conseguir los regalos

Método de las ecuaciones:

Empleamos de nuevo el método de las ecuaciones sabiendo que:

  • Partimos del primer vértice (0,0,0)
  • Llegamos al último vértice (2,2,2)

Observamos la matriz del sistema:

show(matriz_fundamental(g)) 
       
show(vector_fundamental(g)) 
       
ls = zip(g.vertices(), matriz_fundamental(g)\vector_fundamental(g)) print ls[0][0],':',ls[0][1] 
       
(0, 0, 0) : 347/36
(0, 0, 0) : 347/36
 
       

Permutación sin puntos fijos

Un grupo de amigos quieren hacer un  regalo invisible, seleccionando los regalados de la siguiente forma:

  • Cada amigo escoge un nombre al azar
  • En caso de que haya coincidencias, se realiza lo siguiente
    • Si coincide un único nombre, dos o tres, se reparten todos otra vez
    • Si coinciden más, se reparten solamente aquellos que coincidan

Por otro lado, otro grupo de amigos (con el mismo número de personas que el anterior) realiza otro sorteo con las siguientes normas:

  • Cada amigo escoge un nombre al azar
  • En caso de que haya coincidencias, se reparten todos los nombres de nuevo

Queremos saber qué grupo acabará antes de repartir los nombres

numAmigos = 10 
       

GRUPO A

def probabilidadesCoincidencias(k): dicti = {} for i in range(0,k+1): dicti[i] = 0 for p in Permutations(range(0,k)): counter = 0 for i in range(len(p)): if p[i] == i: counter += 1 dicti[counter] = dicti[counter] + 1 for i in range(0,k+1): dicti[i] = dicti[i] / factorial(k) for i in range(1,4): dicti[k] += dicti[i] return dicti 
       

Otra forma de calcular las probabilidades:

http://en.wikipedia.org/wiki/Rencontres_numbers

probabilidadesCoincidencias(numAmigos) 
       
{0: 16481/44800, 1: 16687/45360, 2: 2119/11520, 3: 103/1680, 4: 53/3456,
5: 11/3600, 6: 1/1920, 7: 1/15120, 8: 1/80640, 9: 0, 10: 123607/201600}
{0: 16481/44800, 1: 16687/45360, 2: 2119/11520, 3: 103/1680, 4: 53/3456, 5: 11/3600, 6: 1/1920, 7: 1/15120, 8: 1/80640, 9: 0, 10: 123607/201600}
def probPermutacionesA(size): if size < 4: raise Exception, 'size debe ser >= 4' g = DiGraph() g.allow_loops(True) g.add_vertices(range(size+1)) for i in range(4, size+1): dicc = probabilidadesCoincidencias(i) for j in range(i+1): if dicc[j] != 0: g.add_edge(i,j,dicc[j]) for i in range(1,4): g.delete_vertex(i) g.add_edge(0,0,1) g.delete_vertex(size-1) return g 
       

 

g = probPermutacionesA(numAmigos) g.show(edge_labels = True, layout = 'circular', figsize = 10) 
       

Los números representan las coincidencias. En el caso de 1, 2 o 3 coincidencias, vuelve directamente al nodo 10.

  • $p_k(v_k) = p(\text{1 coincidencia}) + p(\text{2 coincidencias}) + p(\text{3 coincidencias}) + p(\text{k coincidencias})$
show(matriz_probs(g)) 
       

Calculamos la media mediante la matriz elevada a k

M = matriz_probs(g) var('k') D,P = M.jordan_form(transformation = True) show(potJordan(D,k)) matriz = (P*potJordan(D,k)*~P) show(matriz) 
       

pcum = matriz[-1,0] p = pcum(k=k) - pcum(k=k-1) p 
       
-1145/8954*(11/18)^k + 1145/8954*(11/18)^(k - 1) -
105065/1302444*(1177/1920)^k + 105065/1302444*(1177/1920)^(k - 1) +
23348506643631/33686935376116*(123607/201600)^k -
23348506643631/33686935376116*(123607/201600)^(k - 1) +
6960/205931*(619/1008)^k - 6960/205931*(619/1008)^(k - 1) -
1654425/1210858*(5/8)^k + 1654425/1210858*(5/8)^(k - 1) -
7027624/46224477*(19/30)^k + 7027624/46224477*(19/30)^(k - 1)
-1145/8954*(11/18)^k + 1145/8954*(11/18)^(k - 1) - 105065/1302444*(1177/1920)^k + 105065/1302444*(1177/1920)^(k - 1) + 23348506643631/33686935376116*(123607/201600)^k - 23348506643631/33686935376116*(123607/201600)^(k - 1) + 6960/205931*(619/1008)^k - 6960/205931*(619/1008)^(k - 1) - 1654425/1210858*(5/8)^k + 1654425/1210858*(5/8)^(k - 1) - 7027624/46224477*(19/30)^k + 7027624/46224477*(19/30)^(k - 1)
print 'Probabilidad total: ', mi_suma(p,k,1,oo) m = mi_suma(p*k,k,1,oo) print 'Media: ', m, '=', m.n() v = mi_suma(p*(k-m)^2,k,1,oo) print 'Varianza: ', v, ' = ', v.n() 
       
Probabilidad total:  1
Media:  4714793176145/1735740376447 = 2.71630091695857
Varianza:  14015323199879059714359780/3012794654428373272343809  = 
4.65193443545134
Probabilidad total:  1
Media:  4714793176145/1735740376447 = 2.71630091695857
Varianza:  14015323199879059714359780/3012794654428373272343809  =  4.65193443545134

GRUPO B

Tenemos que, en este caso, se puede resolver el problema mediante las fórmulas de la media y la varianza de una distribución geométrica

probs = probabilidadesCoincidencias(numAmigos) p = probs[0] print 'Media = ', 1/p, ' = ', (1/p).n() print 'Varianza = ', (1-p)/(p^2), ' = ', ((1-p)/(p^2)).n() 
       
Media =  44800/16481  =  2.71828165766640
Varianza =  1268691200/271623361  =  4.67077351273921
Media =  44800/16481  =  2.71828165766640
Varianza =  1268691200/271623361  =  4.67077351273921

Obtenemos que, para 10 amigos, la diferencia es casi inapreciable, si bien el tiempo es mayor en el segundo caso.

 
       

Tartaglia

def Tartaglia(altura): d = {} for n in [0..altura-1]: for k in [0..n]: d[(k,n)] = {(k+1,n+1):1/2, (k,n+1):1/2} for k in [0..altura]: d[(k, altura)] = {(k, altura): 1} g = DiGraph(d) g.set_pos(dict(( (k,n), (2*k-n,-n)) for n in [0..altura] for k in [0..n] )) return g N = 4 g = Tartaglia(N) pl = g.plot(edge_labels=True, vertex_size = 2000, loop_size=0.15) pl.show(figsize=10, ymin=-N-1) 
       
N = 4 g = Tartaglia(N) vs = g.vertices() M = matrix(QQ, len(vs)) vs = g.vertices() for v1,v2,l in g.edges(): j1 = vs.index(v1) j2 = vs.index(v2) M[j1,j2] = l for k in [0..5]: probs = vector([1] + [0]*(len(vs)-1))*M^k probs /= sum(probs) pl = g.plot(edge_labels=True, loop_size=0.15 , vertex_size = 2000, vertex_colors=dict(((1-p,1-p,1-0.001*random()),[v]) for v,p in zip(vs,probs))) print ['%s:%s'%(v,p) for v,p in zip(g.vertices(), probs) if p>0] pl.show(axes=False, figsize=10, ymin=-N-1) 
       
['(0, 0):1']

['(0, 1):1/2', '(1, 1):1/2']

['(0, 2):1/4', '(1, 2):1/2', '(2, 2):1/4']

['(0, 3):1/8', '(1, 3):3/8', '(2, 3):3/8', '(3, 3):1/8']

['(0, 4):1/16', '(1, 4):1/4', '(2, 4):3/8', '(3, 4):1/4', '(4, 4):1/16']

['(0, 4):1/16', '(1, 4):1/4', '(2, 4):3/8', '(3, 4):1/4', '(4, 4):1/16']
['(0, 0):1']

['(0, 1):1/2', '(1, 1):1/2']

['(0, 2):1/4', '(1, 2):1/2', '(2, 2):1/4']

['(0, 3):1/8', '(1, 3):3/8', '(2, 3):3/8', '(3, 3):1/8']

['(0, 4):1/16', '(1, 4):1/4', '(2, 4):3/8', '(3, 4):1/4', '(4, 4):1/16']

['(0, 4):1/16', '(1, 4):1/4', '(2, 4):3/8', '(3, 4):1/4', '(4, 4):1/16']

Tenemos que el tiempo en llegar al nivel inferior ¡es siempre 4!

La probabilidad estable se obtiene de la matriz $M^4$

 

 
       
def probEstableTartaglia(altura): g = Tartaglia(altura) vs = g.vertices() M = matrix(QQ, len(vs)) for v1,v2,l in g.edges(): j1 = vs.index(v1) j2 = vs.index(v2) M[j1,j2] = l ls = [1] + [0]*(len(g.vertices())-1) return vector(ls)*M^altura 
       
prob = probEstableTartaglia(N) for i in range(len(g.vertices())): if prob[i] != 0: print g.vertices()[i],':', prob[i] 
       
(0, 4) : 1/16
(1, 4) : 1/4
(2, 4) : 3/8
(3, 4) : 1/4
(4, 4) : 1/16
(0, 4) : 1/16
(1, 4) : 1/4
(2, 4) : 3/8
(3, 4) : 1/4
(4, 4) : 1/16