The outermost strong part of a planet is called the lithosphere. When loads, such as volcanoes, sediments, or intrusions, are applied to the lithosphere, it flexes. The amount of flexure is controlled by the flexural rigidity of the lithosphere. The elastic thickness Te is the thickness of an equivalent fully elastic spherical shell which flexes in the same way as the real lithosphere. It is an important quantity because it strongly depends on the lithospheric heat flux at the time that the load was applied, which is controlled by the thermal evolution of the planet and varies in space and time. The loads and the associated lithospheric deflections cause gravity anomalies and topographic relief. Observations of these can be used to constrain Te. In this study, a global map of the elastic thickness of Mars is presented. Several recent missions to Mars have provided global spherical harmonics data sets of gravity and topography. These data are inverted for the shape of the crust-mantle boundary, or equivalently the crustal thickness. This is assuming that all gravity anomalies are caused by only two density interfaces which are the surface and the crust-mantle boundary. The large amplitudes of the Martian gravity field necessitate the application of a finite amplitude correction. A simple, but realistic model, which allows loading and compensation at the same two interfaces, is derived. It uses the differential equations for the flexure of a thin elastic shell and depends on six parameters: the elastic thickness Te, the ratio F of the amplitude of the loads at the two interfaces before flexure, the correlation r between these loads, and three parameters of a covariance function of the isotropic Matérn class which describes the topography before flexure. The input data are localized to specific grid points of a map using multitaper spectral estimation. Contrary to most elastic thickness studies which compute observed and modelled admittance or coherence to find a best-fit solution for Te, this study uses maximum likelihood estimation as first proposed by Simons and Olhede (2013). This technique allows to determine the parameter set which is most likely to have produced the localized estimates of the topography and the shape of the crust-mantle boundary. Maps of the most likely parameter sets are presented for different localization window sizes. The results generally agree with previous studies, yielding Te ~10km in the southern uplands and higher values at the large volcanoes. This also corresponds to thermal evolution models predicting a more rigid lithosphere in more recently formed areas. Log-likelihood contours and Monte Carlo simulations with synthetically generated topographies reveal the quality of the results. The elastic thickness is well constrained in the southern uplands and at Elysium and Ascraeus Mons, but poorly constrained in the northern lowlands and at the other volcanoes. While this study shows that it is possible to retrieve Te with maximum likelihood estimation, more research is needed to explain these poor constraints