Let us consider Dirac equation $$(i\gamma^\mu\partial_\mu -m)\psi ~=~0$$ as a classical field equation. Is it possible to introduce Poisson bracket on the space of spinors $\psi$ in such a way that Dirac equation becomes Hamiltonian equation
$$\dot{\psi}~=~\{ \psi,H\}_{PB}~?$$
Of course, such Poisson bracket would be graded (super Poisson bracket), but if it exists this would explain on classical level why $\frac{1}{2}$-spinors correspond to fermions.
Given the Dirac Lagrangian density
$$ \begin{align} {\cal L}~=~&\overline{\psi}(i\sum_{\mu=0}^3\gamma^{\mu}\partial_{\mu}-m)\psi, \cr \overline{\psi}~:=~&\psi^{\dagger}\gamma^0, \cr \{\gamma^{\mu},\gamma^{\nu}\}_{+} ~=~&2\eta^{\mu\nu}{\bf 1}_{4\times 4}, \end{align} \tag{1} $$
with Minkowski signature $(+,-,-,-)$, and $\psi$ is a Grassmann-odd Dirac-spinor, the question is How to find the corresponding Hamiltonian formalism?
The Legendre transformation of (1) is singular. The Dirac-Bergmann analysis of the theory (1) leads to constraints, cf. e.g. Ref. 1 or this Phys.SE post. Here we will instead take a shortcut using the Faddeev-Jackiw method.
I) Complex Grassmann-fields. We first identify the Hamiltonian density ${\cal H}$ as (minus) the terms in (1) that don't involve time derivatives:
$$ \begin{align} {\cal L}~=~&i\psi^{\dagger}\dot{\psi}-{\cal H}, \cr {\cal H}~=~& \overline{\psi} (-i\sum_{j=1}^3\gamma^{j}\partial_{j}+m)\psi. \end{align}\tag{2}$$
The symplectic one-form potential can be transcribed from the kinetic term in (2):
$$ \vartheta(t) ~=~\int\! d^3x~ i\psi^{\dagger}({\bf x},t) ~\mathrm{d}\psi({\bf x},t), \tag{3} $$
where $\mathrm{d}$ denotes the exterior derivative$^1$ on the infinite-dimensional configuration space for the fermion field. The symplectic two-form is then
$$ \begin{align} \omega(t)~=~&\mathrm{d}\vartheta(t) ~=~\int\! d^3x~ i\mathrm{d}\psi^{\dagger}({\bf x},t) \wedge \mathrm{d}\psi({\bf x},t) \cr ~=~&\int\! d^3x~d^3y~ i\mathrm{d}\psi^{\dagger}({\bf x},t) \wedge \delta^3({\bf x}-{\bf y}) ~\mathrm{d}\psi({\bf y},t). \end{align}\tag{4} $$
The equal-time super-Poisson/Dirac bracket on fundamental fields is the inverse supermatrix of the supermatrix for the symplectic two-form (4):
$$ \begin{align} \{\psi_{\alpha}({\bf x},t), \psi^{\dagger}_{\beta}({\bf y},t)\}_{PB}~=~& -i \delta_{\alpha\beta}~\delta^3({\bf x}-{\bf y})\cr ~=~&\{\psi^{\dagger}_{\alpha}({\bf x},t), \psi_{\beta}({\bf y},t)\}_{PB}, \end{align}\tag{5} $$
and other fundamental super-Poisson brackets vanish, cf. e.g. my Phys.SE answer here. Due to the QM correspondence principle, the canonical anticommutation relations (CARs) are the super-Poisson brackets (5) multiplied with $i\hbar$:
$$ \begin{align} \{\hat{\psi}_{\alpha}({\bf x},t), \hat{\psi}^{\dagger}_{\beta}({\bf y},t)\}_{+} ~=~& \hbar\delta_{\alpha\beta}~\delta^3({\bf x}-{\bf y})\hat{\bf 1}\cr ~=~&\{\hat{\psi}^{\dagger}_{\alpha}({\bf x},t), \hat{\psi}_{\beta}({\bf y},t)\}_{+}, \end{align}\tag{6} $$ and other CARs vanish.
I) Real Grassmann-fields. Alternatively, let us decompose the complex Dirac spinor
$$\psi_{\alpha}~\equiv~(\psi^1_{\alpha}+i\psi^2_{\alpha})/\sqrt{2} \quad\text{and}\quad \psi^{\dagger}_{\alpha}~\equiv~(\psi^1_{\alpha}-i\psi^2_{\alpha})/\sqrt{2}, \tag{7} $$
in real and imaginary parts. The Lagrangian density (2) reads up to total derivative terms$^2$
$$\begin{align} {\cal L}~=~&\frac{i}{2}\left(\psi^{\dagger}\dot{\psi}- \dot{\psi}^{\dagger}\psi\right)-{\cal H}\cr ~=~&\frac{i}{2}\sum_{a=1}^2(\psi^a)^T\dot{\psi}^a-{\cal H}. \end{align}\tag{2'} $$
The corresponding symplectic one-form potential is
$$ \vartheta(t) ~=~\sum_{a=1}^2\int\! d^3x~ \frac{i}{2}\psi^a({\bf x},t)^T ~\mathrm{d}\psi^a({\bf x},t). \tag{3'}$$
The symplectic two-form is
$$\begin{align} \omega(t)~=~&\mathrm{d}\vartheta(t) ~=~\sum_{a=1}^2\int\! d^3x~ \frac{i}{2}\mathrm{d}\psi^a({\bf x},t)^T \wedge \mathrm{d}\psi^a({\bf x},t) \cr ~=~&\sum_{a,b=1}^2\int\! d^3x~d^3y~ \frac{i}{2}\mathrm{d}\psi^a({\bf x},t)^T \wedge \delta_{ab}~\delta^3({\bf x}-{\bf y}) ~\mathrm{d}\psi^b({\bf y},t). \end{align}\tag{4'}$$
The equal-time super-Poisson is
$$ \{\psi^a_{\alpha}({\bf x},t), \psi^b_{\beta}({\bf y},t)\}_{PB}~=~ -i \delta^{ab}~\delta_{\alpha\beta}~\delta^3({\bf x}-{\bf y}).\tag{5'} $$
The CARs are
$$ \{\hat{\psi}^a_{\alpha}({\bf x},t), \hat{\psi}^b_{\beta}({\bf y},t)\}_{+} ~=~ \hbar\delta^{ab}~\delta_{\alpha\beta}~\delta^3({\bf x}-{\bf y})\hat{\bf 1} . \tag{6'} $$
References:
--
$^1$ In our super-conventions, the exterior derivative $\mathrm{d}$ is Grassmann-even and carries form-degree +1.
$^2$ Note that adding a total time derivative
$$i\psi^{\dagger}\dot{\psi}~\longrightarrow~i\psi^{\dagger}\dot{\psi}+ \frac{d}{dt}(\alpha\psi^{\dagger}\psi)\tag{8} $$
to the kinetic term (2) corresponds to adding an exact term
$$ \vartheta(t)~\longrightarrow~\vartheta(t)+ \mathrm{d} \int\! d^3x~ \alpha\psi^{\dagger}({\bf x},t) \psi({\bf x},t) \tag{9} $$
to the symplectic one-form potential (3), which has no effect on the symplectic 2-form (4).
Now, I don't know what the word rigorous means, but here is a straight off the bat naive answer. Given $$ H = \int d^3 x \, \bar{\Psi}(i \gamma_i \partial_i +m)\Psi $$ from $$ \mathcal{L} = i\bar{\Psi}\gamma^{\mu}\partial_{\mu}\Psi - \bar{\Psi} m \Psi \quad \text{and}\quad H = \int d^3x \, (\pi \dot{\Psi}-\mathcal{L}) $$ with $(+,-,-,-)$. Let's use the possoin bracket $$ \{A,B\} = \int d^3 x \, (\delta_{\Psi}A \, \delta_{\pi}B-\delta_{\pi}A \, \delta_{\Psi}B) $$ and let's remember that $\pi_{\Psi} = i \bar{\Psi}\gamma_0 \implies \bar{\Psi} = -i \pi \gamma_0$ so that $$ H = \int d^3x \, -i\gamma_0\pi(i\gamma_i\partial_i +m)\Psi $$ consider $$ \delta_{\pi}H = \int d^3x \, \gamma_0 (\delta_{\pi}\pi) \gamma_i \partial_i \Psi + \gamma_0\pi\gamma_i\partial_i(\delta_{\pi}\Psi) - i\gamma_0 (\delta_{\pi}\pi) m\Psi - i\gamma_0 \pi m (\delta_{\pi}\Psi) $$ $$ =\gamma_0 \gamma_i \partial_i \Psi - i\gamma_0 m \Psi $$ Then looking at $$ \{H,\Psi\} = -\gamma_0 \gamma_i \partial_i \Psi + i\gamma_0 m \Psi = \dot{\Psi} $$ God knows if this is right, but this sounds like what you are looking for.
In my opinion, one cannot rigorously [*] define the bracket. Suppose that you use the Dirac field equation for arriving to the ordinary Lagrangian density
$$ \mathcal{L} = c \bar{\psi} \left( i\hbar \gamma^\mu \frac{\partial}{\partial x^\mu} \right) \psi $$
This is a function of the spinor components $\psi_i$ and their adjoints $\bar{\psi_i}$. The problem begins when you try to obtain the conjugate momentum for the adjoints (the dot denotes time derivative)
$$\bar{\pi_i} = \frac{\partial \mathcal{L}}{\partial \dot{\bar{\psi_i}}} = 0$$
which implies that not all the canonical variables are independent and not true 'phase-space' structure exists.
You could try to formally define Poisson brackets in the usual fashion,
$$ \{ A, B \} \equiv \sum_i \frac{\partial A}{\partial \psi_i} \frac{\partial B}{\partial \pi_i} - \frac{\partial A}{\partial \pi_i} \frac{\partial B}{\partial \psi_i} + \sum_j \frac{\partial A}{\partial \bar{\psi}_j} \frac{\partial B}{\partial \bar{\pi}_j} - \frac{\partial A}{\partial \bar{\pi}_j} \frac{\partial B}{\partial \bar{\psi}_j} $$
but note that this is only formally valid, because the variables are not all independent. The equations of motion would be written somewhat as
$$ \dot{A} \approx \{ A, \mathcal{H} \} $$
using Dirac weak equality sign, because this is a constrained dynamics. The Hamiltonian density is obtained from
$$ \mathcal{H} \approx \sum_i \pi_i \dot{\psi}_i + \sum_j \bar{\pi}_j \dot{\bar{\psi}}_j - \mathcal{L} $$
Notice that all of this is a quantum treatment. There is not classical spinor theory.
[*] I suppose that all depends on what are you trying to do.
It is possible to construct a Hamiltonian. In fact this is the way how Dirac initially wrote his equation. For that the time and space coordinates have to be treated differently. In the Schrödinger picture, the Hamiltonian generates the time dynamics via ($\hbar =0$) $$i \partial_t \psi = H \psi.$$ We see that we can obtain this structure from the Dirac equation by multiplying it by $\gamma^0$ (and using $\gamma_0^2=1$). With that we obtain the result $$ H = {\bf\alpha} \cdot {\bf p} + \beta m$$ where we introduced the conventional notation $\beta =\gamma^0$, $\alpha^k =\gamma^0 \gamma^k$, and $p_k = -i\partial_k$.
If you want to write this down in terms of a classical field theory, with the field $\psi$ evolving as $$\partial_t \psi({\bf r}) = \{\mathcal{H},\psi ({\bf r})\}, $$ the Hamiltonian is given by $$\mathcal{H} = \int\!d^3r \,\psi^*({\bf r})H \psi({\bf r}).$$
Edit:
I define the Poisson bracket to be $$\{ A, B\} = \int\!d^3r\left[\frac{\delta A}{\delta \psi({\bf r})}\frac{\delta B}{\delta \psi^*({\bf r})} -\frac{\delta B}{\delta \psi({\bf r})}\frac{\delta A}{\delta \psi^*({\bf r})} \right]$$ where the derivatives are functional derivatives and we have assumed (as usual) that $\psi$ and $\psi^*$ are independent variables with the defining relations $$ \frac{ \delta \Psi({\bf r})}{\delta \Psi({\bf r}')} = \frac{ \delta \Psi^*({\bf r})}{\delta \Psi^*({\bf r}')} = \delta^3({\bf r}-{\bf r}'), \quad\frac{ \delta \Psi({\bf r})}{\delta \Psi^*({\bf r}')} = \frac{ \delta \Psi^*({\bf r})}{\delta \Psi({\bf r}')} = 0.$$
