¿Se puede vectorizar este código numpy?

Question

Escribí la siguiente función para producir n realizaciones de un proceso CIR para un conjunto dado de parámetros:

def cir_simulations(alpha, mu, sigma, delta_t, n, num_sims):

    x = np.reshape(np.array([mu] * num_sims), (-1, 1))

    for i in range(0, n):
        x = np.concatenate((x, np.reshape(x[:, -1], (-1, 1)) + alpha * (
                    mu - np.reshape(x[:, -1], (-1, 1))) * delta_t + sigma * np.sqrt(
            np.reshape(x[:, -1], (-1, 1))) * np.sqrt(delta_t) * np.random.normal(0, 1, size=(num_sims, 1))), axis=1)

    return x

Este código funciona, pero ahora me pregunto si sería posible eliminar el bucle y vectorizarlo por completo. He tenido problemas para encontrar una manera de hacer esto, ya que la operación que se realiza en el ciclo es recursiva (los valores en la siguiente columna de la matriz no dependen linealmente de los valores de las columnas anteriores).

Además, ¿podría simplificarse este código? En particular, siento que puede haber una complejidad innecesaria en la forma en que accedo a la última columna de la matriz usando

np.reshape(x[:, -1], (-1, 1))

Answer 1

Dado que la salida de cada paso depende del paso anterior, es imposible eliminar por completo el ciclo a menos que la fórmula se pueda simplificar para permitir un cálculo de varios pasos. Sin embargo, el código definitivamente se puede mejorar.

• La llamada np.concatenatees generalmente ineficaz ya que requiere reasignación de memoria y copia de datos. Siempre que sea posible, se debe preasignar toda la memoria necesaria y actualizarla iterativamente.

• Algunos cálculos se pueden sacar del ciclo para mejorar la eficiencia. Por ejemplo, la np.random.normalllamada se puede realizar solo una vez. Lo mismo con alpha * delta_t.

• Se podrían elegir nombres de variables más significativos.

• Siempre es una buena práctica tener cadenas de documentos y sugerencias de tipo en su código.

Aquí hay una versión del código mejorado:

def cir_simulations(alpha: float, mu: float, sigma: float, delta_t: float, sim_steps: int, num_sims: int):
    """
    Simulate the CIR process.

    Parameters:
    Input
    ...

    Output
    ...
    """
    output_shape = (num_sims, sim_steps + 1)
    sim_results = np.empty(output_shape)
    sigma_dW = np.random.normal(0, sigma * np.sqrt(delta_t), size=output_shape) 
    alpha_dt = alpha * delta_t
    sim_results[0, :] = r_prev = mu

    for r_t, sigma_dWt in zip(sim_results[1:], sigma_dW):
        r_t[:] = r_prev = (mu - r_prev) * alpha_dt + np.sqrt(r_prev) * sigma_dWt

    return sim_results

No estoy muy familiarizado con el significado de todos los parámetros de la fórmula. Hay dos cosas que tengo dudas en tu cálculo:

• La tasa \$r\$en la fórmula se inicializa con el mismo parámetro muutilizado en el cálculo iterativo. ¿Es esperado?
• El incremento \$dW_t\$utilizó una varianza fija de 1. ¿Se supone que este es el caso o debería ser otro parámetro de función?