Based on the fermionic mean field theory for spins, we develop the Gutzwiller projection method for the BCS-type mean field states of S = 1 systems. Through variational Monte Carlo (VMC) approach, we can use

the projected wavefunctions to approximate the ground state. We study the S = 1 bilinear-biquadratic model $H=\sum_i J\mathbf S_i\cdot\mathbf S_{i+1} +K(\mathbf S_i\cdot\mathbf S_{i+1})^2$  and find that the optimized projected wavefunction is very close to the ground state when there is gap. The topology of the mean field state is important. After projection the topologically non-trivial states (weak pairing states) correspond to the Haldane phase, while the topologically trivial states (strong pairing states) fall in the dimerized phase. The transition point obtained from VMC is very close to the exact result. The correspondence between weak/strong pairing and Haldane/dimerized phase can be explained by the Z2 gauge uctuations around the mean field ground state.