This function accomodates missing elements for the data of simple
    Stochastic distribution subclasses. The masked_values argument is an
    object of type numpy.ma.MaskedArray, which contains the raw data and
    a boolean mask indicating missing values. The resulting list contains
    a list of stochastics of type dist_class, with the extant values as data
    stochastics and the missing values as variable stochastics.

    :Arguments:
      - name : string
        Name of the data stochastic
      - dist_class : Stochastic
        Stochastic subclass such as Poisson, Normal, etc.
      - imputable : numpy.ma.core.MaskedArray or iterable
        A masked array with missing elements (where mask=True, value
        is assumed missing), or any iterable that contains None
        elements that will be imputed.
      - parents (optional): dict
        Arbitrary keyword arguments.

def Impute(name, dist_class, imputable, **parents):
    """
    This function accomodates missing elements for the data of simple
    Stochastic distribution subclasses. The masked_values argument is an
    object of type numpy.ma.MaskedArray, which contains the raw data and
    a boolean mask indicating missing values. The resulting list contains
    a list of stochastics of type dist_class, with the extant values as data
    stochastics and the missing values as variable stochastics.

    :Arguments:
      - name : string
        Name of the data stochastic
      - dist_class : Stochastic
        Stochastic subclass such as Poisson, Normal, etc.
      - imputable : numpy.ma.core.MaskedArray or iterable
        A masked array with missing elements (where mask=True, value
        is assumed missing), or any iterable that contains None
        elements that will be imputed.
      - parents (optional): dict
        Arbitrary keyword arguments.
    """

    dims = np.shape(imputable)
    masked_values = np.ravel(imputable)

    if not isinstance(masked_values, np.ma.core.MaskedArray):
        # Generate mask

        mask = [v is None or np.isnan(v) for v in masked_values]
        # Generate masked array
        masked_values = np.ma.masked_array(masked_values, mask)

    # Initialise list
    vars = []
    for i in xrange(len(masked_values)):

        # Name of element
        this_name = name + '[%i]' % i
        # Dictionary to hold parents
        these_parents = {}
        # Parse parents
        for key, parent in six.iteritems(parents):

            try:
                # If parent is a PyMCObject
                shape = np.shape(parent.value)
            except AttributeError:
                shape = np.shape(parent)

            if shape == dims:
                these_parents[key] = Lambda(key + '[%i]' % i,
                                            lambda p=np.ravel(parent),
                                            i=i: p[i])
            elif shape == np.shape(masked_values):
                these_parents[key] = Lambda(key + '[%i]' % i, lambda p=parent,
                                            i=i: p[i])
            else:
                these_parents[key] = parent

        if masked_values.mask[i]:
            # Missing values
            vars.append(dist_class(this_name, **these_parents))
        else:
            # Observed values
            vars.append(dist_class(this_name, value=masked_values[i],
                                   observed=True, **these_parents))
    return np.reshape(vars, dims)