Why are anticommutators needed in quantization of Dirac fields? Why is the anticommutator actually needed in the canonical quantization of free Dirac field?
 A: Assuming the question is about why anticommutators rather than commutators,as per David's comment:
If I denote the Dirac particle properties by a general index $\alpha$ - these can include spins, momenta - then if I create a two particle state, with particle 1 having $\alpha_1$ and particle 2 $\alpha_2$ then the state is given by applying creation operators in the order
$$|\alpha_{1} \alpha_{2}\rangle = b_{\alpha_1}^{\dagger}b_{\alpha_2}^{\dagger}|0\rangle>$$  If we create them the other way round:$$|\alpha_{2} \alpha_{1}\rangle = b_{\alpha_2}^{\dagger}b_{\alpha_1}^{\dagger}|0\rangle>$$ then if the b's obey anticommutation relations, it is easy to see that $$|\alpha_{2} \alpha_{1}\rangle = -|\alpha_1\alpha_2\rangle$$ i.e. the theory naturally reproduces fermi statistics, as you would want for spin 1/2 Dirac particles.
A: The most elementary reason is that the Dirac field Hamiltonian is bounded below only when you use anticommutation relations on the creation/annihilation operators instead of commutators. A free quantum field theory with energy unbounded below has no stable vacuum.
It is easiest to demonstrate this in two dimensions, where there are no polarization issues.
Instructive 2d example
In two dimensions (one space one time), there is a nice dimensionally reduced analog, which is the right-moving (necessarily massless) Majorana-Weyl Fermion (the argument also works with 2d Dirac fermions with two components, but this is the simplest case). This is a single component field $\psi$ which obeys the equation of motion
$$ (\partial_t -\partial_x) \psi = 0 $$
This simple equation is derived from the 2d Dirac equation using the (real convention, explicitly real) 2d Dirac matrices (0,1;-1,0) an (0,1;1,0), which are $\gamma_0 = \sigma_x$ and $\gamma_1 = i\sigma_y$. They square to 1 and -1 respectively, and they anticommute, so they reproduce the 1+1 dimensional metric tensor. The $\gamma_5$ analog, which I'll call $\Gamma$ to accommodate different dimensions, is diagonal in this explicit representation, and $\Gamma=\sigma_z$.
The two eigenvectors of $\Gamma$ propagate independently by the 2d massless equation of motion
$$ \gamma_i \partial_i \psi = 0 $$
And further, because the $\gamma$ matrices are real, this is a Majorana representation (most physicists write the dirac equation with an i factor in front of the derivative, so that the Dirac matrices for a Majorana representation are purely imaginary. I'm using a mathematician's convention for this, because I like the equations of motion to be real. Others like the k-space propagator to not have factors of i in the k part. Unfortunately, physicists never settled on a unique sensible convention--- everyone has their own preferred way to write Dirac matrices). So it is sensible in the equation of motion to restrict $\psi$ to be Hermitian, since its Hermitian conjugate obeys the exact same equation.
So that the field has a k decomposition
$$ \psi(x) = \int a_k e^{ikx - ikt }dk$$
And the reality condition (Hermiticity) tells you that $a^{\dagger}(-k) = a(k)$ (one should say that the normalization of the $a$ operators expansion is not completely conceptually trivial--- the $a$'s are both relativistically and nonrelativistically normalized, because the spinor polarization $\sqrt{w}$ factor cancels the mass-shell hyperbola factor, so that the dk integration is not weighted by anything, it's just the normal calculus integral with uniform measure)
An operator with definite frequency, which (Heisenberg picture) evolves in time according to
$$ \partial_t O = i\omega O$$
Has the property that it is a raising operator--- acting with this operator adds $\omega$ to the energy. If $\omega$ is negative, $a$ is an annihilation operator. The condition that the vacuum is stable says that all the annihilation operators give 0 when acting on the vacuum state.
But notice that the frequency in the expansion of $\psi$ changes sign at $k=0$. This came from the linearity of the Dirac Hamiltonian in the momenta. It means that the operator $a_k$ acts to raise the energy for k>0, but acts to lower the energy for $k<0$. This means that the $k>0$ operators create, and the $k<0$ operators annihilate, so that the right way to $a^{\dagger}(-k)$ are creation operators, while the $k<0$ operators are annihilation operators.
The energy operator counts the number of particles of momentum k, and multiplies by their energy:
$$ H = \int_{k>0} k a^{\dagger}(k) a(k) dk $$
And this is manifestly not a local operator, it is defined only integrated over k>0. To make it a local operator, you need to extend the integration to all k, but then the negative k and positive k contributions have opposite sign, and they need to be equal. To arrange this, you must take anticommutation relations
$$ \{ a^{\dagger}(k),a(k)\} = i\delta(k-k') $$
And then
$$ H = {1\over 2} \int k a^{\dagger}(k) a(k) = \int \psi(x) i \partial_x \psi(x) dx $$
Note that this looks like it is a perfect derivative, and it would be if $\psi$ weren't anticommuting quantity. For anticommuting quantities, 
$$ \partial_x \psi^2 = \psi \partial_x \psi + \partial_x \psi \psi $$
Which is zero, because of the anticommutation.
Deeper reasons
Although this looks like an accidental property, that the energy was negative without anticommutators, it is not. The deeper reason is explained with Euclidean field theory using a Feynman-Schwinger formalism, but this requires understanding of the Euclidean and path integral versions of anticommuting fields, which requires being comfortable with anticommuting quantities, which requires a motivation. So it is best to learn the shallow reason first.
