<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def poissreg(x0, rho, x, y):
 """ Proximal operator for Poisson regression Computes the proximal operator of the negative log-likelihood loss assumping a Poisson noise distribution. Parameters x0 : array_like The starting or initial point used in the proximal update step rho : float Momentum parameter for the proximal step (larger value -> stays closer to x0) x : (n, k) array_like A design matrix consisting of n examples of k-dimensional features (or input). y : (n,) array_like A vector containing the responses (outupt) to the n features given in x. Returns ------- theta : array_like The parameter vector found after running the proximal update step """


    # objective and gradient
    n = float(x.shape[0])
    f = lambda w: np.mean(np.exp(x.dot(w)) - y * x.dot(w))
    df = lambda w: (x.T.dot(np.exp(x.dot(w))) - x.T.dot(y)) / n

    # minimize via BFGS
    return bfgs(x0, rho, f, df)