Distribution Implementations

solidor: 

apaszke: We didn’t plan on adding that, but it seems like a useful thing to have. It’s not going to be a priority for us, but if someone wants to send a PR, then we’ll be happy to merge it in.

junpenglao: I would love to contribute as well. However, I am still getting familiar with the infrastructure of PyTorch. It would be great if someone from the dev could provide a template for implementation (e.g., a Normal distribution), ideally one for the continuous and one for the discrete distribution, as a reference for best practice using PyTorch.

apaszke: 

junpenglao: 

apaszke: <SCODE>def normal_pdf(x, mean, std):
    # compute the pdf here using x, mean and std Variables
<ECODE>

junpenglao: 

apaszke: How do you compute the gradients in that case? Right now StochasticFunctions will use REINFORCE.

alexis-jacq: A distribution could be implemented as a class containing one Function “self.pdf” and a method sample() <SCODE>class Distribution():
    def __init__(self):
        self.pdf = # autograd Function

    def sample():
        # use self.pdf for Monte Carlo algorithm
<ECODE> ( by the way, how do you use python coloration when providing code in this forum ?)

apaszke: When opening up the block append “python” right after the opening backticks:

alexis-jacq: 

alexis-jacq: <SCODE>W = torch.rand(nv, nh)); a = torch.rand(nh)); b = torch.rand(nv));
net_forward = nn.Sequential(nn.Linear(wieght=W, bias=a), nn.Sigmoid())
net_backward = nn.Sequential(nn.Linear(wieght=W, bias=b), nn.SoftMax())

# training:
for input in training_set:
    v1 = input # from training set

    x = net_forward(v1)
    H1 = distrib.estimation(x, distrib.Bernouilli())
    h1 = distrib.sample(H1)

    y = net_backward(h)
    V2 = distrib.estimation(y, distrib.Multinomial())
    v2 = distrib.sample(V2)

    x = net_forward(v2)
    H2 =distrib. estimation(x, distrib.Bernouilli())
    h2 = distrib.sample(H2)

    # gradient descent step:
    W -= lr * (v1.h1 - v2.h2)
    a -= lr * (h1 - h2)
    b -= lr * (v1 - v2)
<ECODE>

junpenglao: <SCODE>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.
<ECODE>

chenyuntc: I think distribution could be implemented as a collection of functions as below: <SCODE>class Distribution:

    @staticmethod
    def pdf(self):
        raise NotImplementedError

    @staticmethod
    def logpdf(self):
        raise NotImplementedError
    
    @staticmethod
    def cdf(self):
        raise NotImplementedError
    
    @staticmethod
    def qtl(self):
        raise NotImplementedError

    @staticmethod
    def rnd(self):
        raise NotImplementedError

class Gauss(Distribution):
    #scalar version
    @staticmethod
    def logpdf(x, mu, sigma):
        return -0.5 * (x-mu) * (x-mu)/(sigma*sigma) - 0.5 * torch.log(2*math.pi*sigma*sigma)

<ECODE> and it could be used like this: <SCODE>Gauss.logpdf(x, mu, sigma)
<ECODE>

kirk86: 

chenyuntc: 

kirk86: 

stepelu: 

Andres_Masegosa: Any advances on that? Highly interested.

chssch: