Abstract This computational work deals with the optimal design of the thermal ablation treatment of skin cancer, by considering uncertainties in the model parameters. The tumor and other tissues were heated by a laser. Nanoparticles were used to improve the effects of the heating procedure and to promote thermal damage localized in the region containing the tumor. Treatment protocols examined in this work involved one single heating session with different prespecified durations, where the design variables were considered as the volume fraction of nanoparticles in the epidermis and tumor, as well as the time variation of the incident laser fluence rate. The optimal design problems were solved with the Markov Chain Monte Carlo method, by applying a modified version of the Metropolis-Hastings algorithm with sampling by blocks of parameters. The two parameter blocks were given by the properties of the tissues and by the design variables. The prior for the volume fraction of nanoparticles was given by a truncated Gaussian distribution, while a noninformative Gaussian Markov random field prior was used for the time variation of the laser fluence rate. The posterior distributions of the design variables were estimated by taking into account uncertainties in the model parameters and the desired statistical distribution of the thermal damage in the region of interest. The stochastic simulations resulted in optimal thermal damages with small uncertainties, which closely followed their desired statistical distribution functions.