Numerical calculation of two-loop electroweak corrections to the muon anomalous magnetic moment ($g$-2) is done based on, on shell renormalization scheme (OS) and free quark model (FQM). The GRACE-FORM system is used to generate Feynman diagrams and corresponding amplitudes. Total 1780 two-loop diagrams and 70 one-loop diagrams composed of counter terms are calculated to get the renormalized quantity. As for the numerical calculation, we adopt trapezoidal rule with Double Exponential method (DE). Linear extrapolation method (LE) is introduced to regularize UV- and IR-divergences and to get finite values. The reliability of our result is guaranteed by several conditions. The sum of one and two loop electroweak corrections in this renormalization scheme becomes $a_mu^{EW:OS}[1{rm+}2{rm -loop}]= 151.2 (pm 1.0)times 10^{-11}$, where the error is due to the numerical integration and the uncertainty of input mass parameters and of the hadronic corrections to electroweak loops. By taking the hadronic corrections into account, we get $a_mu^{EW}[1{rm+}2 {rm -loop}]= 152.9 (pm 1.0)times 10^{-11}$. It is in agreement with the previous works given in PDG within errors.