Rendimiento de iteración

Hice una función para evaluar experimentalmente el siguiente problema, tomado de un Manual básico para las Matemáticas de Ingeniería Financiera .

Problema : Sea X el número de veces que debes lanzar una moneda justa hasta que caiga cara. ¿Qué son E [X] (valor esperado) y var (X) (varianza)?

Siguiendo la solución del libro de texto, el siguiente código da la respuesta correcta:

from sympy import * k = symbols('k') Expected_Value = summation(k/2**k, (k, 1, oo)) # Both solutions work Variance = summation(k**2/2**k, (k, 1, oo)) - Expected_Value**2 

Para validar esta respuesta, decidí intentar hacer una función para simular este experimento. El siguiente código es lo que se me ocurrió.

 def coin_toss(toss, probability=[0.5, 0.5]): """Computes expected value and variance for coin toss experiment""" flips = [] # Collects total number of throws until heads appear per experiment. for _ in range(toss): # Simulate n flips number_flips=[] # Number of flips until heads is tossed while sum(number_flips) == 0: # Continue simulation while Tails are thrown number_flips.append(np.random.choice(2, p=probability)) # Append result to number_flips flips.append(len(number_flips)) #Append number of flips until lands heads to flips Expected_Value, Variance = np.mean(flips), np.var(flips) print('E[X]: {}'.format(Expected_Value), '\nvar[X]: {}'.format(Variance)) # Return expected value 

El tiempo de ejecución si simulo experimentos 1e6, utilizando el siguiente código es de aproximadamente 35.9 segundos.

  from timeit import Timer t1 = Timer("""coin_toss(1000000)""", """from __main__ import coin_toss""") print(t1.timeit(1)) 

En el interés de desarrollar mi comprensión de Python, ¿es esta una forma particularmente eficiente / pythonica de abordar un problema como este? ¿Cómo puedo utilizar las bibliotecas existentes para mejorar la eficiencia / ejecución de flujo?

Para poder codificar de forma eficiente y pirónica, debes echar un vistazo a PythonSpeed y NumPy . Un ejemplo de un código más rápido con numpy se puede encontrar a continuación.

El abc de la optimización en python + numpy es vectorizar las operaciones, que en este caso es bastante difícil porque hay un while que podría ser infinito, la moneda se puede voltear las colas 40 veces seguidas. Sin embargo, en lugar de hacer un for con iteraciones de toss , el trabajo se puede hacer en trozos . Esa es la principal diferencia entre coin_toss de la pregunta y este enfoque coin_toss_2d .

coin_toss_2d

La principal coin_toss_2d de coin_toss_2d es trabajar con trozos, el tamaño de estos trozos tiene algunos valores predeterminados, pero se pueden modificar (ya que afectarán la velocidad). Por lo tanto, solo iterará durante el curso while current_toss varias veces toss%flips_at_a_time . Esto se logra con numpy, que permite generar una matriz con los resultados de repetir flips_at_a_time veces el experimento de lanzar una moneda flips_per_try veces. Esta matriz contendrá 0 (colas) y 1 (cabezas).

 # ie doing only 5 repetitions with 3 flips_at_a_time flip_events = np.random.choice([0,1],size=(repetitions_at_a_time,flips_per_try),p=probability) # Out [[0 0 0] # still no head, we will have to keep trying [0 1 1] # head at the 2nd try (position 1 in python) [1 0 0] [1 0 1] [0 0 1]] 

Una vez obtenido este resultado, se llama argmax . Este encuentra el índice correspondiente al máximo (que será 1, una cabeza) de cada fila (repetición) y, en caso de ocurrencias múltiples, devuelve la primera, que es exactamente lo que se necesita, la primera cabeza después de una secuencia de colas. .

 maxs = flip_events.argmax(axis=1) # Out [0 1 0 0 2] # The first position is 0, however, flip_events[0,0]!=1, it's not a head! 

Sin embargo, el caso donde toda la fila es 0 debe ser considerado. En este caso, el máximo será 0, y su primera aparición también será 0, la primera columna (bash). Por lo tanto, verificamos que todos los máximos encontrados en el primer bash correspondan a una cabeza en el primer bash.

 not_finished = (maxs==0) & (flip_events[:,0]!=1) # Out [ True False False False False] # first repetition is not finished 

