Short answer: Modern implementations of these two methods lead to similar accuracies.
Longer answer: The calculation of phonons requires the calculation of the Hessian of the potential energy surface $V(\mathbf{R})$, also known as the matrix of force constants:
$$
\frac{\partial^2 V(\mathbf{R})}{\partial \mathbf{R}_i\partial\mathbf{R}_j}=-\frac{\partial \mathbf{F}_j}{\partial\mathbf{R}_i},
$$
where $\mathbf{R}$ is a collective coordinate of all atomic positions, $i$ and $j$ label atoms $i$ and $j$ in your system, and $\mathbf{F}_j=-\partial V(\mathbf{R})/\partial\mathbf{R}_j$ is the force felt when displacing atom $j$.
Finite displacement. This is what you call "frozen phonon", and in this method you calculate the forces in DFT, and then calculate the derivative of the forces by finite difference methods. Therefore, the numerical approximation is that of approximating a derivative with a finite difference formula. In principle, you can make this calculation as accurate as you want by using increasingly accurate approximations to the numerical derivative, but in practice even low-order approximations lead to very accurate answers. The advantages of this method are that it is very simple to implement, and therefore it is in fact available using any underlying electronic structure method that can calculate forces, which includes semilocal DFT, hybrid DFT, or other non-DFT methods, like force fields or dynamical mean-field theory. The disadvantage of this method is that it requires the construction of supercells to capture long-wavelength phonons, which can make the calculations expensive. Most finite displacement codes use "diagonal" supercells, which lead to poor scaling, but the recently introduced "nondiagonal" supercells here [disclaimer: I am a co-author of this work], significantly reduce the computational cost.
DFPT. In density functional perturbation theory, the calculation of the Hessian is specialized to DFT, and the second derivative of the energy is calculated as
$$
\frac{\partial^2 E}{\partial\lambda_i\lambda_j}=\int\frac{\partial^2V(\mathbf{r})}{\partial\lambda_i\lambda_j}n(\mathbf{r})d\mathbf{r}+\int\frac{\partial n(\mathbf{r})}{\partial\lambda_i}\frac{V(\mathbf{r})}{\partial\lambda_j}d\mathbf{r}.
$$
This expression is general for parameters $\lambda$, and in the case of phonons they simply become the atomic coordinate $\lambda_i=\mathbf{R}_i$. This expression requires the calculation of the derivative of the density, which in turn requires the derivative of the Kohn-Sham states. These terms can all be calculated within DFPT with the usual numerical approximations of finite basis sets, etc. Again, in practice, modern implementations are relatively easy to converge. The advantage of DFPT is that it does not require the construction of supercells, one can build a finite wave vector response within the primitive cell, so the computational cost is smaller than in the finite displacement method. The disadvantage of this method is that it is restricted to DFT (so no DMFT for example), and furthermore, the algorithmic implementation is not trivial, so it is only widely available for semilocal DFT (so no hybrid DFT either).
In summary, these two methods lead to comparable accuracies. If DFPT is available, then the calculations will be cheaper and DFPT should be the method of choice. However, DFPT is only widely available with semilocal DFT, which means that if you want phonons at the hybrid functional level, or using beyond-DFT methods, then you have to use the finite displacement method.
