The objective of this work is to quantify the uncertainty in probability of failure estimates resulting from incomplete knowledge of the probability distributions for the input random variables. We propose a framework that couples the widely used Subset simulation (SuS) with Bayesian/information theoretic multi-model inference. The process starts with data used to infer probability distributions for the model inputs. Often such data sets are small. Multi-model inference is used to assess uncertainty associated with the model-form and parameters of these random variables in the form of model probabilities and the associated joint parameter probability densities. A sampling procedure is used to construct a set of equally probable candidate probability distributions and an optimal importance sampling distribution is determined analytically from this set. Subset simulation is then performed using this optimal sampling density and the resulting conditional probabilities are re-weighted using importance sampling. The result of this process are empirical probability distributions of failure probabilities that provide direct estimates of the uncertainty in failure probability estimates that result from inference on small data sets. The method is demonstrated to be both computationally efficient -- requiring only a single subset simulation and nominal cost of sample re-weighting -- and to provide reasonable estimates of the uncertainty in failure probabilities.