Starspots are thought to be regions of locally strong magnetic fields, similar to sunspots, and they can generate photometric brightness modulations. To deduce stellar and spot properties, such as spot emergence and decay rates, we implement a computational code for starspot modeling. It is implemented with an adaptive parallel tempering algorithm and an importance sampling algorithm for parameter estimation and model selection in the Bayesian framework. For evaluating the performance of the code, we apply it to synthetic light curves produced with three spots. The light curves are specified in the spot parameters, such as the radii, intensities, latitudes, longitudes, and emergence/decay durations. The spots are circular with specified radii and intensities relative to the photosphere, and the stellar differential rotation coefficient is also included in the light curves. As a result, stellar and spot parameters are uniquely deduced, and the number of spots is correctly determined: the three-spot model is preferable because the model evidence is much greater than that of the two-spot model by orders of magnitude and more than that of the four-spot model by a more modest factor, whereas the light curves are produced to have two or one local minimum during one equatorial rotation period by adjusting the values of longitude. The spot emergence and decay rates can be estimated with error less than an order of magnitude, considering the difference of the number of spots.