Creation of large Numpy tables with broadcasting - how to speed-up

In a python code, I need to compute a huge 4-D table of ~5000x50x100x100 elements.

It contains a normal distribution whose parameters span the different dimensions, and the creation essentially reads like :

result = gaussian( val[np.newaxis, :, np.newaxis, np.newaxis], mean[:, np.newaxis, :, :], sigma[:, np.newaxis, np.newaxis, np.newaxis] )

Where gaussian is obviously :

def gaussian(x, mu, sigma): return np.exp(-0.5 * ((x - mu) / sigma) ** 2) / np.sqrt(2 * np.pi ) / sigma

The sole creation of this table takes very long and slows down the code significantly. Is there any obvious way to speed-up this step ?