This computational work deals with the optimal design of the thermal ablation treatment of skin cancer, by considering uncertainties in the model parameters. Heating of the tissues that include a tumour was promoted by a diode laser. Nanoparticles were supposed uniformly distributed in the epidermis and in the tumor, to concentrate the thermal damage in the region of interest. The heating protocols examined in this work involved one single treatment session with pre-specified durations, where the design variables were considered as the volume fraction of nanoparticles in the epidermis and tumour, along with 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, respectively. The prior for the volume fraction of nanoparticles was given by a truncated Gaussian distribution, while a non-informative Gaussian Markov Random Field prior was used for the time variation of the laser fluence rate. The posterior distribution of the design variables was estimated for the protocols examined, taking into account uncertainties in the model parameters and the desired statistical distribution of the thermal damage in the region of interest. The stochastic simulation based on the Metropolis-Hastings algorithm with sampling by blocks of parameters resulted in optimal thermal damages that followed the desired probability distribution function with small uncertainties.