In his classic book, Lectures on Quantum Mechanics [1], Dirac introduces a generic theory of many particles through the theory's Lagrangian $L$, which, assuming (as usual) that $L$ depends only on the particle positions $q_n$ and velocities $\dot{q}_n$, leads directly to the Euler-Lagrange equations of motion, obtained by minimizing the action $S=\int dt L$:
 
(1)\begin{align}
\frac{d}{dt}\left(\frac{\partial L}{\partial\dot{q}_n}\right)=\frac{\partial L}{\partial q_n}.\label{eq:EL}
\end{align}
Dirac then switches to the Hamiltonian formalism by introducing the canonical momenta
 
(2)\begin{align}
p^n=\frac{\partial L}{\partial\dot{q}_n}.\label{eq:p}
\end{align}
(My preference is to be fussy and always distinguish contravariant and covariant indices. Dirac used lower indices in all expressions.)

Dirac then moves on to introduce, what he calls, the primary constraints of the theory, given in the form
 
(3)\begin{align}
\phi_m(p^n,q_n)=0.\label{eq:constraint}
\end{align}
Dirac then defines the Hamiltonian the usual way, given by
\begin{align}
H=p^n\dot{q}_n-L,
\end{align}
where the Einstein summation convention ($x^iy_i\equiv\sum_ix^iy_i$) applies.

The variation of the Hamiltonian is given by
\begin{align}
\delta H=\delta p^n\dot{q}_n-\frac{\partial L}{\partial q_n}\delta q_n.
\end{align}

Dirac also observes that, as a consequence of the constraints $\phi_m$, the Hamiltonian is not uniquely determined; any linear combination of the $\phi_m$ can be added to $H$ to go over to another Hamiltonian,
\begin{align}
H^*=H+c^m\phi_m,
\end{align}
where the quantities $c_m$ can be functions of the $p$-s and $q$-s.

Inexplicably, however, Dirac next "deduces" (his word) the Hamiltonian equations of motion in the form,
 
 
 
 
(4)\begin{align}
\dot{q}_n&=\frac{\partial H}{\partial p^n}+u^m\frac{\partial\phi_m}{\partial p^n},\\
-\dot{p}^n=-\frac{\partial L}{\partial q_n}&=\frac{\partial H}{\partial q_n}+u^m\frac{\partial\phi_m}{\partial q_n},\label{eq:H2}
\end{align}
where the unknown coefficients $u_m$ are now not functions of the $p$-s and $q$-s. In the paragraph preceding these equations, Dirac actually indicates that the $p$-s and the $q$-s cannot be varied independently, because of the constraints (3).

Notice, however, that on the left-hand side of Eq. (4), Dirac uses the Euler-Lagrange equation (1), which was obtained from the unconstrained Lagrangian $L$. This is clearly self-contradictory: the system is either constrained by the $\phi_m$ or it isn't, but in either case, the constraints should be consistently applied (or not applied).

The Lagrangian of the constrained system is obtained by introducing the constraints in the form of Lagrange-multipliers. Specifically, we replace the Lagrangian $L(q_n,\dot{q}_n)$ with
\begin{align}
L^*(q_n,\dot{q}_n,\lambda^m)=L(q_n,\dot{q}_n)+\lambda^m\psi_m(q_n,\dot{q}_n),\label{eq:L*}
\end{align}
where the $\lambda_m$ are not functions of the $q$-s and $\dot{q}$-s, nor are they constants: they are new (nondynamical) degrees of freedom. Varying $L^*$ with respect to the $\lambda^m$ yields the constraint equations
\begin{align}
\psi_m(q_n,\dot{q}_n)=0,
\end{align}
while variation with respect to the $q$-s yields
\begin{align}
\frac{d}{dt}\left(\frac{\partial L^*}{\partial\dot{q}_n}\right)=\frac{\partial L^*}{\partial q_n}.\label{eq:EL*}
\end{align}

