The hierarchically orthogonal functional decomposition of any measurable function f of a random vector X=(X_1,...,X_p) consists in decomposing f(X) into a sum of increasing dimension functions depending only on a subvector of X. Even when X_1,..., X_p are assumed to be dependent, this decomposition is unique if components are hierarchically orthogonal. That is, two of the components are orthogonal whenever all the variables involved in one of the summands are a subset of the variables involved in the other. Setting Y=f(X), this decomposition leads to the definition of generalized sensitivity indices able to quantify the uncertainty of Y with respect to the dependent inputs X. In this paper, a numerical method is developed to identify the component functions of the decomposition using the hierarchical orthogonality property. Furthermore, the asymptotic properties of the components estimation is studied, as well as the numerical estimation of the generalized sensitivity indices of a toy model. Lastly, the method is applied to a model arising from a real-world problem.