We study mixture of linear regression (random coefficient) models, which capture population heterogeneity by allowing the regression coefficients to follow an unknown distribution $G^*$. In contrast to common parametric methods that fix the mixing distribution form and rely on the EM algorithm, we develop a fully nonparametric maximum likelihood estimator (NPMLE). We show that this estimator exists under broad conditions and can be computed via a discrete approximation procedure inspired by the exemplar method. We further establish theoretical guarantees demonstrating that the NPMLE achieves near-parametric rates in estimating the conditional density of $Y|X$, both for fixed and random designs, when $σ$ is known and $G^*$ has compact support. In the random design setting, we also prove consistency of the estimated mixing distribution in the Lévy-Prokhorov distance. Numerical experiments indicate that our approach performs well and additionally enables posterior-based individualized coefficient inference through an empirical Bayes framework.