  A probabilistic approach is developed for the exact solution $u$ to a
determinist partial differential equation as well as for its associated
approximation $u^{(k)}_{h}$ performed by $P_k$ Lagrange finite element. Two
limitations motivated our approach: on the one hand, the inability to determine
the exact solution $u$ to a given partial differential equation (which
initially motivates one to approximating it) and, on the other hand, the
existence of uncertainties associated with the numerical approximation
$u^{(k)}_{h}$. We thus fill this knowledge gap by considering the exact
solution $u$ together with its corresponding approximation $u^{(k)}_{h}$ as
random variables. By way of consequence, any function where $u$ and
$u_{h}^{(k)}$ are involved as well. In this paper, we focus our analysis to a
variational formulation defined on $W^{m,p}$ Sobolev spaces and the
corresponding a priori estimates of the exact solution $u$ and its
approximation $u^{(k)}_{h}$ to consider their respective $W^{m,p}-$norm as a
random variable, as well as the $W^{m,p}$ approximation error with regards to
$P_k$ finite elements. This will enable us to derive a new probability
distribution to evaluate the relative accuracy between two Lagrange finite
elements $P_{k_1}$ and $P_{k_2}, (k_1 < k_2)$.
