Why should spin "degeneracy" of electron not be infinite?

Good question! Really what you care about when calculating the density of states is the number of orthogonal states, whether you are talking about spin, spatial states, or whatever. In other words, you want to know density of Hilbert space dimension. To give a non-spin example if you are looking at the density of states in an $L\times L\times L$ box of Fermi gas, you get the momentum states by enforcing periodic boundary conditions... but by your argument, as soon as you have two momentum eigenfunctions, you should have infinite density of states, because each single particle could be in any superposition of those two states. As you know, this is not actually the case: the same applies for spins. Count the number of orthogonal states, not the total number of vectors in the Hilbert space (usually infinite!)

I like Will's answer, but I thought I would also add one showing where the factor of $2$ (and, in general, the dimension of the degenerate Hilbert space) actually comes from in the formalism.

In the canonical formulation, if we decompose the Hilbert space in question as $\mathcal{H}=\mathcal{H}_1\otimes\mathcal{H}_2$, where you can think of $\mathcal{H}_1$ as describing the space component and $\mathcal{H}_2$ as describing spin, then the partition function for the system at temperature $\beta=1/kT$ is

$$Z\equiv\sum_{E}e^{-\beta E}=\text{Tr}\left(e^{-\beta H}\right),$$

where $H$ is the Hamiltonian acting on $\mathcal{H}$, and the $\text{Tr}$ trace is acting on the whole Hilbert space as well.

Now, if all spins are degenerate at each energy level, then the Hamiltonian naturally decomposes as $H=H_{s}\otimes\textbf{1}$, where $H_{s}$ can be thought of as the spatial Hamiltonian. Then the partition function will be

$$Z=\text{Tr}\left(e^{-\beta H_s\otimes\textbf{1}}\right)=\text{Tr}\left(e^{-\beta H_s}\otimes\textbf{1}\right)=\text{tr}_1\left(e^{-\beta H_s}\right)\text{tr}_2(\textbf{1})=\text{dim}(\mathcal{H}_2)Z_{s},$$

where $Z_{s}$ is defined to be the partition function for $H_s$. This is schematically where the dimension of the degenerate Hilbert space comes from in these calculations. In particular, since the density of states for a continuous spectrum is defined so that

$$Z=\int\mathrm{d}E\,g(E)\,e^{-\beta E},$$

it is clear that having a degenerate subspace $\mathcal{H}_2$ corresponds to multiplying the density of states by $\dim(\mathcal{H}_2)$, or, in the case of spin, multiplying by $2$ (or $2s+1$ for a massive spin-$s$ particle).

I hope this helps!