I am just wondering your help which R function or code should I use to model a three stage hierarchical model for a count data (a Poisson Process With Change Point)?
Is MCMC the best choice to estimate the time when the rate of an event changed? What should be specified in the MCMC model?

It will be great if I can get some worked examples too.
Regards

I think Bayesian inference is the best way to perform inference for your model. MCMC samplers are often used to perform Bayesian inference, in cases like this where you don't have a conjugate model. This example uses a Gibbs sampler:

However, I'm not sure a Gibbs sampler is completely appropriate for this problem. In case you want to test a more general MCMC sampler such as Hamiltonian Monte Carlo, have a look at the brms package: