Unfortunately, there is no implementation of Mixture density networks (MDN) in MxNet yet. And, since MxNet is a community effort, you are more than welcome to contribute!
Migrating code from Keras/TF should be quite straightforward in your case. R bindings for MxNet are quite limited at the moment in terms that it is impossible for now to create custom operations, but looking into the example, I don't see that any custom operations would be required.
I haven't run this code, but here is how MDN model from your example would look like using MxNet Python Symbol API:
def mapping(self, X):
"""pi, mu, sigma = NN(x; theta)"""
hidden1 = mx.sym.FullyConnected(data=X, num_hidden=15) # fully-connected layer with 15 hidden units
act1 = mx.sym.Activation(data=hidden1, act_type='relu')
hidden2 = mx.sym.FullyConnected(data=act1, num_hidden=15) # fully-connected layer with 15 hidden units
act2 = mx.sym.Activation(data=hidden2, act_type='relu')
self.mus = mx.sym.FullyConnected(data=act2, num_hidden=self.K) # fully-connected layer with 15 hidden units
sigma_fc = mx.sym.FullyConnected(data=act2, num_hidden=self.K)
self.sigmas = mx.sym.exp(data=sigma_fc) # the variance
pi_fc = mx.sym.FullyConnected(data=act2, num_hidden=self.K)
self.pi = mx.sym.SoftmaxActivation(data=pi_fc) # the mixture components