########################Fengjuan
Xuan######################################
########################here is the data set named
x#######################
x_c(1,0,6,5,5,2,5,0,3,2,0,6,6,0,4,7,6,3,5,4,5,0,5,7,1,0,6,6,6,2,
4,5,5,0,6,1,7,6,6,5,4,6,4,6,1,2,7,5,5,1,8,3,2,6,2,5,7,5,8,1,
2,3,0,1,0,3,2,0,5,6,5,6,6,4,3,0,6,2,4,3,5,1,1)
#########summarize the data set named y.eg y[1]=11 indicates there
are#####
#########11 x that are equal to 0.and so
on################################
y_1:9
for(i in 0:8){
y[(i+1)]_sum(x==i)
}
##############this very complicated program is written
for some other #course##########
#####it can fit any components mixture model including 2###############################
#####you might ignor the last part ,since that is to calculate the standard errors#####
binomfunc_function(yobs,lambda0,p0,k=length(p0),epsi=10^(-4)){
diff_10
m_length(yobs)-1
oldvalue_logfunc(yobs,lambda0,p0)
oldoldvalue_2*oldvalue
z_matrix(0,(m+1),k)
lambda_lambda0
p_p0
g_0
while(diff>epsi|g<10
){
for(i in
0:m){
z[(i+1),]_lambda*p^i*(1-p)^(m-i)/sum(lambda*p^i*(1-p)^(m-i))
}
temp_t(z)%*%matrix(yobs,(m+1),1)
templambda_apply(temp,1,mean)/sum(yobs)
temp1_t(z)%*%matrix(yobs*(0:m),(m+1),1)
tempp_temp1/templambda/m/sum(yobs)
lambda_templambda/sum(templambda)
p_tempp
newvalue_logfunc(yobs,lambda,p,m)
rate_(newvalue-oldvalue)/(oldvalue-oldoldvalue)
diff_(oldvalue-oldoldvalue)/(1-rate)
oldoldvalue_oldvalue
oldvalue_newvalue
if(g/100==floor(g/100)){
cat("(diff,rate)=",c(diff,rate),"\n")
cat("lambda=",lambda,"\n")
cat("p=",p,"\n")}
g_g+1
}
MD_matrix(0,(m+1),1)
for(j in 1:k){
MD_MD+lambda[j]*dbinom((0:m),m,p[j])
}
expected_sum(yobs)*MD
##calculate
the information matrix##
Info_matrix(0,(2*k-1),(2*k-1))
for(i in
1:(2*k-1)){
for(j in
1:(2*k-1)){
if(i<=(k-1)&
j<=(k-1)& i<=j){
Info[i,j]_-sum((dbinom((0:m),m,p[i])-dbinom((0:m),m,p[k]))*
(dbinom((0:m),m,p[j])-dbinom((0:m),m,p[k]))/MD^2*yobs)
}
if(i>(k-1)&
j>(k-1) & i<=j){
tl_lambda[(i-k+1)]
tp_p[(i-k+1)]
if(i==j){
Info[i,j]_sum(tl*ddfunc(0:m,m,tp)/MD*yobs)-
sum(tl^2*dfunc(0:m,m,tp)^2/MD^2*yobs)
}
else{
Info[i,j]_-sum(tl*lambda[j-k+1]*dfunc(0:m,m,tp)*
dfunc(0:m,m,p[j-k+1])/MD^2*yobs)
}
}
if(j>(k-1)&
i<=(k-1)){
tl_lambda[i]
tp_p[(j-k+1)]
Info[i,j]_-sum((dbinom((0:m),m,p[i])-dbinom((0:m),m,p[k]))*
lambda[j-k+1]*dfunc(0:m,m,tp)/MD^2*yobs)+
sum(dfunc(0:m,m,p[i])/MD*(i==j-k+1)*yobs)
}
}
}
Info_-Info
for(j in
1:(2*k-2)){
for(i in
(j+1):(2*k-1)){
Info[i,j]_Info[j,i]
}
}
bias_expected-yobs
goodness_sum(bias^2/expected)
df_13-2*k
serror_sqrt(diag(solve(Info)))
cat("number
of iteration=",g,"\n")
cat("rate
of convergence=",(1-rate),"\n")
cat("rate
of precision=",sqrt(epsi),"\n")
cat("the
log-likelihood value=",newvalue,"\n")
cat("the chi-square is",goodness,"with
DOF",df,"\n")
list(lambda=lambda,p=p,logvalue=newvalue,bias=bias,SE=serror,
Information=Info,Cov=solve(Info),expected=expected)
}
####################################################################
logfunc_function(yobs,lambda,p,m=length(yobs)-1){
S_0
for(i in 0: m){
S_S+log(sum(lambda*dbinom(i,m,p)))*yobs[i+1]
}
S
}
dfunc_function(x,n,p){
y_x
for(i in
1:length(x)){
y[i]_choose(n,x[i])*p^(x[i]-1)*(1-p)^(n-x[i]-1)*(x[i]-n*p)
}
y
}
ddfunc_function(x,n,p){
y_x
for(i in
1:length(x)){
y[i]_choose(n,x[i])*p^(x[i]-2)*(1-p)^(n-x[i]-2)*
((x[i]-n*p)*(x[i]-1-n*p+2*p)-n*p*(1-p))
}
y
}
####################################################################
###############Now here is the
result###############################
$lambda
[1] 0.3645353 0.6354647
$p
[,1]
[1,] 0.1351944
[2,] 0.6547623
$logvalue
[1] -173.7489
$bias
[,1]
[1,] -1.5233766
[2,]
2.9999770
[3,] -1.4505243
[4,] -0.9084670
[5,]
3.0357103
[6,] -1.3243396
[7,] -3.1264674
[8,]
2.5157670
[9,] -0.2182794
$SE
[1] 0.06420369 0.03085614 0.02945034
$Information
[,1] [,2]
[,3]
[1,] 300.3955 -183.6139 -151.4544
[2,] -183.6139 1347.9936 -409.0000
[3,] -151.4544 -409.0000 1432.9154
$Cov
[,1] [,2] [,3]
[1,] 0.0041221137 0.0007594514 0.0006524655
[2,] 0.0007594514 0.0009521016 0.0003520318
[3,] 0.0006524655 0.0003520318 0.0008673225