This paper presents a probabilistic perspective on iterative methods for approximating the solution $\mathbf{x} \in \mathbb{R}^d$ of a nonsingular linear system $\mathbf{A} \mathbf{x} = \mathbf{b}$. Classically, an iterative method produces a sequence $\mathbf{x}_m$ of approximations that converge to $\mathbf{x}$ in $\mathbb{R}^d$. Our approach, instead, lifts a standard iterative method to act on the set of probability distributions, $\mathcal{P}(\mathbb{R}^d)$, outputting a sequence of probability distributions $\mu_m \in \mathcal{P}(\mathbb{R}^d)$. The output of a probabilistic iterative method can provide both a "best guess" for $\mathbf{x}$, for example by taking the mean of $\mu_m$, and also probabilistic uncertainty quantification for the value of $\mathbf{x}$ when it has not been exactly determined. A comprehensive theoretical treatment is presented in the case of a stationary linear iterative method, where we characterise both the rate of contraction of $\mu_m$ to an atomic measure on $\mathbf{x}$ and the nature of the uncertainty quantification being provided. We conclude with an empirical illustration that highlights the potential for probabilistic iterative methods to provide insight into solution uncertainty.