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}

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}

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}

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

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

M k = P \cdot D k \cdot P - 1

$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)

⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ 00000 2 (1 2) k (1 2) k 000 - (1 2) k + 3 5 (5 6) k - 1 2 (1 2) k + 1 2 (5 6) k 2 5 (5 6) k 3 10 (5 6) k 0 - 2 (1 2) k + 6 5 (5 6) k - (1 2) k + (5 6) k 4 5 (5 6) k 3 5 (5 6) k 0 (1 2) k - 9 5 (5 6) k + 1 1 2 (1 2) k - 3 2 (5 6) k + 1 - 6 5 (5 6) k + 1 - 9 10 (5 6) k + 1 1 ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrr} 0 & 2 \, \left(\frac{1}{2}\right)^{k} & -\left(\frac{1}{2}\right)^{k} + \frac{3}{5} \, \left(\frac{5}{6}\right)^{k} & -2 \, \left(\frac{1}{2}\right)^{k} + \frac{6}{5} \, \left(\frac{5}{6}\right)^{k} & \left(\frac{1}{2}\right)^{k} - \frac{9}{5} \, \left(\frac{5}{6}\right)^{k} + 1 \\ 0 & \left(\frac{1}{2}\right)^{k} & -\frac{1}{2} \, \left(\frac{1}{2}\right)^{k} + \frac{1}{2} \, \left(\frac{5}{6}\right)^{k} & -\left(\frac{1}{2}\right)^{k} + \left(\frac{5}{6}\right)^{k} & \frac{1}{2} \, \left(\frac{1}{2}\right)^{k} - \frac{3}{2} \, \left(\frac{5}{6}\right)^{k} + 1 \\ 0 & 0 & \frac{2}{5} \, \left(\frac{5}{6}\right)^{k} & \frac{4}{5} \, \left(\frac{5}{6}\right)^{k} & -\frac{6}{5} \, \left(\frac{5}{6}\right)^{k} + 1 \\ 0 & 0 & \frac{3}{10} \, \left(\frac{5}{6}\right)^{k} & \frac{3}{5} \, \left(\frac{5}{6}\right)^{k} & -\frac{9}{10} \, \left(\frac{5}{6}\right)^{k} + 1 \\ 0 & 0 & 0 & 0 & 1 \end{array}\right)$

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)

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 :

⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ 00000 1 1 2 000 0 1 6 1 3 1 4 0 0 1 3 2 3 1 2 0 000 1 4 1 ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrr} 0 & 1 & 0 & 0 & 0 \\ 0 & \frac{1}{2} & \frac{1}{6} & \frac{1}{3} & 0 \\ 0 & 0 & \frac{1}{3} & \frac{2}{3} & 0 \\ 0 & 0 & \frac{1}{4} & \frac{1}{2} & \frac{1}{4} \\ 0 & 0 & 0 & 0 & 1 \end{array}\right)$

⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ 00000 1 1 2 000 0 1 6 1 3 1 4 0 0 1 3 2 3 1 2 0 000 1 4 1 ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟

0 :

⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ 1000001000001000001000001 ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrr} 1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \end{array}\right)$

⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ 0000021000 - 2 5 0 2 5 3 10 0 - 4 5 0 4 5 3 5 0 1 5 0 - 1 5 1 10 1 ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrr} 0 & 2 & -\frac{2}{5} & -\frac{4}{5} & \frac{1}{5} \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & \frac{2}{5} & \frac{4}{5} & -\frac{1}{5} \\ 0 & 0 & \frac{3}{10} & \frac{3}{5} & \frac{1}{10} \\ 0 & 0 & 0 & 0 & 1 \end{array}\right)$

-1 :
Traceback (click to the left of this block for traceback)
...
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)

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

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

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

E X i = 1 + p (Ir al nodo a) \cdot E X a + p (Ir al nodo b) \cdot E X b . . . .

$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:

E X i = 0

$EX_i = 0$

para un estado $i$ es absorbente.

El sistema en forma matricial

Sea $E$ el vector:

(E X 1, \dots, E X N)

$(EX_1, \dots, EX_N)$

