The most rigorous way to derive corrections to the real gas equation is the virial expansion. At second order this gives Van der Waals law, but this derivation is completely general and can be pushed to higher order.
My answer assumes a basic knowledge of classical mechanics, statistical physics and analysis.
Hamiltonian of the system
We consider a system of identical particles described by the following Hamiltonian:
$$H=\sum_i\frac{\vec{p}_i^2}{2m}+\sum_{i\lt j}U(|\vec{r}_i-\vec{r}_j|)$$
If you don't know what a Hamiltonian is, you can consider it as the energy of the system.
Canonical partition function for the perfect gas
For a perfect gas, you consider that the particles are non-interacting, so that
$$U(|\vec{r}_i-\vec{r}_j|)=0,\quad\forall i,j.$$
In this particular case, life is easy! In fact you can easily compute the Canonical partition function (this is a definition)
$$Z_\mathrm C=\frac1{N!h^{3N}}\int\mathrm d^N\vec{p}\,\mathrm d^N\vec{r}\,\mathrm e^{-\beta H}$$
which, in the case of a non-interacting gas, is simply
$$Z_\mathrm{PG}=\frac1{N!h^{3N}}\int\mathrm d^N\vec{p}\,\mathrm d^N\vec{r}\,\mathrm e^{-\beta\sum_i\frac{\vec{p}_i^2}{2m}}=\frac{1}{N!}\left(\frac{2 \pi mk_\mathrm BT}{h^2}\right)^{\frac{3N}{2}}$$
In order to compute this integral, you can simply use the definition of Gaussian integrals. This is the partition function of a perfect gas and you can now derive all the properties of this system.
Canonical partition function for an interacting gas
If you want to consider an interacting gas, your partition function becomes
$$Z_\mathrm I=\frac1{N!h^{3N}}\int\mathrm d^N\vec p\,\mathrm d^N\vec r\,\mathrm e^{-\beta\left(\sum_i \frac{\vec{p}_i^2}{2m}+\sum_{i\lt j} U(|\vec{r}_i-\vec{r}_j|)\right)}$$
The integration over momenta can be easily computed as before (using Gaussian integrals). In contrast, the integral
$$Q_N=\int\mathrm d^N\vec r\,\mathrm e^{-\beta\sum_{i\lt j}U(|\vec{r}_i-\vec{r}_j|)}$$
is particulary difficult (even for simple forms of the interacting potential $U$). We can write temporarily the partition function of the system as
$$Z_\mathrm I=Z_\mathrm{PG}Q_N=\frac1{N!}\left(\frac{2\pi mk_\mathrm BT}{h^2}\right)^{\frac{3N}{2}}Q_N$$
Grand Canonical partition function for the interacting gas
Now, as you sais, the critical parameter is the density $\rho$: at low density you recover the perfect gas behavior. Therefore, it seems a good idea to develop our formalism in powers of $\rho$. It is easier to work in the Grand Canonical ensemble (you don't have to worry much about this), where the partition function is given by (this is a definition)
$$Z_\mathrm{GC}=\sum_N^\infty\mathrm e^{\alpha N}Z_\mathrm C$$
This is just a generalization of the Canonical partition function $Z_\mathrm C$. With this definition we have, for the interacting gas,
$$Z_\mathrm{II}=\sum_N^\infty\mathrm e^{\alpha N}\frac1{N!}\left(\frac{2\pi mk_\mathrm BT}{h^2}\right)^{\frac{3N}2}Q_N=\sum_N^\infty\frac{z^N}{N!}Q_N,$$
where we defined for simplicity
$$z=\mathrm e^{\alpha}\left(\frac{2\pi mk_\mathrm BT}{h^2}\right)^{\frac32}$$
Relation between the Grand Canonical partition function and the equation of state
From $Z_\mathrm I$ in the Canonical ensemble we can now compute the ratio $p/(k_\mathrm BT)$. The Grand Canonical potential is defined by
$$\phi=-k_\mathrm BT\ln(Z_\mathrm{II})$$
and this Grand Canonical potential is linked to volume and pressure by
$$\phi=-pV,$$
so that you finally have a link between the Gran Canonical partition function and the pressure:
$$\frac p{k_\mathrm BT}=\frac{1}{V}\ln(Z_\mathrm{II}).$$
Since $Z_\mathrm{II}$ is a development in powers of $z$, the same is valid for $\ln(Z_\mathrm{II})/V$. This gives
$$\frac1V\ln(Z_\mathrm{II})=z+\frac{a_2}{2!}z^2+\frac{a_3}{3!}z^3+\cdots$$
Density
Using an easy argument that I will not show here (unless explicitly required), you can then write the density as
$$\rho=\frac{\langle N\rangle}V=z\frac{\partial}{\partial z}\left(\frac1V\ln(Z_\mathrm{II})\right),$$
where $\langle N\rangle$ is the average number of particles in your gas.
Therefore, for the density we can write
$$\rho=z+a_2z^2+\frac{a_3}{2!}z^3+\cdots$$
Developement of the state equation in powers of $\rho$
Since we said that $\rho$ is our critical parameter we can write
$$\frac p{k_\mathrm BT}=\rho+b_2\rho^2+b_3\rho^3+\cdots.$$
In fact the perfect gas law was the first-order expansion with respect to the density. For a diluite gas, $\rho \ll 1$ and therefore we simply get
$$\frac p{k_\mathrm BT}=\rho.$$
Virial coefficients
It is now easy to link the coefficients $a$ to the coefficients $b$, by equating the powers of $z$. In our case you find:
$$b_2=-\frac{a_2}2,\qquad b_3=-a_2^2-\frac{a_3}3.$$
It now remains to link these coefficients to $Q_N$. By definition of $Z_I$ we can write
$$\ln(Z_\mathrm I)=\ln\left(1+Q_1z+Q_2\frac{z^2}{2!}+Q_3\frac{z^3}{3!}+\cdots\right)$$
By Taylor expansion of the logarithm we get
$$\ln(Z_\mathrm I)=Q_1+Q_2\frac{z^2}{2!}+Q_3\frac{z^3}{3!}-\frac{1}{2}\left(Q_1z+Q_2\frac{z^2}{2}\right)^2+Q_1^3\frac{z^3}{3}+\cdots$$
For $a_1$ and $a_2$ we have
$$a_1=\frac{Q_1}V\qquad a_2=\frac1V\left(\frac{Q_2}2-\frac{Q_1^2}2\right)$$
which gives
$$b_2=\frac1{2V}\left(Q_1^2-Q_2\right).$$
Computation of $Q_i$
Now, if we compute $Q_1$ and $Q_2$ explicitly (which is much easier than computing $Q_N$ for a general $N$) we get
$$Q_1=\int\mathrm d^3\vec r=V$$
$$Q_2=\int\mathrm d^3\vec{r}_1\,\mathrm d^3\vec{r}_2\,\mathrm e^{-\beta U(|\vec{r}_1-\vec{r}_2|)}=Va+V^2$$
Van der Waals law
With the expression for the virial coefficients and the analytical expressions of $Q_1$ and $Q_2$ we finally write
$$b_2=-\frac a2$$
Therefore we finally get Van der Waals law
$$\frac p{k_\mathrm BT}=\rho-\frac a2\rho^2.$$
At the second order in $\rho$ we get Van der Waals law, but this procedure is completely general and can be used to compute higher order corrections to the perfect gas law. Computing the coefficients $Q_N$ is quite cumbersome, but can be performed using the cluster expansion.
