A method for calculating the EPR g-tensor based on coupled perturbed Hartree-Fock (HF) and density functional theory (DFT) is presented. The one-electron molecular orbitals of a spin-unrestricted Slater determinant are calculated up to first order in the applied magnetic field. The g-tensor is evaluated as a mixed second derivative property with respect to the applied field and the electron magnetic moment. Thus, spin-polarization and spin-orbit coupling are simultaneously included in the calculation. The treatment focuses on orbitally nondegenerate molecules but is valid for a general ground state spin S and, for the first time, it is possible to include hybrid density functionals in the treatment. The relativistic mass and diamagnetic gauge corrections are also considered. An implementation of the theory is described. Extensive numerical calculations for a series of small molecules are reported with the Hartree-Fock (HF) method, the local density approximation (LSD), the generalized gradient approximation (GGA) and hybrid density functionals such as B3LYP and PBE0 and large Gaussian basis sets. Detailed comparison with available ab initio and DFT calculations are made. The results indicate that the hybrid functionals offer little or no improvement over the GGA functionals for small radicals made of light atoms. For transition metal complexes the situation is different. The hybrid functionals give, on average, better results than the GGA functionals but significant disagreement between theoretical and experimental g-shifts still remain. Overall, the results indicate that the present method is among the most accurate so far developed models for the prediction of g values.