Average as a constant approximation in $L^p$ Let $I=[0,1]$. For $p\in[1,\infty]$ define $C_p$ as the best constant such that for all $f\in L^p(I)$
$$
\left\|f-\int_If\,\right\|_{L^p(I)}\leq C_p\inf_{c\in\mathbb{R}}\left\|f-c\,\right\|_{L^p(I)}.
$$
The following properties are fairly elementary:

*

*$1\leq C_p\leq 2$. The first inequality is trivial, the second follows by convexity: for any $c$
$$
\left\|f-\int_If\,\right\|_{L^p(I)}\leq \left\|f-c\,\right\|_{L^p(I)}+\left\|c-\int_If\,\right\|_{L^p(I)}\leq 2\left\|f-c\,\right\|_{L^p(I)}.
$$

*$C_2=1$. This is a standard argument: to minimize $\left\|f-c\,\right\|_{L^2(I)}^2$ we differentiate with respect to $c$:
$$
\frac{\partial}{\partial c}\left\|f-c\,\right\|_{L^2(I)}^2=2c-2\int_I f,
$$
so the average is the minimizer.

*$C_\infty=2$. No better constant exists for the sequence $f_n=n\mathbf{1}_{[0,1/n]}$. Indeed, $\|f_n-1\|_{L^\infty(I)}=n-1$ and the optimal choice is $\|f_n-n/2\|_{L^\infty(I)}=n/2$.

*$C_1=2$. For the same sequence $f_n=n\mathbf{1}_{[0,1/n]}$, one has $\|f_n-1\|_{L^1(I)}=2(n-1)/n$ and the optimal choice is $\|f_n-0\|_{L^1(I)}=1$.

What can be said about $p\neq 1,2,\infty$? The topic of best constant (or more generally, polynomial) approximation is very classical but I could not find results on this particular question.
Based on the above, I computed numerically the suggested constant by taking $f_{100}$, which gives the corresponding approximate lower bounds

*

*$C_3\geq 1.0553$

*$C_4\geq  1.1465$

*$C_5\geq 1.2340$

*$C_6\geq  1.3096$

*$C_{10}\geq 1.5114$

*$C_{15}\geq   1.6425$

*$C_{25}\geq   1.7645$
 A: Upper bounds
This is not a full answer. But provides some ideas.
Duality
First, I'd like to spell out the duality argument and show that if $\frac1p + \frac1q = 1$, then $C_p = C_q$.
For this we will need the well-known dual characterization of $L^p$ norms:
$$ \|f \|_p = \sup_{g: \|g\|_q = 1} \int fg $$
Now, on the left we want to estimate $\| f- \int f\|_p$, so we compute, for an arbitrary $g\in L^q$ with unit norm the quantity $\int (f - \int f)g$. Observe that since $f-\int f$ has zero mean, its integrals against constants vanish. So we have
$$ \int (f-\int f)g = \int (f-\int f)(g - \int g) $$
Now modifying $f$ using the same reason, we find
$$ = \int (f-c)(g-\int g) \leq \|f-c\|_p \|g - \int g\|_q $$
where we applied, for the inequality, Holder's inequality.
To conclude, we have that for any $f\in L^p$ and $g\in L^q$ that
$$ \int (f-\int f)g \leq \|f - c\|_p \|g - \int g\|_q \leq \|f - c\|_p \cdot C_q \|g - c'\|_q $$
for any $c'\in \mathbb{R}$.
This in particular applies to $c' = 0$.
Now optimizing over all $g$ with $\|g\|_q = 1$ you find
$$ \|f - \int f\|_p \leq C_q \|f - c\|_p $$
for every $f\in L^p$ and every $c\in \mathbb{R}$. This shows that $C_p \leq C_q$.
But the argument is symmetric in $p$ and $q$, thus we conclude that $C_p = C_q$.
Interpolation bounds
First note that the question can be rephrased as looking for the operator norm of the linear operator $T: f \mapsto f - \int f$, since given $g = T(f)$, then $T^{-1}(g) = \{ f + c: c\in \mathbb{R}\}$. So the problem of finding the operator norm of $T$ can be written as finding
$$ \sup_{f\in L^p\setminus \{0\}} \frac{\|Tf\|_p}{\|f\|_p} = 
\sup_{g\in L^p\setminus \{0\}, \int g = 0} \sup_{f \in T^{-1}(g)} \frac{\|g\|_p}{\|f\|_p} = \sup_{g\in L^p\setminus \{0\}, \int g = 0} \frac{\|g\|_p}{\inf_{c} \|g+c\|_p}$$
The fact that $C_p$ can be found explicitly for $p = 1, 2, \infty$ means that you can apply Riesz-Thorin-Stein interpolation.
For $p > 2$, we can write it as
$$ \frac1p = \frac{\theta}{2} + \frac{1-\theta}{\infty}$$
(or more simply, $\theta = 2/p$). Complex interpolation tells us that our operator having operator norm $C_2 = 1$ and $C_\infty = 2$ must have operator norm
$$ C_p \leq C_2^\theta C_\infty^{1-\theta} = 2^{1-\theta} = 2^{1-2/p}$$
In particular, $C_p < 2$ for all $p \neq 1,\infty$.
For comparison, here are the numerical values of the function above:




