System properties

k = 2;
m = 1;
x0 = 1;
ES = (1/2) k (x0^2);

Function to be averaged

F[p_, q_] = (p^2);

Time average

ω = (Sqrt[(k/m)]); 
    T = ((2 π)/ω);
p[
          t_] = -m ω x0 Sin[ω t];
q[t_] = x0 Cos[ω t];
Favt = (1/T) Integrate[F[p[t], q[t]], {t, 0, T}]

Phase space average using a volume intregral

pE[
          q_] = (Sqrt[2  m (Es - k (q^2)/2)]);
V[Es_] = 
        Integrate[
          Integrate[
            1, {p, -pE[q], 
              pE[q]}], {q, -(Sqrt[2 Es/k]), (Sqrt[2 Es/k])}];
Ω[Es_] = V'[Es];
Fav = (1/Ω)[Es] (D[
              Integrate[
                Integrate[
                  F[p, q], {p, -pE[q], 
                    pE[q]}], {q, -(Sqrt[2 Es/k]), (Sqrt[2 Es/k])}], Es]) /. 
        Es -> ES