We use very deep near-infrared (NIR) imaging data obtained in MOIRCS Deep Survey (MODS) to investigate the evolution of the galaxy stellar mass function back to z~3. The MODS data reach J=24.2, H=23.1, K=23.1 (5sigma, Vega magnitude) over 103 arcmin^2 (wide) and J=25.1, H=23.7, K=24.1 over 28 arcmin^2 (deep) in the GOODS-North region. The wide and very deep NIR data allow us to measure the number density of galaxies down to low stellar mass (10^9-10^10 Msun) even at high redshift with high statistical accuracy. The normalization of the mass function decreases with redshift and the integrated stellar mass density becomes ~ 8-18% of the local value at z~2 and ~ 4-9% at z~3, which are consistent with results of previous studies in general fields. Furthermore, we found that the low-mass slope becomes steeper with redshift from alpha ~- 1.3 at z~1 to alpha ~- 1.6 at z~3, and that the evolution of the number density of low-mass (10^9-10^10 Msun) galaxies is weaker than that of M* (~10^11 Msun) galaxies. This indicates that the contribution of low-mass galaxies to the total stellar mass density has been significant at high redshift. The steepening of the low-mass slope with redshift is opposite trend expected from the stellar mass dependence of the specific star formation rate reported in previous studies. The present result suggests that the hierarchical merging process overwhelmed the effect of the stellar mass growth by star formation and was very important for the stellar mass assembly of these galaxies at 1<~z<~3.