$p$
$2^{1-2/p}$




3
1.26


4
1.414


5
1.516


6
1.587


8
1.682


10
1.741


15
1.823


20
1.866


25
1.892


50
1.945




So you see that the numerical lower bounds are compatible with the theoretical upper bounds, but there is a gap.
A: Lower Bounds
Let $f_\sigma$ be the function that equals $1/\sigma$ on $[0,\sigma]$ and $0$ otherwise (generalizing your $f_n$), you can actually compute by hand that $\inf \|f_\sigma - c\|_{p}$ is achieved with
$$ c = \frac{1}{\sigma + \sigma (1/\sigma - 1)^{1/(p-1)}} $$
This yields the ratio (and lower bound for $C_p$)
$$ \frac{\|f_\sigma - \int f_\sigma\|_p}{\|f_\sigma - c\|_p} = \left( \sigma^{p-1} + (1-\sigma)^{p-1}\right)^{1/p}\left( \sigma^{1/(p-1)} + (1-\sigma^{1/(p-1)}\right)^{1-1/p} $$
from which you see that your numerical values based on $f_{100}$ are underestimates (rather significantly so when $p\in (2,10]$).
Numerically optimizing the above expression gives you slightly more refined bounds




$p$
lower bound for $C_p$




3
1.095


4
1.212


5
1.302


6
1.374


8
1.478


10
1.55


15
1.662


20
1.727


25
1.769


50
1.867




Additional Numerical Results
The computations above maximized a quotient for $f$ that can be expressed as a step function with at most two distinct output values.
For what it is worth, I wrote up a numerical optimization routine that computed the lower bound for $C_p$ when $f$ can be expressed as a step function with no more than three distinct output values.
Interestingly the optimal values found for $C_p$ are identical to those listed in the table above, and the optimal function found are those that only has two distinct output values (so can be identified as one of the $f_\sigma$ above). So there's a chance that the values above are sharp.
For those interested: the code is below (in julia using JuMP)
The function is normalized to have one (the largest in absolute value) of the three output values to be $-1$ (since the problem is scale invariant). The variables w1, w2, w3 are the widths of the three segments, and h1, h2 are the other two heights of the step function. The variable c is the constant that is to be subtracted. Execute the findCp function with the value of p you want. The last variable in the output is the value of $C_p$ found.
I tried with a few different initial seed function values, and they all seem to converge to the same answer.
using JuMP, Ipopt

function findCp(p)
        model = Model(Ipopt.Optimizer)
        set_silent(model)
        @variable(model, -1<= h1<=1, start = 0.0)
        @variable(model, -1<= h2<= 1, start = 1/2)
        @variable(model, 0<= w1 <= 1, start = 1/4)
        @variable(model, 0<= w2 <= 1, start = 1/2)
        @variable(model, 0<= w3 <= 1, start = 1/4)
        @variable(model, c)
        @constraint(model,  w1+w2+w3 == 1)
        @constraint(model, h1*w1 + h2*w2 - w3 == 0)

        @NLobjective(model, Max, (abs(h1)^p*w1 + abs(h2)^p*w2 + w3) / (abs(h1-c)^p*w1 + abs(h2-c)^p*w2 + abs(1+c)^p*w3) )

        optimize!(model)
        return value(h1), value(h2), value(w1), value(w2), value(w3), value(c), objective_value(model)^(1/p)
end

findCp(3)
````

