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 [1,2]. 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 (GES) with a statistical parameter $alpha =1-2/gamma$, i.e. the quasimomenta of $N$ strongly interacting bosons map to the momenta of $N$ free fermions via $k_iapprox alpha k_i^F $ with $i=1,ldots, N$. Here $gamma$ is the dimensionless interaction strength within the Lieb-Liniger model. We rigorously prove that such a statistical parameter $alpha$ 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 $gamma$ and $alpha$ 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 [3-5]). We also express the leading order of the short distance emph{non-local} correlation functions $langlePsi^dagger(x_1)cdotsPsi^dagger(x_M)Psi(y_M)cdotsPsi(y_1)rangle$ of the strongly repulsive Bose gas in terms of the wave function of $M$ bosons at zero collision energy and zero total momentum. Here $Psi(x)$ 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.