Although we have a new Lagrangian, it does not follow automatically that we can perform a Legendre transformation unambiguously and obtain a corresponding new Hamiltonian. However, we can express the modified action $S^*=\int dt L^*$ using the original Hamiltonian and canonical momenta, which are defined through a Legendre transformation of $L$ and their existence is not dependent on the Euler-Lagrange equations (1):
\begin{align}
S^*=\int~dt\left[p^n\dot{q}_n-H+\lambda^m\psi_m\right].
\end{align}

We note that Eq. (2) can always be solved for the $\dot{q}_n$:
\begin{align}
\dot{q}_n=v_n(p^n,q_n),
\end{align}
provided the Hessian metric,
\begin{align}
H^{ij}=\frac{\partial^2L}{\partial\dot{q}_i\partial\dot{q}_j},
\end{align}
is non-singular [2]. This allows us to write
\begin{align}
\psi_m(q_n,\dot{q}_n)=\psi_m(q_n,v_n(p^n,q_n))=\phi_m(p^n,q_n),
\end{align}
and write the action in the form,
\begin{align}
S^*=\int~dt\left[p^n\dot{q}_n-H+\lambda^m\phi_m\right].
\end{align}
The variation of this action is given by
\begin{align}
\delta S^*=\int&~dt\left[\delta p^n\dot{q}_n+p^n\delta\dot{q}_n-\frac{\partial H}{\partial p^n}\delta p^n-\frac{\partial H}{\partial q_n}\delta q_n+\delta\lambda^m\psi_m+\lambda^m\frac{\partial \phi_m}{\partial q_n}\delta q_n+\lambda^m\frac{\partial \phi_m}{\partial p^n}\delta p^n\right].\nonumber
\end{align}
Proceeding the usual way, we observe that as the operators $d/dt$ and $\delta$ commute, we can eliminate $\delta\dot{q}$ at the cost of introducing a surface term, which does not contribute to the variation of $S^*$:
\begin{align}
p^n\delta\dot{q}_n=p^n\frac{d}{dt}\delta q_n=\frac{d}{dt}(p^n\delta q_n)-\dot{p}^n\delta q_n,
\end{align}
leading us to
\begin{align}
\delta S^*=\int&~dt\left[\delta p^n\dot{q}_n-\dot{p}^n\delta q_n-\frac{\partial H}{\partial p^n}\delta p^n-\frac{\partial H}{\partial q_n}\delta q_n+\delta\lambda^m\psi_m+\lambda^m\frac{\partial \phi_m}{\partial q_n}\delta q_n+\lambda^m\frac{\partial \phi_m}{\partial p^n}\delta p^n\right],\nonumber
\end{align}
or, after collecting like terms,
\begin{align}
\delta S^*=\int~dt&\left\{\delta p^n\left[\dot{q}_n-\frac{\partial H}{\partial p^n}+\lambda^m\frac{\partial \phi_m}{\partial p^n}\right]-\left[\dot{p}^n+\frac{\partial H}{\partial q_n}-\lambda^m\frac{\partial \phi_m}{\partial q_n}\right]\delta q_n+\delta\lambda^m\psi_m\right\},\nonumber
\end{align}
leading to the Hamiltonian equations of motion,
\begin{align}
\dot{q}_n&=\frac{\partial H}{\partial p^n}-\lambda^m\frac{\partial \phi_m}{\partial p^n},\\
-\dot{p}^n&=\frac{\partial H}{\partial q_n}-\lambda^m\frac{\partial \phi_m}{\partial q_n},\\
\phi_m(p^n,q_n)&=0.
\end{align}
These are the same equations (after equating $\lambda^m=-u^m$) that Dirac "deduced". Clearly, Dirac's instincts were in the right place (unsurprisingly, given his reputation.) However, we obtained this result without resorting to ill-defined concepts such as "strong" and "weak" equations, and without any heuristic appeals to deduction, using rigorous mathematics instead. Furthermore, the role of the $\lambda_m$ as Lagrange-multipliers is now clear, and no use is made of the unconstrained Euler-Lagrange equation (1), which is clearly not satisfied in the constrained system.

[1] P. A. M. Dirac, Lectures on Quantum Mechanics (Dover Publications, 2001).
[2] H. Kleinert, Path Integrals in Quantum Mechanics, Statistics, Polymer Physics, and Financial Markets (World Scientific, 2006), 4th ed.