There are several long-term statistical researches using the Molchan diagram (MD) to prove the relation between thermal infrared (TIR) anomalies and earthquakes in different regions, however, these studies are flawed: 1) the original MD is based on the spatially uniform Poisson model, but it will offer wrong evaluations for inhomogenous systems with uneven spatial-temporal distribution of earthquakes; 2) some of these studies have de-clustered earthquake catalogs before applying MD, however, de-clustering will introduce ambiguities to final scores and it may also conceal the potential relation between precursors and the ‘removed’ events. Until now, we have to admit that the long-term statistical evidence for proving the correspondence between earthquake and TIR anomalies is still absent and the power of TIR anomalies for predicting earthquakes is limited. In this study, we use daily nighttime Outgoing Longwave Radiation (OLR) data provided by the National Oceanic and Atmospheric Administration (NOAA) to extract the TIR anomalies of Chinese Mainland (20°-54°N, 73°-135°E). The data from Jan. 2007 to Dec. 2010 is the training dataset to obtain the best parameters for extracting TIR anomalies and the best time-distance-magnitude (TDM) windows for determining the correspondence between earthquakes and TIR anomalies, and the data from Jan. 2011 to Dec. 2017 is the testing dataset. The new 3D Molchan diagram offers a score for each model with different parameters. Unlike the original MD that only deals with the rate of missed events and the size of warning space, the 3D Molchan diagram quantifies the errors including false alarms and missed predictions. We assume that the best parameters and TDM windows are spatially variable for different sub-regions, because the Signal/Noise ratio is spatially variable due to the different geological and meteorological backgrounds. Moreover, we construct a new probability prediction model based on non-seismic binary alarms. Results show that the TIR anomalies is strongly related to normal or reverse earthquakes with magnitude≥ 4.0, and the TIR anomalies caused by earthquakes should be persistent in space and time. Moreover, the spatially variable model is superior to the global invariant one. We succeed in transforming the non-seismic binary alarms into probabilistic predictions based on the TIR anomalies and Relative Intensity index, which is defined as the rate of past earthquakes occurring in each spatial cell. Moreover, our new probabilistic model is superior to the spatially inhomogeneous Poisson model. However, this new probability model is still naïve and weak, and needs to be improved in the future.