I had implemented a Normal distribution as a torch.autograd.Function:
import torch
from torch.autograd import Variable,Function
import torch.optim as optim
from scipy import stats
import numpy as np

dtype = torch.FloatTensor
def get_tau_sd(tau=None, sd=None):
    """
    Function port from PyMC3

    Find precision and standard deviation
    .. math::
        \tau = \frac{1}{\sigma^2}
    Parameters
    ----------
    tau : array-like, optional
    sd : array-like, optional
    Results
    -------
    Returns tuple (tau, sd)
    Notes
    -----
    If neither tau nor sd is provided, returns (1., 1.)
    """
    if tau is None:
        if sd is None:
            sd = 1.
            tau = 1.
        else:
            tau = sd**-2.

    else:
        if sd is not None:
            raise ValueError("Can't pass both tau and sd")
        else:
            sd = tau**-.5

    # cast tau and sd to float in a way that works for both np.arrays
    # and pure python
    tau = 1. * tau
    sd = 1. * sd

    return (tau, sd)

class Normal(Function):
    
    def __init__(self, mu=0, sd=None, tau=None, **kwargs):
        """Construct Normal distributions with mean and stddev `loc` and `scale`.
        The parameters `loc` and `scale` must be shaped in a way that supports
        broadcasting (e.g. `loc + scale` is a valid operation).
        Args:
          mu: Floating point tensor; the means of the distribution(s).
          sd: Floating point tensor; the stddevs of the distribution(s).
            Must contain only positive values.
            Raises:
          TypeError: if `loc` and `scale` have different `dtype`.
        """
        tau, sd = get_tau_sd(tau=tau, sd=sd)
        self.sd = sd
        self.tau = tau

        self.mean = self.median = self.mode = self.mu = mu
        self.variance = 1. / self.tau

        super(Normal, self).__init__(**kwargs)

    def random(self, size=None):
        standard_normal = torch.randn(size)
        sd = self.sd
        mu = self.mu
        
        muvec = mu.repeat(size)
        sdvec = sd.repeat(size)
        
        return standard_normal.mul_(sdvec).add_(muvec)

    def logp(self, value):
        tau = self.tau
        mu = self.mu
        
        muvec = mu.repeat(value.size(0), 1)
        tauvec = tau.repeat(value.size(0), 1)
        
        return (-tauvec * (value - muvec)**2 + torch.log(tauvec / np.pi / 2.)) / 2.

I think it is not easy to implement a distrib.estimation(tensor x, distribution A) function, especially if the distribution contains multiple parameters - most of the case you need to formulate the backward gradient carefully to make it work.
It would be lovely to have a distribution class, even with just simple random number generation and a log-likelihood function. But currently, there are also some crucial functions missing 31 for fast implementation.