Let $\lambda_n$ be the $n$-dimensional Lebesgue measure, and let $f\in L^1_{loc}(\mathbb R^n)$ and $Mf$ be the maximal function of $f$. For fixed $R>0$ and $x \in \mathbb R^n$ and we define $$ I(x)= \int_{B(x,R)}\frac{ f(y)}{|x-y|^{n-1}} d\lambda_n(y).$$
Prove that $$ |I(x)|\leq CR \cdot (Mf)(x)$$ for each $x\in \mathbb R^n$, where $C=C(n)>0$ is a constant only depending on $n$ and where $Mf$ is the Hardy-Littlewood maximal function of $f$.
My first step was to express $\int_{B(x,R)} \frac{|f(y)|}{|x-y|^{n-1}} d\lambda_n(y)$ as a double integral and using Fubini-Tonelli to change the order of integration. However, I was not able to get further. I know that $\lambda_n(\{Mf > \alpha \}) \leq \frac{5^n}{\alpha} ||f||_1$ and that for $||Mf||_p \leq C ||f||_p$ where $C > 0$ depends on $n$ and $p$ and where $1 < p \leq \infty$, but I don't see how I can apply these results for the above problem.
Also, I just wanted to ask if the order of inequality we are trying to prove ($I(x)= \int_{B(x,R)}\frac{ f(y)}{|x-y|^{n-1}} d\lambda_n(y)$) is correct. I got that, if $f$ is nonnegative, then the reverse inequality holds for some $C$.
There are two useful principles we can use here.
If you're trying to bound a quantity by the maximal function of $f$, you should try to bound your quantity by averages of $f$, since the maximal function is just a supremum over averages.
When dealing with integrals with singularities, a reasonable strategy is to break up the integral into pieces so you can leverage different behavior of your integrand near and away from the singularity.
In this case, the singularity is at $y=x$. To apply principle 2, let's consider breaking up the domain of integration by $$ B(x,R) = \{x\}\cup\left(\bigcup_{m\geq 0} A_m(x,R)\right) $$ where $A_m(x,R) = \{y\in\mathbb{R}^n: 2^{-(m+1)}R \leq |x-y| \leq 2^{-m}R\}$. For $y\in A_m(x,R)$, we have $$ |x-y|^{-(n-1)} \leq \left(\frac{2^{m+1}}{R}\right)^{n-1}. $$ Thus on each annulus we have $$ \int_{A_m(x,R)} \frac{|f(y)|}{|x-y|^{n-1}}dy \leq \int_{A_m(x,R)} |f(y)|\left(\frac{2^{m+1}}{R}\right)^{n-1}dy. $$ We can now bound this by expanding the region of integration to $B(x,2^{-m}R)$: $$ \int_{A_m(x,R)} |f(y)|\left(\frac{2^{m+1}}{R}\right)^{n-1}dy \leq \int_{B(x,2^{-m}R)} |f(y)|\left(\frac{2^{m+1}}{R}\right)^{n-1}dy. $$ Next, to use principle 1, we work in the average of $|f|$ over this ball, by multiplying and dividing by its volume $|B(x,2^{-m}R)| = C_n (2^{-m}R)^n$: $$ \int_{B(x,2^{-m}R)} |f(y)|\left(\frac{2^{m+1}}{R}\right)^{n-1}dy \leq \frac{C_n}{|B(x,2^{-m}R)|}\int_{B(x,2^{-m}R)} |f(y)|(2^{-m}R)^n\left(\frac{2^{m+1}}{R}\right)^{n-1}dy. $$ Writing $$ (2^{-m}R)^n\left(\frac{2^{m+1}}{R}\right)^{n-1} = 2^{n-1}\cdot 2^{-m}R, $$ absorbing the constant $2^{n-1}$ into $C_n$, and taking the supremum over $x$, we find that $$ \int_{A_m(x,R)} \frac{|f(y)|}{|x-y|^{n-1}}dy \leq \frac{C_nR2^{-m}}{|B(x,2^{-m}R)|}\int_{B(x,2^{-m}R)} |f(y)|dy \leq C_nR2^{-m}Mf(x). $$ Now taking the sum over all $m\geq 0$, we find that (still absorbing constants) $$ \int_{B(x,R)} \frac{|f(y)|}{|x-y|^{n-1}}dy \leq \sum_{m\geq 0}\int_{A_m(x,R)} \frac{|f(y)|}{|x-y|^{n-1}}dy \leq C_nRMf(x)\cdot\sum_{m\geq 0} 2^{-m} \leq C_nRMf(x). $$
Incidentally, the choice of a dyadic decomposition for the ball was not merely a lucky guess. It is often a good idea to decompose regions of integration so that the integrand varies by a constant factor (independent of the summation index) within the decomposed regions, as it sometimes lets you pull out geometric factors. $|x-y|^{-(n-1)}$ varies within $A_m(x,R)$ by a factor of $2^{n-1}$, i.e. independent of $m$; and in fact, imposing such a requirement on the decomposition quickly leads you to consider concentric annuli with geometrically decaying radii.