In this paper, we deal with the challenging computational issue of interpolating large data sets, with eventually non-homogeneous densities. To such scope, the Radial Basis Function Partition of Unity (RBF-PU) method has been proved to be a reliable numerical tool. However, there are not available techniques enabling us to efficiently select the sizes of the local PU subdomains which, together with the value of the RBF shape parameter, greatly influence the accuracy of the final fit. Thus here, by minimizing an \emph{a priori} error estimate, we propose a RBF-PU method by suitably selecting variable shape parameters and subdomain sizes. Numerical results and applications show performaces of the interpolation technique.