Fix a probability space $(\Omega, \mathcal{A}, \Bbb P),$ a continuously differentiable function $f:\Bbb R^n \rightarrow \Bbb R,$ and a random vector $X: \Omega \rightarrow \Bbb R^n.$ Furthermore, we define the random variable $Y: \Omega \rightarrow \Bbb R$ as $Y:= f \circ X.$
Note that, if $\Bbb P_X$ denotes the distribution of $X$ on $(\Bbb R^n, \mathcal{B}),$ then $f$ can be seen as a random variable defined in the probability space $(\Bbb R^n, \mathcal{B}, \Bbb P_X).$ Similarly, the function $\nabla f: \Bbb R^n \rightarrow \Bbb R^n$ is random vector defined on that space. I am interested in the expectation of this random vector, i.e., $E_{X \sim \Bbb P_X} [\nabla f].$
Suppose now that an i.i.d random sample $X_1, \ldots, X_n$ of $\mathbb{P}_X$ is observed, together with the values $Y_i := f \circ X_i.$ Thus, in particular, the function $f$ remains unknown. Is there a way to construct a good (unbiased, consistent) estimator for $E_{X \sim \Bbb P_X} [\nabla f]$?
Without further assumptions on $\mathbb{P}_X$ you can't do anything. As a counterexample, consider $\mathbb{P}_X = \delta_0$ (i.e., $\mathbb{P}_X\big[\{0\}\big] =1$). Now, if your sample size is $m \in \mathbb{N}$, any statistics $h_m :(\mathbb{R}^{n}\times\mathbb{R})^m \to \mathbb{R}^n$ is such that \begin{equation*} \mathbb{E}_\mathbb{P}\bigg[h_m\Big(\big(X_1,f(X_1)\big),\dots ,\big(X_m,f(X_m)\big)\Big)\bigg] = h_m\Big(\big(0,f(0)\big),\dots,\big(0,f(0)\big)\Big). \end{equation*} In particular, any estimator of the form $h_m\Big(\big(X_1,f(X_1)\big),\dots ,\big(X_m,f(X_m)\big)\Big)$ has expected value that is a function of $f(0)$ alone, which tells you nothing about $\nabla f(0) = \mathbb{E}_{\mathbb{P}}[\nabla f(X)]$.