$E$ resuelve el sistema:

E \cdot A = v

$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

⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ 10000 - 1 1 2 000 0 - 1 6 2 3 - 1 4 0 0 - 1 3 - 2 3 1 2 0 000 - 1 4 1 ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrr} 1 & -1 & 0 & 0 & 0 \\ 0 & \frac{1}{2} & -\frac{1}{6} & -\frac{1}{3} & 0 \\ 0 & 0 & \frac{2}{3} & -\frac{2}{3} & 0 \\ 0 & 0 & -\frac{1}{4} & \frac{1}{2} & -\frac{1}{4} \\ 0 & 0 & 0 & 0 & 1 \end{array}\right)$

(1, 1, 1, 1, 0)

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(1,\,1,\,1,\,1,\,0\right)$

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) \cdot p (Llegar a v 0) + p X (v 1) \cdot p (Llegar a v 1) + . . .

$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_X(v_1), \dots, p_X(v_N))$

$P$ resuelve el sistema:

P \cdot A = v

$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

⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ 1 - 1 3 00 - 1 3 000 - 1 3 10 - 1 3 0000 - 1 3 01 - 1 3 00 - 1 3 0 0 - 1 3 01000 - 1 3 - 1 3 00010 - 1 3 0 0 - 1 3 00 - 1 3 10 - 1 3 0000 - 1 3 01 - 1 3 000 - 1 3 00 - 1 3 1 ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrrrr} 1 & -\frac{1}{3} & -\frac{1}{3} & 0 & -\frac{1}{3} & 0 & 0 & 0 \\ -\frac{1}{3} & 1 & 0 & -\frac{1}{3} & 0 & -\frac{1}{3} & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & -\frac{1}{3} & -\frac{1}{3} & 1 & 0 & 0 & 0 & -\frac{1}{3} \\ -\frac{1}{3} & 0 & 0 & 0 & 1 & -\frac{1}{3} & -\frac{1}{3} & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & -\frac{1}{3} & 0 & -\frac{1}{3} & 0 & 1 & -\frac{1}{3} \\ 0 & 0 & 0 & -\frac{1}{3} & 0 & -\frac{1}{3} & -\frac{1}{3} & 1 \end{array}\right)$

Terminos independientes del sistema

(0, 0, 1, 0, 0, 0, 0, 0)

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(0,\,0,\,1,\,0,\,0,\,0,\,0,\,0\right)$

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 V e n 2 (v i) = 1 - p V e n 1 (v i)

$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

⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ 1 - 1 3 00 - 1 3 000 - 1 3 10 - 1 3 0000 - 1 3 01 - 1 3 00 - 1 3 0 0 - 1 3 01000 - 1 3 - 1 3 00010 - 1 3 0 0 - 1 3 00 - 1 3 10 - 1 3 0000 - 1 3 01 - 1 3 000 - 1 3 00 - 1 3 1 ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟

Terminos independientes del sistema

(1, 1, 0, 1, 1, 0, 1, 1)

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(1,\,1,\,0,\,1,\,1,\,0,\,1,\,1\right)$

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, 13/28), (2, 5/14), (3, 13/28), (4, 15/28), (5, 1/2),
(6, 15/28), (7, 1/2), (8, 1/2), (9, 9/14), ('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

⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ 1 - 1 3 00 - 1 3 0000000 - 1 3 10 - 1 3 00000000 - 1 3 01 - 1 3 00 - 1 3 00000 0 - 1 3 01000 - 1 3 0000 - 1 3 00010 - 1 3 00000 0 - 1 3 00 - 1 3 10 - 1 3 0000 0000 - 1 3 01 - 1 3 000 - 1 2 000 - 1 3 00 - 1 3 10000 0000000010 - 1 - 1 2 00000000 - 1 100 000000000 - 1 10 000000000001 ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrrrrrrrr} 1 & -\frac{1}{3} & -\frac{1}{3} & 0 & -\frac{1}{3} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ -\frac{1}{3} & 1 & 0 & -\frac{1}{3} & 0 & -\frac{1}{3} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & -\frac{1}{3} & -\frac{1}{3} & 1 & 0 & 0 & 0 & -\frac{1}{3} & 0 & 0 & 0 & 0 \\ -\frac{1}{3} & 0 & 0 & 0 & 1 & -\frac{1}{3} & -\frac{1}{3} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & -\frac{1}{3} & 0 & -\frac{1}{3} & 0 & 1 & -\frac{1}{3} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & -\frac{1}{3} & 0 & -\frac{1}{3} & -\frac{1}{3} & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & -1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & -1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & -\frac{1}{2} & 0 & -\frac{1}{2} & 0 & 0 & 1 \end{array}\right)$

