We introduce a robust fitting method into pulsar timing analysis to cope with the non-Gaussian noise. The general maximum likelihood estimator (M-estimator) can resist the impact of non-Gaussian noise by employing convex and bounded loss functions. Three loss functions, including the Huber function, the Bisquare function and the Welsch function, are investigated. The Shapiro–Wilk test is employed to test whether the uncertainty in the observed times of arrival is drawn from a non-Gaussian distribution. Two simulations, where the non-Gaussian distribution is modelled as contaminated Gaussian distributions, are performed. It is found that M-estimators are unbiased and could achieve a root-mean-square error smaller than that obtained by the least square (LS) at the cost of a slightly higher computation complexity in a non-Gaussian environment. M-estimators are also applied to the real timing data of PSR J1713+0747. The results have shown that the fitting results of M-estimators are more accurate than those of LS and are closer to the result of very long baseline interferometry.