Abstract One of the big questions in quarkonia production is how heavy-quark pairs are produced in the initial hard scattering and hadronize into final quarkonia states. This question also affects the estimation of quarkonia breakup cross-sections and color screening in heavy ion collisions. Many existing models describe well general features of quarkonia production like cross-section and transverse momentum distributions. In order to differentiate between various models one has to study quarkonia production in more details. One of such tools is a study of angular distributions of decay leptons, usually called polarization. The PHENIX experiment has measured inclusive J / ψ polar and azimuthal angular decay coefficients in the mid ( | y | < 0.35 ) and forward ( 1.2 < | y | < 2.2 ) rapidity regions in p+p collisions at 200 GeV and 510 GeV. This talk will present the analysis strategy as well as the results in different polarization frames.