Terminos independientes del sistema

(0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0)

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(0,\,0,\,1,\,0,\,0,\,0,\,0,\,0,\,0,\,0,\,0,\,0\right)$

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 V e n 1 (8) = p V e n 1 (9)

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

p V e n 1 (9) = p V e n 1 (10)

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

p V e n 1 (10) = p V e n 1 (8)

$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]

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

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

$p_X(X)=1$

Para otro estado absorbente $v_i$ :

p X (v i) = 0

$p_X(v_i)=0$

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

p X (v i) = 0

$p_X(v_i)=0$

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

p X (v i) = p X (v 0) \cdot p (Llegar a v 0) + p X (v 1) \cdot p (Llegar a v 1) + . . .

$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

⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ 1 - 1 3 00 - 1 3 0000000 - 1 3 10 - 1 3 00000000 - 1 3 01 - 1 3 00 - 1 3 00000 0 - 1 3 01000 - 1 3 0000 - 1 3 00010 - 1 3 00000 0 - 1 3 00 - 1 3 10 - 1 3 0000 0000 - 1 3 01 - 1 3 000 - 1 2 000 - 1 3 00 - 1 3 10000 00000000100 - 1 2 000000000100000000000010000000000001 ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrrrrrrrr} 1 & -\frac{1}{3} & -\frac{1}{3} & 0 & -\frac{1}{3} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ -\frac{1}{3} & 1 & 0 & -\frac{1}{3} & 0 & -\frac{1}{3} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & -\frac{1}{3} & -\frac{1}{3} & 1 & 0 & 0 & 0 & -\frac{1}{3} & 0 & 0 & 0 & 0 \\ -\frac{1}{3} & 0 & 0 & 0 & 1 & -\frac{1}{3} & -\frac{1}{3} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & -\frac{1}{3} & 0 & -\frac{1}{3} & 0 & 1 & -\frac{1}{3} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & -\frac{1}{3} & 0 & -\frac{1}{3} & -\frac{1}{3} & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & -\frac{1}{2} & 0 & -\frac{1}{2} & 0 & 0 & 1 \end{array}\right)$

Terminos independientes

(0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0)

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(0,\,0,\,1,\,0,\,0,\,0,\,0,\,0,\,0,\,0,\,0,\,0\right)$

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

⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ 1 - 1 3 00 - 1 3 0000000 - 1 3 10 - 1 3 00000000 - 1 3 01 - 1 3 00 - 1 3 00000 0 - 1 3 01000 - 1 3 0000 - 1 3 00010 - 1 3 00000 0 - 1 3 00 - 1 3 10 - 1 3 0000 0000 - 1 3 01 - 1 3 000 - 1 2 000 - 1 3 00 - 1 3 10000 00000000100 - 1 2 000000000100000000000010000000000001 ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrrrrrrrr} 1 & -\frac{1}{3} & -\frac{1}{3} & 0 & -\frac{1}{3} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ -\frac{1}{3} & 1 & 0 & -\frac{1}{3} & 0 & -\frac{1}{3} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & -\frac{1}{3} & -\frac{1}{3} & 1 & 0 & 0 & 0 & -\frac{1}{3} & 0 & 0 & 0 & 0 \\ -\frac{1}{3} & 0 & 0 & 0 & 1 & -\frac{1}{3} & -\frac{1}{3} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & -\frac{1}{3} & 0 & -\frac{1}{3} & 0 & 1 & -\frac{1}{3} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & -\frac{1}{3} & 0 & -\frac{1}{3} & -\frac{1}{3} & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & -\frac{1}{2} & 0 & -\frac{1}{2} & 0 & 0 & 1 \end{array}\right)$