Si ese no es el caso, repetimos en bucle ese mismo proceso, pero solo para las repeticiones donde no hubo ninguna pista en ninguno de los bashs.

 n = np.sum(not_finished) while n!=0: # while there are sequences without any head flip_events = np.random.choice([0,1],size=(n,flips_per_try),p=probability) # number of experiments reduced to n (the number of all tails sequences) maxs2 = flip_events.argmax(axis=1) maxs[not_finished] += maxs2+flips_per_try # take into account that there have been flips_per_try tries already (every iteration is added) not_finished2 = (maxs2==0) & (flip_events[:,0]!=1) not_finished[not_finished] = not_finished2 n = np.sum(not_finished) # Out # flip_events [[1 0 1]] # Now there is a head # maxs2 [0] # maxs [3 1 0 0 2] # The value of the still unfinished repetition has been updated, # taking into account that the first position in flip_events is the 4th, # without affecting the rest 

Luego se almacenan los índices correspondientes a la primera aparición de cabecera (tenemos que agregar 1 porque la indexación de Python comienza en cero en lugar de 1). Hay un bloque de try ... except ... para hacer frente a los casos en que el toss no es un múltiplo de repetitions_at_a_time .

 def coin_toss_2d(toss, probability=[.5,.5],repetitions_at_a_time=10**5,flips_per_try=20): # Initialize and preallocate data current_toss = 0 flips = np.empty(toss) # loop by chunks while current_toss 

coin_toss_map

Aquí el código es básicamente el mismo, pero ahora, intrinsec while se realiza en una función separada, a la que se llama desde la función wrapper coin_toss_map usando map .

 def toss_chunk(args): probability,repetitions_at_a_time,flips_per_try = args # repeat repetitions_at_a_time times experiment "flip coin flips_per_try times" flip_events = np.random.choice([0,1],size=(repetitions_at_a_time,flips_per_try),p=probability) # store first head ocurrence maxs = flip_events.argmax(axis=1) # Check for all tails sequences not_finished = (maxs==0) & (flip_events[:,0]!=1) n = np.sum(not_finished) while n!=0: # while there are sequences without any head flip_events = np.random.choice([0,1],size=(n,flips_per_try),p=probability) # number of experiments reduced to n (the number of all tails sequences) maxs2 = flip_events.argmax(axis=1) maxs[not_finished] += maxs2+flips_per_try # take into account that there have been flips_per_try tries already (every iteration is added) not_finished2 = (maxs2==0) & (flip_events[:,0]!=1) not_finished[not_finished] = not_finished2 n = np.sum(not_finished) return maxs+1 def coin_toss_map(toss,probability=[.5,.5],repetitions_at_a_time=10**5,flips_per_try=20): n_chunks, remainder = divmod(toss,repetitions_at_a_time) args = [(probability,repetitions_at_a_time,flips_per_try) for _ in range(n_chunks)] if remainder: args.append((probability,remainder,flips_per_try)) flips = np.concatenate(map(toss_chunk,args)) # Once all values are obtained, average and return them Expected_Value, Variance = np.mean(flips), np.var(flips) return Expected_Value, Variance 

Comparación de rendimiento

En mi computadora, tengo el siguiente tiempo de cómputo:

 In [1]: %timeit coin_toss(10**6) # Out # ('E[X]: 2.000287', '\nvar[X]: 1.99791891763') # ('E[X]: 2.000459', '\nvar[X]: 2.00692478932') # ('E[X]: 1.998118', '\nvar[X]: 1.98881045808') # ('E[X]: 1.9987', '\nvar[X]: 1.99508631') # 1 loop, best of 3: 46.2 s per loop In [2]: %timeit coin_toss_2d(10**6,repetitions_at_a_time=5*10**5,flips_per_try=4) # Out # 1 loop, best of 3: 197 ms per loop In [3]: %timeit coin_toss_map(10**6,repetitions_at_a_time=4*10**5,flips_per_try=4) # Out # 1 loop, best of 3: 192 ms per loop 

Y los resultados para la media y la varianza son:

 In [4]: [coin_toss_2d(10**6,repetitions_at_a_time=10**5,flips_per_try=10) for _ in range(4)] # Out # [(1.999848, 1.9990739768960009), # (2.000654, 2.0046035722839997), # (1.999835, 2.0072329727749993), # (1.999277, 2.001566477271)] In [4]: [coin_toss_map(10**6,repetitions_at_a_time=10**5,flips_per_try=4) for _ in range(4)] # Out # [(1.999552, 2.0005057992959996), # (2.001733, 2.011159996711001), # (2.002308, 2.012128673136001), # (2.000738, 2.003613455356)]