We observe $n$ pairs of independent (but not necessarily i.i.d.) random variables $X_{1}=(W_{1},Y_{1}),\ldots,X_{n}=(W_{n},Y_{n})$ and tackle the problem of estimating the conditional distributions $Q_{i}^{\star}(w_{i})$ of $Y_{i}$ given $W_{i}=w_{i}$ for all $i\in\{1,\ldots,n\}$. Even though these might not be true, we base our estimator on the assumptions that the data are i.i.d.\ and the conditional distributions of $Y_{i}$ given $W_{i}=w_{i}$ belong to a one parameter exponential family $\bar{\mathscr{Q}}$ with parameter space given by an interval $I$. More precisely, we pretend that these conditional distributions take the form $Q_{{\boldsymbolθ}(w_{i})}\in \bar{\mathscr{Q}}$ for some ${\boldsymbolθ}$ that belongs to a VC-class $\bar{\boldsymbolΘ}$ of functions with values in $I$. For each $i\in\{1,\ldots,n\}$, we estimate $Q_{i}^{\star}(w_{i})$ by a distribution of the same form, i.e.\ $Q_{\hat{\boldsymbolθ}(w_{i})}\in \bar{\mathscr{Q}}$, where $\hat {\boldsymbolθ}=\hat {\boldsymbolθ}(X_{1},\ldots,X_{n})$ is a well-chosen estimator with values in $\bar{\boldsymbolΘ}$. We show that our estimation strategy is robust to model misspecification, contamination and the presence of outliers. Besides, we provide an algorithm for calculating $\hat{\boldsymbolθ}$ when $\bar{\boldsymbolΘ}$ is a VC-class of functions of low or moderate dimension and we carry out a simulation study to compare the performance of $\hat{\boldsymbolθ}$ to that of the MLE and median-based estimators.