  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 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 3 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. The number of spots is correctly determined: the 3-spot model
is preferable because the model evidence is much greater than that of 2-spot
models by orders of magnitude and more than that of 4-spot model by a more
modest factor, whereas the light curves are produced to have 2 or 1 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.
