Magnetic responses of superparamagnetic nanoparticles to high-frequency AC magnetic fields with sufficiently large amplitudes are numerically simulated to exactly clarify the phenomena occurring in magnetic particle imaging. When the magnetic anisotropy energy inevitable in actual nanoparticles is taken into account in considering the magnetic potential, larger nanoparticles exhibit a delayed response to alternations of the magnetic fields. This kind of delay is rather remarkable in the lower-amplitude range of the field, where the assistance by the Zeeman energy to thermally activated magnetization reversal is insufficient. In some cases, a sign inversion of the third-order harmonic response was found to occur at some specific amplitude, despite the lack in DC bias magnetic field strength. Considering the attenuation of the AC magnetic field generated in the human body, it is possible that the phases of the signals from nanoparticles deep inside the body and those near the body surface are completely different. This may lead to artifacts in the reconstructed image. Furthermore, when the magnetic/thermal torque-driven rotation of the anisotropic nanoparticles as well as the magnetic anisotropy energy are taken into account, the simulated results show that, once the easy axes are aligned toward the direction of the DC bias magnetic field, it takes time to randomize them at the field-free point. During this relaxation, the third-order harmonic response depends highly upon the history of the magnetic field. This is because non-linearity of the anhysteretic magnetization curve for the superparamagnetic nanoparticles varies with the orientations of the easy axes. This history dependence may also lead to another artifact in magnetic particle imaging, when the scanning of the field-free point is faster than the Brownian relaxations.