Numerical simulations of realistic non-vacuum magnetospheres of isolated neutron stars have shown that pulsar spin-down luminosities depend weakly on the magnetic obliquity α. In particular, L ∝ B 2 (1 + sin 2 α), where B is the magnetic field strength at the star surface. Being the most accurate expression to date, this result provides the opportunity to estimate B for a given radiopulsar with quite a high accuracy. In the current work, we present a refinement of the classical 'magneto-dipolar' formula for pulsar magnetic fields B md = (3.2 × 10 19 G) PṖ , where P is the neutron star spin period. The new, robust timing-based estimator is introduced as log B = log B md + ∆ B (M, α), where the correction ∆ B depends on the equation of state (EOS) of dense matter, the individual pulsar obliquity α and the mass M . Adopting state-of-theart statistics for M and α we calculate the distributions of ∆ B for a representative subset of 22 EOSs that do not contradict observations. It has been found that ∆ B is distributed nearly normally, with the average in the range −0.5 to −0.25 dex and standard deviation σ[∆ B ] ≈ 0.06 to 0.09 dex, depending on the adopted EOS. The latter quantity represents a formal uncertainty of the corrected estimation of log B because ∆ B is weakly correlated with log B md . At the same time, if it is assumed that every considered EOS has the same chance of occurring in nature, then another, more generalized, estimator B * ≈ 3B md /7 can be introduced providing an unbiased value of the pulsar surface magnetic field with ∼30 per cent uncertainty with 68 per cent confidence. Finally, we discuss the possible impact of pulsar timing irregularities on the timing-based estimation of B and review the astrophysical applications of the obtained results.