用sas的nlp过程求两阶段的负二项分布的最大释然估计-经管之家官网！

用sas的nlp过程求两阶段的负二项分布的最大释然估计。数据是模拟出来的。
程序如下：
data part1(type=est);
keep _type_expb0 expb1 theta;
_type_='parms'; expb0 = 0 ; expb1=1 ; theta=0.5; output;
_type_='lb'; expb0 = 1.0e-6 ; expb1 =1.0e-6 ;theta = 1.0e-6;output;
run;
proc nlp data=sasuser.com tech=tr inest=part1 outest=opart1 /*b0=0.05,b1=2,gamma=0.2,lamda=5,n=100*/
outmodel=model1 cov=2 vardef=n pcov phes;
max logf;
parms expb0 expb1 theta;
profile expb0 expb1 theta / alpha = 0.05;
y_th = y; x_th = x;
lamda = 5.3600 ; /* lamda is the mle of the second stage model of possion distribution */
b=log(gamma(y_th+theta))-log(gamma(theta))+theta*log(theta) ;
if x_th=1 then c=lamda*expb0*expb1;
if x_th=0 then c=lamda*expb0;
a=y*log(c)-(y+theta)*log(theta+c);
s=a+b;
logf = s;
run;
最后算出来的Theta不对呢？请各位高手指点迷津！！
数据产生程序如下（R软件）：
simudata<-function(b0,b1,gamma,lamda,n){
t<-rgamma(n,1/gamma,1,gamma)
m<-rpois(n,lamda)
x<-rbinom(n,1,0.5)
u<-c(1:2)
y<-c(1:2)
p<-c(1:2)
for(i in 1:n)
{ b<-x
if(b<=0)
u<-b0
else
u<-b0*b1
}
for(i in 1:n){
p<-t*u
o<-m
b<-rbinom(1,o,p)
y<-b
}
data<-data.frame(Y=y,M=m,X=x)
print(data)
}
simudata(0.05,2,0.2,5,100)