I have a function
$f_n(x) = \frac{1}{3} + 2x \left( \frac{1}{3x} \right)^n$
Where n is a natural number larger or equal than 3 and a constant.
I am sampling the value of x numerically in the interval $x \in [\frac{1}{3}, 1]$ a gazillion times and thus can numerically compute the average $<f>_n$ after a sufficiently large number of evaluations. This way, I can construct the function $f(n) = \langle f\rangle_n$ which is the average value of $f$ for a given $n$ when sampling $x$ sufficiently often. Now, I can compute the average of $f(x)$ numerically, for example, for n = 3 I get $\langle f\rangle_3 = 0.629608$.
However, I should be able to compute this analytically too, which is where I'm failing at. Initially, I tried the approach:
$\langle f\rangle_n = f_n(\langle x\rangle) = f_n(\frac{2}{3}) = 0.5 \quad$ (for n = 3)
Since the average value of x in the above given interval is $\frac{2}{3}$. As you can see, this does not match the numerical result.
Lastly, I tried the mean value theorem:
$\langle f\rangle_n = \frac{ \int_{\frac{1}{3}}^{1} f_n(x) \,dx }{1-\frac{1}{3}} = \frac{(n-1)*3^n-9}{(3n-6)*3^n} = 0.5555 \quad $ (for $n = 3$)
which also doesn't match the numerical result. Do you know where I am wrong-headed?
The numerical result has been obtained by evaluating $f_n(x)$ a sufficiently large number of times and simply computing the average value obtained when sampling $x$ in the above given interval.
By LOTUS, if $x$ has density $p$, the mean of a function $f$ is
$$E[f(x)]=\int_{\mathbb{R}} f(x)p(x)dx,$$
yours being the case where $$f(x)=\frac{1}{3} + 2x \left( \frac{1}{3x} \right)^n\\ p(x)=\frac{1}{1-1/3}\mathbf{1}_{x\in [1/3,1]}$$
assuming you draw $x$ uniformly on $[1/3,1].$
Your computation using the integral is exactly this, and the law of large numbers predicts you should get close to this by simulation. Check your code.