The correlation function is an important quantity in the physics of ultracold quantum gases because it provides information about the quantum many-body wave function beyond the simple density profile. In this paper we first study the M-body local correlation functions, g M , of the one-dimensional (1D) strongly repulsive Bose gas within the Lieb–Liniger model using the analytical method proposed by Gangardt and Shlyapnikov (2003 Phys. Rev. Lett. 90 010401; 2003 New J. Phys. 5 79). In the strong repulsion regime the 1D Bose gas at low temperatures is equivalent to a gas of ideal particles obeying the non-mutual generalized exclusion statistics with a statistical parameter , i.e. the quasimomenta of N strongly interacting bosons map to the momenta of N free fermions via with . Here γ is the dimensionless interaction strength within the Lieb–Liniger model. We rigorously prove that such a statistical parameter α solely determines the sub-leading order contribution to the M-body local correlation function of the gas at strong but finite interaction strengths. We explicitly calculate the correlation functions g M in terms of γ and α at zero, low, and intermediate temperatures. For M = 2 and 3 our results reproduce the known expressions for g 2 and g 3 with sub-leading terms (see for instance (Vadim et al 2006 Phys. Rev. A 73 051604(R); Kormos et al 2009 Phys. Rev. Lett. 103 210404; Wang et al 2013 Phys. Rev. A 87 043634). We also express the leading order of the short distance non-local correlation functions of the strongly repulsive Bose gas in terms of the wave function of M bosons at zero collision energy and zero total momentum. Here is the boson annihilation operator. These general formulas of the higher-order local and non-local correlation functions of the 1D Bose gas provide new insights into the many-body physics.