Terminos independientes

(0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0)

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(0,\,0,\,0,\,0,\,0,\,1,\,0,\,0,\,0,\,0,\,0,\,0\right)$

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]

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

⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ 1 - 1 3 00 - 1 3 0000000 - 1 3 10 - 1 3 00000000 - 1 3 01 - 1 3 00 - 1 3 00000 0 - 1 3 01000 - 1 3 0000 - 1 3 00010 - 1 3 00000 0 - 1 3 00 - 1 3 10 - 1 3 0000 0000 - 1 3 01 - 1 3 000 - 1 2 000 - 1 3 00 - 1 3 10000 00000000100 - 1 2 000000000100000000000010000000000001 ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrrrrrrrr} 1 & -\frac{1}{3} & -\frac{1}{3} & 0 & -\frac{1}{3} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ -\frac{1}{3} & 1 & 0 & -\frac{1}{3} & 0 & -\frac{1}{3} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & -\frac{1}{3} & -\frac{1}{3} & 1 & 0 & 0 & 0 & -\frac{1}{3} & 0 & 0 & 0 & 0 \\ -\frac{1}{3} & 0 & 0 & 0 & 1 & -\frac{1}{3} & -\frac{1}{3} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & -\frac{1}{3} & 0 & -\frac{1}{3} & 0 & 1 & -\frac{1}{3} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & -\frac{1}{3} & 0 & -\frac{1}{3} & -\frac{1}{3} & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & -\frac{1}{2} & 0 & -\frac{1}{2} & 0 & 0 & 1 \end{array}\right)$

Terminos independientes

(1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1)

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(1,\,1,\,0,\,1,\,1,\,0,\,1,\,1,\,1,\,1,\,1,\,1\right)$

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)

show(M) print(M.jordan_form())

⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ 1 2 0000 1 2 1 2 000 0 1 2 1 2 00 00 1 2 1 2 0 000 1 2 1 ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrr} \frac{1}{2} & \frac{1}{2} & 0 & 0 & 0 \\ 0 & \frac{1}{2} & \frac{1}{2} & 0 & 0 \\ 0 & 0 & \frac{1}{2} & \frac{1}{2} & 0 \\ 0 & 0 & 0 & \frac{1}{2} & \frac{1}{2} \\ 0 & 0 & 0 & 0 & 1 \end{array}\right)$

[  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))

⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ 10000 0 (1 2) k 000 0 k (1 2) (k - 1) (1 2) k 00 0 1 2 (k - 1) k (1 2) (k - 2) k (1 2) (k - 1) (1 2) k 0 0 1 6 (k - 2) (k - 1) k (1 2) (k - 3) 1 2 (k - 1) k (1 2) (k - 2) k (1 2) (k - 1) (1 2) k ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrr} 1 & 0 & 0 & 0 & 0 \\ 0 & \left(\frac{1}{2}\right)^{k} & k \left(\frac{1}{2}\right)^{{\left(k - 1\right)}} & \frac{1}{2} \, {\left(k - 1\right)} k \left(\frac{1}{2}\right)^{{\left(k - 2\right)}} & \frac{1}{6} \, {\left(k - 2\right)} {\left(k - 1\right)} k \left(\frac{1}{2}\right)^{{\left(k - 3\right)}} \\ 0 & 0 & \left(\frac{1}{2}\right)^{k} & k \left(\frac{1}{2}\right)^{{\left(k - 1\right)}} & \frac{1}{2} \, {\left(k - 1\right)} k \left(\frac{1}{2}\right)^{{\left(k - 2\right)}} \\ 0 & 0 & 0 & \left(\frac{1}{2}\right)^{k} & k \left(\frac{1}{2}\right)^{{\left(k - 1\right)}} \\ 0 & 0 & 0 & 0 & \left(\frac{1}{2}\right)^{k} \end{array}\right)$

matriz = (P*potJordan(D,k)*~P) show(matriz)

⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ (1 2) k 0000 1 2 k (1 2) (k - 1) (1 2) k 000 1 8 (k - 1) k (1 2) (k - 2) 1 2 k (1 2) (k - 1) (1 2) k 00 1 48 (k - 2) (k - 1) k (1 2) (k - 3) 1 8 (k - 1) k (1 2) (k - 2) 1 2 k (1 2) (k - 1) (1 2) k 0 - 1 48 (k - 2) (k - 1) k (1 2) (k - 3) - 1 8 (k - 1) k (1 2) (k - 2) - 1 2 k (1 2) (k - 1) - (1 2) k + 1 - 1 8 (k - 1) k (1 2) (k - 2) - 1 2 k (1 2) (k - 1) - (1 2) k + 1 - 1 2 k (1 2) (k - 1) - (1 2) k + 1 - (1 2) k + 1 1 ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrr} \left(\frac{1}{2}\right)^{k} & \frac{1}{2} \, k \left(\frac{1}{2}\right)^{{\left(k - 1\right)}} & \frac{1}{8} \, {\left(k - 1\right)} k \left(\frac{1}{2}\right)^{{\left(k - 2\right)}} & \frac{1}{48} \, {\left(k - 2\right)} {\left(k - 1\right)} k \left(\frac{1}{2}\right)^{{\left(k - 3\right)}} & -\frac{1}{48} \, {\left(k - 2\right)} {\left(k - 1\right)} k \left(\frac{1}{2}\right)^{{\left(k - 3\right)}} - \frac{1}{8} \, {\left(k - 1\right)} k \left(\frac{1}{2}\right)^{{\left(k - 2\right)}} - \frac{1}{2} \, k \left(\frac{1}{2}\right)^{{\left(k - 1\right)}} - \left(\frac{1}{2}\right)^{k} + 1 \\ 0 & \left(\frac{1}{2}\right)^{k} & \frac{1}{2} \, k \left(\frac{1}{2}\right)^{{\left(k - 1\right)}} & \frac{1}{8} \, {\left(k - 1\right)} k \left(\frac{1}{2}\right)^{{\left(k - 2\right)}} & -\frac{1}{8} \, {\left(k - 1\right)} k \left(\frac{1}{2}\right)^{{\left(k - 2\right)}} - \frac{1}{2} \, k \left(\frac{1}{2}\right)^{{\left(k - 1\right)}} - \left(\frac{1}{2}\right)^{k} + 1 \\ 0 & 0 & \left(\frac{1}{2}\right)^{k} & \frac{1}{2} \, k \left(\frac{1}{2}\right)^{{\left(k - 1\right)}} & -\frac{1}{2} \, k \left(\frac{1}{2}\right)^{{\left(k - 1\right)}} - \left(\frac{1}{2}\right)^{k} + 1 \\ 0 & 0 & 0 & \left(\frac{1}{2}\right)^{k} & -\left(\frac{1}{2}\right)^{k} + 1 \\ 0 & 0 & 0 & 0 & 1 \end{array}\right)$

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

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

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 k - r) (1 - p) r p k - r for k = r, r + 1, r + 2, \dots,

$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)

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)

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)]

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))

⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ 100000000000000000000000000 - 1 3 10000000000000000000000000 0 - 1 3 2 3 000000000000000000000000 - 1 3 00100000000000000000000000 0 - 1 3 0 - 1 3 10000000000000000000000 00 - 1 3 0 - 1 3 2 3 000000000000000000000 000 - 1 3 00 2 3 00000000000000000000 0000 - 1 3 0 - 1 3 2 3 0000000000000000000 00000 - 1 3 0 - 1 3 1 3 000000000000000000 - 1 3 00000000100000000000000000 0 - 1 3 0000000 - 1 3 10000000000000000 00 - 1 3 0000000 - 1 3 2 3 000000000000000 000 - 1 3 00000 - 1 3 00100000000000000 0000 - 1 3 00000 - 1 3 0 - 1 3 10000000000000 00000 - 1 3 00000 - 1 3 0 - 1 3 2 3 000000000000 000000 - 1 3 00000 - 1 3 00 2 3 00000000000 0000000 - 1 3 00000 - 1 3 0 - 1 3 2 3 0000000000 00000000 - 1 3 00000 - 1 3 0 - 1 3 1 3 000000000 000000000 - 1 3 00000000 2 3 00000000 0000000000 - 1 3 0000000 - 1 3 2 3 0000000 00000000000 - 1 3 0000000 - 1 3 1 3 000000 000000000000 - 1 3 00000 - 1 3 00 2 3 00000 0000000000000 - 1 3 00000 - 1 3 0 - 1 3 2 3 0000 00000000000000 - 1 3 00000 - 1 3 0 - 1 3 1 3 000 000000000000000 - 1 3 00000 - 1 3 00 1 3 00 0000000000000000 - 1 3 00000 - 1 3 0 - 1 3 1 3 0 00000000000000000 - 1 3 00000 - 1 3 0 - 1 3 1 ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrrrrrrrrrrrrrrrrrrrrrrr} 1 & -\frac{1}{3} & 0 & -\frac{1}{3} & 0 & 0 & 0 & 0 & 0 & -\frac{1}{3} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & -\frac{1}{3} & 0 & -\frac{1}{3} & 0 & 0 & 0 & 0 & 0 & -\frac{1}{3} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \frac{2}{3} & 0 & 0 & -\frac{1}{3} & 0 & 0 & 0 & 0 & 0 & -\frac{1}{3} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & -\frac{1}{3} & 0 & -\frac{1}{3} & 0 & 0 & 0 & 0 & 0 & -\frac{1}{3} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & -\frac{1}{3} & 0 & -\frac{1}{3} & 0 & 0 & 0 & 0 & 0 & -\frac{1}{3} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{2}{3} & 0 & 0 & -\frac{1}{3} & 0 & 0 & 0 & 0 & 0 & -\frac{1}{3} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & \frac{2}{3} & -\frac{1}{3} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -\frac{1}{3} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{2}{3} & -\frac{1}{3} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -\frac{1}{3} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{1}{3} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -\frac{1}{3} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & -\frac{1}{3} & 0 & -\frac{1}{3} & 0 & 0 & 0 & 0 & 0 & -\frac{1}{3} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & -\frac{1}{3} & 0 & -\frac{1}{3} & 0 & 0 & 0 & 0 & 0 & -\frac{1}{3} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{2}{3} & 0 & 0 & -\frac{1}{3} & 0 & 0 & 0 & 0 & 0 & -\frac{1}{3} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & -\frac{1}{3} & 0 & -\frac{1}{3} & 0 & 0 & 0 & 0 & 0 & -\frac{1}{3} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & -\frac{1}{3} & 0 & -\frac{1}{3} & 0 & 0 & 0 & 0 & 0 & -\frac{1}{3} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{2}{3} & 0 & 0 & -\frac{1}{3} & 0 & 0 & 0 & 0 & 0 & -\frac{1}{3} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{2}{3} & -\frac{1}{3} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -\frac{1}{3} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{2}{3} & -\frac{1}{3} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -\frac{1}{3} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{1}{3} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -\frac{1}{3} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{2}{3} & -\frac{1}{3} & 0 & -\frac{1}{3} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{2}{3} & -\frac{1}{3} & 0 & -\frac{1}{3} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{1}{3} & 0 & 0 & -\frac{1}{3} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{2}{3} & -\frac{1}{3} & 0 & -\frac{1}{3} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{2}{3} & -\frac{1}{3} & 0 & -\frac{1}{3} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{1}{3} & 0 & 0 & -\frac{1}{3} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{1}{3} & -\frac{1}{3} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{1}{3} & -\frac{1}{3} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{array}\right)$

show(vector_fundamental(g))

(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0)

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(1,\,1,\,1,\,1,\,1,\,1,\,1,\,1,\,1,\,1,\,1,\,1,\,1,\,1,\,1,\,1,\,1,\,1,\,1,\,1,\,1,\,1,\,1,\,1,\,1,\,1,\,0\right)$

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

