1. A computer-implemented method for multiple imputation for retail data sets with missing data values, the method comprising: receiving an original data set including values including a plurality of products, a plurality of stores or chains in which each said product is sold, and a plurality of time-periods indicating when said products were sold; identifying and encoding the missing data values in the original data set with dummy indicator variables corresponding to specific product, store and time-period combinations; obtaining a joint probability distribution for the magnitudes of the missing data values in the original data set, the obtaining the joint probability distribution comprising: specifying a probability model for the entries of the original data set based on a mean value obtained from a tensor-product factorization of dimensions comprising of product, store and time-period, and additionally, comprised of an additive noise term that has a zero mean and non-zero variance, and for obtaining a likelihood function for non-missing values of the original data set based on the probability model; specifying probability models with parameters for latent factors in this tensor-product factorization; specifying a posterior joint conditional distribution for said latent factors, the parameters in the probability models for these latent factors, and the said non-zero variance of the additive noise term, given the non-missing data values in the original data set; and specifying the joint distribution of the missing values in the original data set, based on marginalizing the likelihood function over the known non-missing values, given said posterior joint conditional distribution; generating a plurality of complete data sets corresponding to the original data set, wherein each complete data set in said plurality of complete data sets corresponds to the original data set with its non-missing values intact, and replacing, in each of the complete data sets, missing values indicated by said dummy variables with a sampled set of values from the joint probability distribution for the magnitudes of the missing elements as obtained, wherein a programmed processor device performs one or more of one or more the receiving, identifying and encoding, obtaining, generating and replacing.