Otra manera de acercarse a este, más accesibles para los estudiantes de escuela secundaria que las cadenas de Markov, es a través de la simulación. En lugar de tratar de resolver el problema exactamente, simular la propagación de infecciones, muchas veces en un equipo y calcular el promedio de tiempo para infectar a todos los pueblos. He aquí un pequeño script en python que hace esto.
import random
a,b,c,d,e,f,g,h,i = 'A B C D E F G H I'.split()
edges = {a:{b},b:{a,e},c:{e,f},d:{e,g},e:{b,c,d,h,i},
f:{c,i},g:{d,h},h:{e,g,i},i:{e,f,h}}
for x,y in edges.items():
for z in y: assert x in edges[z]
random.seed()
universe={a,b,c,d,e,f,g,h,i}
trials = 1000
average = 0
for _ in range(trials):
infected = {a}
days = 0
while infected != universe:
days += 1
roads = {(town,x) for town in infected for x in edges[town]}
for road in roads:
if random.random()<.1:
infected.add(road[1])
average += days
print(average/trials)
Me encontré con esto de una vez y consiguió $37.927.$ Hay varias posibles mejoras. En primer lugar, ejecute más de $1000$ ensayos. Segundo, hacer un histograma de los resultados, en lugar de calcular la media. En tercer lugar, se puede calcular el error estándar, para calcular un intervalo de confianza para su estimación.
Ya que el tamaño de la matriz de transición crece exponencialmente con el número de ciudades, me imagino que si algo como esto se hace realmente en epidemiología, es más probable hecho con la simulación de cadenas de Markov, aunque estoy seguro de que la simulación sería mucho más sofisticados que los de este.
EDITAR
Con 100.000 ensayos, yo siempre obtienen un promedio de $38$, muy diferente de fanvacoolt la respuesta. No puedo ver ningún error en mi código, así que me voy a tener que probar la cadena de Markov enfoque-pero no esta noche.
EDITAR
Escribí un guión para calcular el número esperado de días hasta que todos los pueblos que están infectados, el uso de cadenas de Markov y la matriz fundamental. Se produce un resultado en perfecto acuerdo con las simulaciones, por lo que debe ser correcto. Aquí está la secuencia de comandos:
import numpy as np
from itertools import combinations
from itertools import chain
a,b,c,d,e,f,g,h,i = 'A B C D E F G H I'.split()
edges = {a:{b},b:{a,e},c:{e,f},d:{e,g},e:{b,c,d,h,i},
f:{c,i},g: {d,h},h:{e,g,i},i:{e,f,h}}
universe={a,b,c,d,e,f,g,h,i}
for x,y in edges.items():
for z in y: assert x in edges[z]
weight={a:1,b:2,c:4,d:8,e:16,f:32,g:64,h:128,i:256}
def powerset(iterable):
"powerset([1,2,3]) --> () (1,) (2,) (3,) (1,2) (1,3) (2,3) (1,2,3)"
s = set(iterable)
return chain.from_iterable(combinations(s, r) for r in range(len(s)+1))
def index(s):
"convert a set of towns to an array index"
return sum(weight[town] for town in s)
def transitionMatrix():
"construct and return the trnsition matrix"
P=np.zeros((512,512)) # transition matrix
for old in powerset(universe):
idx = index(old)
roads = {(town,x) for town in old for x in edges[town] if x not in old}
r = len(roads)
for t in powerset(roads):
k=len(t)
prob = .1**k*.9**(r-k)
spread = {x[1] for x in t}
new = set(old) | set(spread)
jdx=index(new)
P[idx, jdx] += prob
return P
P = transitionMatrix()
Q= P[1:511,1:511] # transitions between transient states
del(P)
N = np.linalg.inv(np.eye(510)-Q) #fundamental matrix
N=N@np.ones((510))
def steps(s):
"Average number of days until all twons infected starting from set s"
return N[index(s)-1] # row 0 was deleted from transition matrix
print(steps({'A'}))
Ejecuta esto produce 38.0376337969. Incluso con sólo $1000$ ensayos, la simulación de la secuencia de comandos produce consistentemente cerca de los valores de $38$.
UNA ÚLTIMA COSA
Si consideramos sólo a los estados que son realmente alcanzable empezar con sólo Una ciudad infectada, hay 51 estados. Aquí es un script de python que se calcula.
from itertools import combinations
from itertools import chain
a,b,c,d,e,f,g,h,i = 'A B C D E F G H I'.split()
edges = {a:{b},b:{a,e},c:{e,f},d:{e,g},e:{b,c,d,h,i},
f:{c,i},g:{d,h},h:{e,g,i},i:{e,f,h}}
universe={a,b,c,d,e,f,g,h,i}
def powerset(iterable):
"powerset([1,2,3]) --> () (1,) (2,) (3,) (1,2) (1,3) (2,3) (1,2,3)"
s = set(iterable)
return chain.from_iterable(combinations(s, r) for r in range(len(s)+1))
def possible(s):
'''
s is a possible state if it includes A and the induced subgraph is connected
'''
if a not in s: return False
Q = {a}
visited = set()
while Q:
x = Q.pop()
visited.add(x)
for y in edges[x]:
if y in s and y not in visited:
Q.add(y)
return len(s)==len(visited)
possibilities = []
for s in powerset(universe):
if possible(s): possibilities.append(s)
print(len(possibilities))
for p in possibilities:print(''.join(sorted(p)))
