The hamiltonian for the whole system can be given by:
$$H_{total}=\sum_{i=1}^{N-2}\frac{p_i^2}{2m}+\sum_{j=1}^{N-1}\frac{k(\gamma-x_j)^2}{2}$$
Where the terms for momentum come from the masses in the chain and the potential comes from the springs. The $\gamma-x$ term comes from the deviation of each spring from their equilibrium position, with $x=\gamma$ giving the point with $0$ potential.
The probability of the system (in thermal contact with surroundings at temperature $T$) being at energy $E$ is given by:
$$q(E)=\frac{1}{Z}e^{-\beta H}$$
Where $\beta=\frac{1}{k_B T}$. The partition function $Z$ is given by integration over phase space of the total hamiltonian of the system. Luckily this hamiltonian can be factorized quite easily.
$$Z=\int_{p,x}e^{-\beta\sum_{1}^{N-2}\frac{p^2}{2m}}e^{-\beta\sum_{1}^{N-1}\frac{k(\gamma-x)^2}{2}}dp\ dx=\int_{-\infty}^{\infty}e^{\frac{\beta(2-N)}{2m}p^2}dp\int_{0}^{\infty}e^{\frac{\beta(1-N)k}{2}(\gamma-x)^2}dx$$
The first is a gaussian ($\int_{-\infty}^{\infty} e^{-ax^2}=\sqrt{\frac{\pi}{a}}$), and the second needs a little massaging.
After integrating out the momentum we get:
$$Z=\sqrt{\frac{2\pi m}{\beta(N-2)}}\int_0^{\infty}e^{-\frac{\beta(N-1)k}{2}(x-\gamma)^2}dx$$
Change the variable $x-\gamma$ to $q$, we get $dx=dq$ and the limits are $\int_{-\gamma}^{\infty}$. This second integral needs the error function to calculate, because of the non-zero lower limit. Change the constants in the exponential into an easier to handle form - $\frac{\beta(N-1)k}{2}=\alpha$:
$$Z=\sqrt{\frac{2\pi m}{\beta(N-2)}}\int_{-\gamma}^{\infty}e^{-\alpha q^2}dq=\sqrt{\frac{2\pi m}{\beta(N-2)}}\sqrt{\frac{\pi}{4\alpha}}\left(1-\text{erf}(-\gamma\sqrt\alpha)\right)$$
Finally:
$$Z=\frac{\pi}{\beta}\sqrt{\frac{m}{k(N-1)(N-2)}}\left(1-\text{erf}(-\gamma\sqrt{\frac{\beta(N-1)k}{2}})\right)$$
There are some approximations that can be made here. If $N$ is large, then we have $(N-1)(N-2)\approx N^2$ and $N-1\approx N$:
$$Z_{large\ N}\approx \frac{\pi}{\beta N}\sqrt{\frac{m}{k}}\left(1-\text{erf}(-\gamma\sqrt{\frac{\beta N k}{2}})\right)$$
Using an approximation for the error function from wikipedia, we can get $Z$ into analytic functions. Using $x=-\gamma\sqrt{\frac{\beta N k}{2}}$, and assuming $\gamma$ is always positive, we get:
$$Z=\frac{\pi}{\beta N}\sqrt{\frac{m}{k}}(1-\sqrt{1-\exp(-x^2\frac{\frac{4}{\pi}+ax^2}{1+ax^2})})$$
Where $a=\frac{8(\pi-3)}{3\pi(4-\pi)}$.
