If we do a negative binomial GLM with unknown deviance, a frequent strategy (used for example by glm.nb in package MASS in R) is to use a Gibbs sampler:
Isn't the mean the same than that would be guessed by Poisson regression? If so, why bother with the Gibbs sampler? Why not just fit the mean using Poisson, then fit the dispersion?
It's not the same but it can be close. In this example, offset and exponent are less than $10^{-3}$ apart between the poisson and negative binomial estimate. The Poisson dispersion estimate is 1.08, while the negative binomial is 1.09 
R code:
library(ggplot2)
library(data.table)
library(MASS)
data=data.table(x=1:100)
data[,mu:=5*x]
data[,y:=rnbinom(100, mu = mu, size=1)]
lmp=glm(y~log(x),data=data,family=poisson(link="log"))
data[,pred.poisson:=predict(lmp,type="response")]
disp.pois=data[,theta.ml(y,mu,.N)]
lmnb=glm.nb(y~log(x), data=data, link="log")
data[,pred.nb:=predict(lmnb,type="response")]
ggplot(data)+geom_point(aes(x,y))+geom_line(aes(x,mu,colour="true mean"))+
geom_line(aes(x,pred.poisson,colour="poisson"))+
geom_line(aes(x,pred.nb,colour="nb"))