用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)