(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}

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))

⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ 1 3 8 11 30 53 144 103 280 2119 5760 16481 44800 0 5 8 0 1 48 1 72 1 64 53 3456 00 19 30 0 1 240 1 360 11 3600 000 11 18 0 1 1440 1 1920 0000 619 1008 0 1 15120 00000 1177 1920 1 80640 000000 123607 201600 ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrrr} 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ \frac{3}{8} & \frac{5}{8} & 0 & 0 & 0 & 0 & 0 \\ \frac{11}{30} & 0 & \frac{19}{30} & 0 & 0 & 0 & 0 \\ \frac{53}{144} & \frac{1}{48} & 0 & \frac{11}{18} & 0 & 0 & 0 \\ \frac{103}{280} & \frac{1}{72} & \frac{1}{240} & 0 & \frac{619}{1008} & 0 & 0 \\ \frac{2119}{5760} & \frac{1}{64} & \frac{1}{360} & \frac{1}{1440} & 0 & \frac{1177}{1920} & 0 \\ \frac{16481}{44800} & \frac{53}{3456} & \frac{11}{3600} & \frac{1}{1920} & \frac{1}{15120} & \frac{1}{80640} & \frac{123607}{201600} \end{array}\right)$

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)

⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ 1000000 0 (19 30) k 00000 00 (5 8) k 0000 000 (619 1008) k 000 0000 (123607 201600) k 00 00000 (1177 1920) k 0 000000 (11 18) k ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrrr} 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & \left(\frac{19}{30}\right)^{k} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \left(\frac{5}{8}\right)^{k} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \left(\frac{619}{1008}\right)^{k} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & \left(\frac{123607}{201600}\right)^{k} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & \left(\frac{1177}{1920}\right)^{k} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & \left(\frac{11}{18}\right)^{k} \end{array}\right)$

⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ 1 - (5 8) k + 1 - (19 30) k + 1 1 2 (11 18) k - 3 2 (5 8) k + 1 522 1067 (619 1008) k - 14 11 (5 8) k - 21 97 (19 30) k + 1 - 2 11 (11 18) k + 21013 29601 (1177 1920) k - 32 23 (5 8) k - 16 117 (19 30) k + 1 - 1145 8954 (11 18) k - 105065 1302444 (1177 1920) k + 23348506643631 33686935376116 (123607 201600) k + 6960 205931 (619 1008) k - 1654425 1210858 (5 8) k - 7027624 46224477 (19 30) k + 1 0 (5 8) k 0 - 3 2 (11 18) k + 3 2 (5 8) k - 14 11 (619 1008) k + 14 11 (5 8) k 6 11 (11 18) k - 490 253 (1177 1920) k + 32 23 (5 8) k 3435 8954 (11 18) k + 1225 5566 (1177 1920) k - 23349838925 12406187838 (123607 201600) k - 560 6369 (619 1008) k + 1654425 1210858 (5 8) k 00 (19 30) k 0 - 21 97 (619 1008) k + 21 97 (19 30) k - 16 117 (1177 1920) k + 16 117 (19 30) k 20 1287 (1177 1920) k - 1319668 8646979 (123607 201600) k - 280 18721 (619 1008) k + 7027624 46224477 (19 30) k 000 (11 18) k 0 - 4 11 (11 18) k + 4 11 (1177 1920) k - 1145 4477 (11 18) k - 5 121 (1177 1920) k + 1330 4477 (123607 201600) k 0000 (619 1008) k 0 - 40 579 (123607 201600) k + 40 579 (619 1008) k 00000 (1177 1920) k - 5 44 (1177 1920) k + 5 44 (123607 201600) k 000000 (123607 201600) k ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrrr} 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ -\left(\frac{5}{8}\right)^{k} + 1 & \left(\frac{5}{8}\right)^{k} & 0 & 0 & 0 & 0 & 0 \\ -\left(\frac{19}{30}\right)^{k} + 1 & 0 & \left(\frac{19}{30}\right)^{k} & 0 & 0 & 0 & 0 \\ \frac{1}{2} \, \left(\frac{11}{18}\right)^{k} - \frac{3}{2} \, \left(\frac{5}{8}\right)^{k} + 1 & -\frac{3}{2} \, \left(\frac{11}{18}\right)^{k} + \frac{3}{2} \, \left(\frac{5}{8}\right)^{k} & 0 & \left(\frac{11}{18}\right)^{k} & 0 & 0 & 0 \\ \frac{522}{1067} \, \left(\frac{619}{1008}\right)^{k} - \frac{14}{11} \, \left(\frac{5}{8}\right)^{k} - \frac{21}{97} \, \left(\frac{19}{30}\right)^{k} + 1 & -\frac{14}{11} \, \left(\frac{619}{1008}\right)^{k} + \frac{14}{11} \, \left(\frac{5}{8}\right)^{k} & -\frac{21}{97} \, \left(\frac{619}{1008}\right)^{k} + \frac{21}{97} \, \left(\frac{19}{30}\right)^{k} & 0 & \left(\frac{619}{1008}\right)^{k} & 0 & 0 \\ -\frac{2}{11} \, \left(\frac{11}{18}\right)^{k} + \frac{21013}{29601} \, \left(\frac{1177}{1920}\right)^{k} - \frac{32}{23} \, \left(\frac{5}{8}\right)^{k} - \frac{16}{117} \, \left(\frac{19}{30}\right)^{k} + 1 & \frac{6}{11} \, \left(\frac{11}{18}\right)^{k} - \frac{490}{253} \, \left(\frac{1177}{1920}\right)^{k} + \frac{32}{23} \, \left(\frac{5}{8}\right)^{k} & -\frac{16}{117} \, \left(\frac{1177}{1920}\right)^{k} + \frac{16}{117} \, \left(\frac{19}{30}\right)^{k} & -\frac{4}{11} \, \left(\frac{11}{18}\right)^{k} + \frac{4}{11} \, \left(\frac{1177}{1920}\right)^{k} & 0 & \left(\frac{1177}{1920}\right)^{k} & 0 \\ -\frac{1145}{8954} \, \left(\frac{11}{18}\right)^{k} - \frac{105065}{1302444} \, \left(\frac{1177}{1920}\right)^{k} + \frac{23348506643631}{33686935376116} \, \left(\frac{123607}{201600}\right)^{k} + \frac{6960}{205931} \, \left(\frac{619}{1008}\right)^{k} - \frac{1654425}{1210858} \, \left(\frac{5}{8}\right)^{k} - \frac{7027624}{46224477} \, \left(\frac{19}{30}\right)^{k} + 1 & \frac{3435}{8954} \, \left(\frac{11}{18}\right)^{k} + \frac{1225}{5566} \, \left(\frac{1177}{1920}\right)^{k} - \frac{23349838925}{12406187838} \, \left(\frac{123607}{201600}\right)^{k} - \frac{560}{6369} \, \left(\frac{619}{1008}\right)^{k} + \frac{1654425}{1210858} \, \left(\frac{5}{8}\right)^{k} & \frac{20}{1287} \, \left(\frac{1177}{1920}\right)^{k} - \frac{1319668}{8646979} \, \left(\frac{123607}{201600}\right)^{k} - \frac{280}{18721} \, \left(\frac{619}{1008}\right)^{k} + \frac{7027624}{46224477} \, \left(\frac{19}{30}\right)^{k} & -\frac{1145}{4477} \, \left(\frac{11}{18}\right)^{k} - \frac{5}{121} \, \left(\frac{1177}{1920}\right)^{k} + \frac{1330}{4477} \, \left(\frac{123607}{201600}\right)^{k} & -\frac{40}{579} \, \left(\frac{123607}{201600}\right)^{k} + \frac{40}{579} \, \left(\frac{619}{1008}\right)^{k} & -\frac{5}{44} \, \left(\frac{1177}{1920}\right)^{k} + \frac{5}{44} \, \left(\frac{123607}{201600}\right)^{k} & \left(\frac{123607}{201600}\right)^{k} \end{array}\right)$

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)

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

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

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']

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