def product(p_x, p_y):
    """
    Compute the product distribution p(x,y) = p(x)p(y) given distributions 
    p(x) and p(y)
    """
    p_x = p_x.reshape(1,-1)
    p_y = p_y.reshape(-1,1)
    Pxy = p_y@p_x
    assert Pxy.shape == (p_y.shape[0], p_x.shape[1]) 
    return Pxy