Let $f$ be a strictly positive function such that $\int_{0}^{\infty}f(x)dx=\int_{0}^{\infty}xf(x)=1$ (i.e., a probability density function with expectation one). Let also $g$ be a nonnegative nonconstant function which satisfies $\int_{0}^{\infty}g(ax)f(x)dx=a$, $\forall a>0$. Does this imply that $g(x)=x$ a.e.?

The answer is no. Let $h$ be a positive function on $R$ with the following properties: a) Fourier transform has a pair of nonreal zeros $\lambda, \overline{\lambda}$, symmetric with respect to the imaginary axis. b) $\int_{\infty}^\infty h(t)dt=1,$ c) $\int_{\infty}^\infty e^t h(t)dt=1.$ Condition a) means $$\int_{\infty}^\infty e^{i\lambda t}h(t)dt=0.$$ Multiplying this by $a^{i\lambda}$, $a>0$ we obtain $$\int_{\infty}^\infty e^{i\lambda(t+\log a)}h(t)dt=0,$$ for every $a>0$. Making change of the variable $t=\log x$, we obtain $$\int_0^\infty e^{i\lambda\log (ax)}h(\log x)\frac{dx}{x}=0.$$ Putting $f(x)=h(\log x)/x$, we obtain a function with the properties you stated; they follow from b), c). Thus $$\int_0^\infty (ax)^{i\lambda} f(x)dx=0.$$ Now $g(x)=x+x^{i\lambda}$ is a function $g$ which satisfies your identity. If a positive function is required, we need $\lambda=\sigmai$, and $$g(x)=x+k(x^{i\lambda}+x^{i\overline{\lambda}})=x(1+k\cos(\sigma\log x)),$$ which is positive if $0<k<1$. Existence of such $h$ is pretty evident. Take any probability density whose Fourier transform is analytic in some strip $Im z<c$ and has some zeros with negative imaginary part. Scaling will give you a zero whose imaginary part is $1$. To achieve c), shift $h$, replacing it with $h(t+c)$. This does not affect the zeros of Fourier transform. Taking $h$ with Fourier transform having infinitely many zeros on a horizontal line, you obtain an infinite dimensional space of $g$ with fixed $f$. Edit. Sorry, I did not notice that you also need $g$ positive. I modified the example to achieve this additional property. 

