Expanding on physicsphile's answer, there is an alternative way of computing the expectation value of $f(q)$, and that is to sum the possible values of $f(q)$ weighted by their respective probabilities as follows.
If $f(q)$ represents a physical dynamical variable, then it is a Hermitian (self-adjoint) operator and can therefore be diagonalized by some set of eigenfunctions $\phi_i(q)$ such that
$$f(q)\phi_i(q) = \lambda_i\phi_i(q),$$
where $\lambda_i$ is the eigenvalue. This set of eigenfunctions is an orthonormal (read: they are orthogonal and normalized) basis for the Hilbert space, and therefore any wave function can be expanded as a linear combination of these eigenfunctions,
$$\psi(q) = \sum_i a_i \phi_i(q).$$
Using your formula for the expected value, we have
$$\langle f(q)\rangle = \int dq~\psi(q)^*f(q)\psi(q) = \int \sum_i a_i^*\phi_i(q)^*f(q)\sum_j a_j\phi_j(q)=\sum_{i,j}a_i^*a_j\int \phi_i(q)^*f(q)\phi_j(q)$$
Now, $f(q)$ acting on $\phi_i(q)$ yields $\lambda_i\phi_i(q)$, and since this eigenvalue is just a number, it can be pulled out of the sum, yielding
$$\langle\phi_i(q)\rangle =\sum_{i,j}a_i^*a_j\lambda_i\int \phi_i(q)^*\phi_j(q).$$
Now, since the eigenfunctions are orthonormal, the integral evaluates to the Kronecker-delta
$$\delta_{ij} = \begin{cases} 1 & i=j\\ 0 & i\neq j \end{cases},$$
in which case the sum collapses to a single sum, yielding
$$\langle\phi_i(q)\rangle =\sum_{i}|a_i|^2\lambda_i.$$
When we recognize $|a_i|^2$ as exactly the expression in physicsphile's answer, you can see that you can interpret the integral you've written as a sum over the eigenvalues (i.e. the possible measured values of $f(q)$!) weighted by their probabilities $|a_i|^2$ in the state $\psi(q)$.
