Rewrite this code in more optimal way and add comments!

from scipy.special import expit, logsumexp

# because you wrote docstrings in english therefore my comments for code and graphics was in english
def posterior(alpha, beta, L):
    """ Posterior over true labels z p(z|l, \alpha, \beta)
    Args:
        alpha: ndarray of shape (n_experts).
        beta: ndarray of shape (n_problems).
        L: ndarray of shape (n_problems, n_experts).
    Returns:
        ndarray of shape (2, n_problems)
        (2 -- for z = 0 and for z = 1, n_problems -- for each problem in data)
    Tip:
        You may use function log_likelihood here
        See page 7 of seminar, the last formula for details
    """

    z = np.zeros(L.shape[0])
    # compute log_likelihood for z = 0 and z = 1 and then vertical stack
    log_likelihood_z = np.vstack([log_likelihood(alpha, beta, L, z + z_) for z_ in [0,1]])

    # try to explain return code. return log(q(z)) instead of q(z) as it is recommended at help material
    # return log q(z) = log (exp(log_likelihood_z) / sum(exp(log_likelihood_0) + exp(log_likelihood_1))) =
    # = log_likelihood_z - log(sum(exp(log_likelihood_0) + exp(log_likelihood_1)))

    return np.exp(log_likelihood_z - logsumexp(log_likelihood_z,axis=0))


def log_likelihood(alpha, beta, L, z):
    """ p(l=z|z, \alpha, \beta)
    Args:
        alpha: ndarray of shape (n_experts).
        beta: ndarray of shape (n_problems).
        L: ndarray of shape (n_problems, n_experts).
        z: ndarray of shape (n_problems).
    Returns:
        ndarray of shape (n_problems,)
    Tips:
        See page 7 of seminar, the second formula for details
    """
    # make all necessary reshape for broadcast
    z = z.reshape(-1,1)
    beta = beta.reshape(-1,1)
    # compute indicators matrix where element = [l_ij == z_i]
    indicators = L == z
    return np.log(0.5) + np.sum(indicators * -softplus(-beta * alpha) + ~indicators * -softplus(beta * alpha),axis=1)


def alpha_grad_lb(alpha, beta, L, q):
    """ Gradient of lower bound wrt alpha
    Args:
        alpha: ndarray of shape (n_experts).
        beta: ndarray of shape (n_problems).
        L: ndarray of shape (n_problems, n_experts).
        q: ndarray of shape (2, n_problems).
    Returns:
        ndarray of shape (n_experts,)
    Tips:
        See pages 8-9 of seminar for details
    """
    # make all necessary reshape for broadcast
    beta = beta.reshape(-1,1)
    # compute for t = 0 and t = 0
    indicators = L == 0
    # it can be modified for case when q.shape[0] > 2 for example by simple cycle, but here we have requirement for shape q
    t_0 = q[0].reshape(-1,1) * beta * (indicators * expit(-beta * alpha) - ~indicators * expit(beta * alpha))
    t_1 = q[1].reshape(-1,1) * beta * (~indicators * expit(-beta * alpha) - indicators * expit(beta * alpha))

    return np.sum(t_0 + t_1, axis=0)


def logbeta_grad_lb(alpha, beta, L, q):
    """ Gradient of lower bound wrt alpha
    Args:
        alpha: ndarray of shape (n_experts).
        beta: ndarray of shape (n_problems).
        L: ndarray of shape (n_problems, n_experts).
        q: ndarray of shape (2, n_problems).
    Returns:
        ndarray of shape (n_problems,)
    Tips:
        See pages 8-9 of seminar for details
    """
    # make all necessary reshape for broadcast
    beta = beta.reshape(-1,1)
    # compute for t = 0 and t = 0
    indicators = L == 0
    # it can be modified for case when q.shape[0] > 2 for example by simple cycle, but here we have requirement for shape q
    t_0 = q[0].reshape(-1,1) * alpha * (indicators * expit(-beta * alpha) - ~indicators * expit(beta * alpha))
    t_1 = q[1].reshape(-1,1) * alpha * (~indicators * expit(-beta * alpha) - indicators * expit(beta * alpha))
    # * beta - because we need gradient for log(beta)
    return np.sum(t_0 + t_1, axis=1) * beta.reshape(-1)


def lower_bound(alpha, beta, L, q):
    """ Lower bound
    Args:
        alpha: ndarray of shape (n_experts).
        beta: ndarray of shape (n_problems).
        L: ndarray of shape (n_problems, n_experts).
        q: ndarray of shape (2, n_problems).
    Returns:
        single value, number
    Tips:
        You may use function log_likelihood here
        See page 8 of seminar, the fourth formula for details
    """
    # compute for t = 0 and t = 0
    # it can be modified for case when q.shape[0] > 2 for example by simple cycle, but here we have requirement for shape q
    t_0 = q[0] * (log_likelihood(alpha, beta, L, np.zeros(L.shape[0])) - np.log(0.5))
    t_1 = q[1] * (log_likelihood(alpha, beta, L, np.ones(L.shape[0])) - np.log(0.5))

    return np.sum(t_0 + t_1)