Why should spin "degeneracy" of electron not be infinite? When we calculate the density of states for an electron, in the standard way as done in statistical mechanics textbooks ( by integrating once over space, then integrating over $\theta$ and $\phi$ of the k-space,etc), we finally multiply a '2' as the spin degeneracy. 
My question is, since the electron can be in an infinite number of states on the Bloch sphere, why are we just taking two states? I know that only two of the infinite number of states are orthogonal, but even if the huge number of states are non-orthogonal, why should not they contribute to an increase in the number of states?
I did not want to use the term degeneracy in the title, because degeneracy is defined as the dimension of the eigenspace of the degenerate eigenvalue, so that is of course two. What I want to know is, why are we taking only the orthogonal states?
 A: 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!)
A: 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!
