We study the usefulness of the permanent state as variational wave functions for bosons, which is the bosonic counterpart of the Slater determinant state for fermions. For a system of $N$ identical bosons, a permanent state is constructed by taking a set of $N$ arbitrary (not necessarily orthonormal) single-particle orbitals, forming their product and then symmetrizing it. It is found that for the one-dimensional Bose-Hubbard model with the periodic boundary condition and at unit filling, the exact ground state can be very well approximated by a permanent state, in that the permanent state has high overlap (at least 0.96 for 12 particles and 12 sites) with the exact ground state and can reproduce both the ground state energy and the single-particle correlators to high precision. For more general models, we have devised an optimization algorithm to find the optimal set of single-particle orbitals to minimize the variational energy or maximize the overlap with a given state. It turns out that quite often the ground state of a bosonic system can be well approximated by a permanent state by all the criterions of energy, overlap, and correlation functions. And even if the error is apparent, it can often be remedied by including more configurations, i.e., by allowing the variational wave function to be the superposition of multiple permanent states.