Interesting stuff. I have no solution but here are some things to play around with that I think at least gives some component to a solution; how to calculate ecdf and how to autograd through the unautogradable 
model = torch.nn.Sequential(torch.nn.Linear(1,10),
                            torch.nn.Tanh(),
                            torch.nn.Linear(10,10),
                            torch.nn.Tanh(),
                            torch.nn.Linear(10,1))

losses = []
optimizer = torch.optim.RMSprop(model.parameters(), lr=1e-4)
for i in range(10000):
    def closure():
        x = Variable(torch.randn([500,1]), requires_grad = False)
        y = model(x)

        x_sorted, x_argsort = x.sort(0)
        y_sorted, y_argsort = y.sort(0)

        ecdf_actual = x_argsort.float()/x_argsort.max().float()
        ecdf_pred = y_argsort.float()/y_argsort.max().float()

        ecdf_pred = Variable(ecdf_pred.data-y.data)+y

        loss = ((ecdf_actual-ecdf_pred)**2).mean()
        optimizer.zero_grad()
        loss.backward()
        losses.append(loss.data.numpy())
        return loss
    optimizer.step(closure)

plt.semilogy(losses)
losses = np.array(losses).flatten().cumsum()
plt.semilogy(losses/np.ones_like(losses).cumsum())

It actually seems to learn to map x s.t the mappings ecdf corresponds to that of the ecdf of x (granted, which a single linear layer would) but I imagine your target value is different.
This is probably not a case for Distributions unless you think you’ll estimate the ecdf with some parametric distribution