def estimate_gk(y, theta, delta, ck):
    """Helper function to estimate gk from SPSA"""
    # Generate Delta vector
    delta_k = delta()

    # Get the two perturbed values of theta
    # list comprehensions like this are quite nice
    ta = [t + ck * dk for t, dk in zip(theta, delta_k)]
    tb = [t - ck * dk for t, dk in zip(theta, delta_k)]

    # Calculate g_k(theta_k)
    ya, yb = y(ta), y(tb)
    gk = [(ya - yb) / (2 * ck * dk) for dk in delta_k]

    return gk