  We determine the nuclear modifications of parton distribution functions of
bound protons at scales $Q^2\ge 1.69$ GeV$^2$ and momentum fractions
$10^{-5}\le x\le 1$ in a global analysis which utilizes nuclear hard process
data, sum rules and leading-order DGLAP scale evolution. The main improvements
over our earlier work {\em EKS98} are the automated $\chi^2$ minimization,
simplified and better controllable fit functions, and most importantly, the
possibility for error estimates. The resulting 16-parameter fit to the N=514
datapoints is good, $\chi^2/{\rm d.o.f}=0.82$. Within the error estimates
obtained, the old {\em EKS98} parametrization is found to be fully consistent
with the present analysis, with no essential difference in terms of $\chi^2$
either. We also determine separate uncertainty bands for the nuclear gluon and
sea quark modifications in the large-$x$ region where they are not stringently
constrained by the available data. Comparison with other global analyses is
shown and uncertainties demonstrated. Finally, we show that RHIC-BRAHMS data
for inclusive hadron production in d+Au collisions lend support for a stronger
gluon shadowing at $x<0.01$ and also that fairly large changes in the gluon
modifications do not rapidly deteriorate the goodness of the overall fits, as
long as the initial gluon modifications in the region $x\sim 0.02-0.04$ remain
small.
