Experimental data obtained for the polarized Bjorken sum rule (BSR) Gamma _1^{p-n}(Q^2) are fitted by using predictions derived within a truncated operator product expansion (OPE) approach to QCD. Four QCD versions are considered: perturbative QCD (pQCD) in the {overline{mathrm{MS}}} scheme, Analytic Perturbation Theory (APT), and 2delta and 3delta analytic QCD versions. In contrast to pQCD, these QCD variants do not have Landau singularities at low positive Q^2, which facilitates the fitting procedure significantly. The fitting procedure is applied first to the experimental data of the inelastic part of BSR, and the known elastic contributions are added after the fitting. In general, when 2delta and 3delta QCD coupling is used the fitted curves give the best results, within the Q^2-range of the fit as well as in extended Q^2-intervals. When the fitting procedure is applied to the total BSR, i.e., to the sum of the experimental data and the elastic contribution, the quality of the results deteriorates significantly.