To improve the accuracy of Monte Carlo estimation of expectations, a set of zero-mean feature functions, known as control variates, can be used. They can be used as feature functions for linear regression of the target function, and we can obtain an unbiased and variance-reduced estimate using its residual. One known way to construct such functions is a method using an equality called Stein's identity, but these functions are not sufficient for the case where the target distribution is multimodal. We propose a different approach to constructing these zero-mean functions based on distribution approximation and the density ratio. We demonstrate that combining the functions constructed by these two strategies can effectively reduce the estimation variance for a bimodal distribution.