library(ggplot2)
#######################define the parameter
c<-0.4
alpha1<- 0.0415
alpha2<- 0.05-alpha1
##################
########################situation for m=1.02
m<-1.02 #mean
b<-m*(1-c)^2
pi<-(1-b-c)/(c*(1-c))
n.hat<-log((1-pi)*875+pi)/log(m) ######MLE estimator
a1<- -log(1-alpha1)
a2<- -log(alpha2)
#get the confidence interval
cilow<-n.hat-log(a2)/log(m)  ##########low boundary of confidence interval
ci.high<-n.hat-log(a1)/log(m) ########high boundary of confidence interval
########################situation for m=1.05
m<-1.05
b<-m*(1-c)^2
pi<-(1-b-c)/(c*(1-c))
n.hat<-log((1-pi)*875+pi)/log(m)
a1<- -log(1-alpha1)
a2<- -log(alpha2)
#get the confidence interval
cilow<-n.hat-log(a2)/log(m)  ##########low boundary of confidence interval
cihigh<-n.hat-log(a1)/log(m)  ########high boundary of confidence interval
########################situation for m=1.1
m<-1.1
b<-m*(1-c)^2
pi<-(1-b-c)/(c*(1-c))
n.hat<-log((1-pi)*875+pi)/log(m)
a1<- -log(1-alpha1)
a2<- -log(alpha2)
#get the confidence interval
cilow<-n.hat-log(a2)/log(m)   ##########low boundary of confidence interval
cihigh<-n.hat-log(a1)/log(m)  ########high boundary of confidence interval
############################situation for m=1
n.hat<- (1-c)/c*(875-1)
a1<- -log(1-alpha1)
a2<- -log(alpha2)
cilow<- n.hat/a2
cihigh<- n.hat/a1
###################################
m<-seq(101,120,1)/100
n.hat<-numeric(20)
cilow<- numeric(20)
cihigh<- numeric(20)
for (i in 1:20) {
  b<-m[i]*(1-c)^2
  pi<-(1-b-c)/(c*(1-c))
  n.hat[i]<-log((1-pi)*4+pi)/log(m)
  a1<- -log(1-alpha1)
  a2<- -log(alpha2)
  #get the confidence interval
  cilow[i]<-n.hat[i]-log(a2)/log(m[i])
  cihigh[i]<-n.hat[i]-log(a1)/log(m[i])
}

## Warning in n.hat[i] <- log((1 - pi) * 4 + pi)/log(m): number of items to
## replace is not a multiple of replacement length

## Warning in n.hat[i] <- log((1 - pi) * 4 + pi)/log(m): number of items to
## replace is not a multiple of replacement length

## Warning in n.hat[i] <- log((1 - pi) * 4 + pi)/log(m): number of items to
## replace is not a multiple of replacement length

## Warning in n.hat[i] <- log((1 - pi) * 4 + pi)/log(m): number of items to
## replace is not a multiple of replacement length

## Warning in n.hat[i] <- log((1 - pi) * 4 + pi)/log(m): number of items to
## replace is not a multiple of replacement length

## Warning in n.hat[i] <- log((1 - pi) * 4 + pi)/log(m): number of items to
## replace is not a multiple of replacement length

## Warning in n.hat[i] <- log((1 - pi) * 4 + pi)/log(m): number of items to
## replace is not a multiple of replacement length

## Warning in n.hat[i] <- log((1 - pi) * 4 + pi)/log(m): number of items to
## replace is not a multiple of replacement length

## Warning in n.hat[i] <- log((1 - pi) * 4 + pi)/log(m): number of items to
## replace is not a multiple of replacement length

## Warning in n.hat[i] <- log((1 - pi) * 4 + pi)/log(m): number of items to
## replace is not a multiple of replacement length

## Warning in n.hat[i] <- log((1 - pi) * 4 + pi)/log(m): number of items to
## replace is not a multiple of replacement length

## Warning in n.hat[i] <- log((1 - pi) * 4 + pi)/log(m): number of items to
## replace is not a multiple of replacement length

## Warning in n.hat[i] <- log((1 - pi) * 4 + pi)/log(m): number of items to
## replace is not a multiple of replacement length

## Warning in n.hat[i] <- log((1 - pi) * 4 + pi)/log(m): number of items to
## replace is not a multiple of replacement length

## Warning in n.hat[i] <- log((1 - pi) * 4 + pi)/log(m): number of items to
## replace is not a multiple of replacement length

## Warning in n.hat[i] <- log((1 - pi) * 4 + pi)/log(m): number of items to
## replace is not a multiple of replacement length

## Warning in n.hat[i] <- log((1 - pi) * 4 + pi)/log(m): number of items to
## replace is not a multiple of replacement length

## Warning in n.hat[i] <- log((1 - pi) * 4 + pi)/log(m): number of items to
## replace is not a multiple of replacement length

## Warning in n.hat[i] <- log((1 - pi) * 4 + pi)/log(m): number of items to
## replace is not a multiple of replacement length

## Warning in n.hat[i] <- log((1 - pi) * 4 + pi)/log(m): number of items to
## replace is not a multiple of replacement length

##############get the coverage plot
m<-1.02 ##########change this value to get the coverage prabability for m=1.05, m=1.1
b<-m*(1-c)^2
pi<-(1-b-c)/(c*(1-c))
#####################################################
###############define the number of generarion
coverage<-numeric(500)
for (n in 1:500) {
  bn<-m^n*((1-pi)/(m^n-pi))^2
  cn<-(m^n-1)/(m^n-pi)
  pn<-numeric(1000)
  for (i in 1:1000) {
    pn[i]<-(1-cn)*cn^(i-1)
    if(pn[i]<0.000000001)
      break
  }
  pn<-pn[pn>0]
  nopn<-length(pn)##############the number of probability
  result10<- numeric(1000)
  for (times in 1:1000) {
    result10[times]<-sample(1:nopn,1,prob = pn)
  }
##############################
n.hat<-log((1-pi)*result10+pi)/log(m)
a1<- -log(1-alpha1)#define the value of a1
a2<- -log(alpha2)#define the value of a2
#get the confidence interval
  cilow<-n.hat-log(a2)/log(m)
  cihigh<-n.hat-log(a1)/log(m)
  coverage[n]<-sum(cilow<n&cihigh>n)/1000
}
df<-data.frame(x=c(1:length(coverage)),y=coverage)
st<-sqrt(0.95*0.05/1000) #binomial standard error
h<- 0.95+1.96*st   #high boundaruy of confidence interval
l<- 0.95-1.96*st    #low boundary of confidence interval
ggplot()+geom_point(data = df,aes(df[,1],df[,2]))+geom_hline(yintercept = h,lty="dashed")+
  geom_hline(yintercept = l,lty="dashed")+
  ylab("Coverage Probability")+theme(axis.title.y = element_text(size = 15))+
  xlab("Age")+theme(axis.title.x = element_text(size = 15))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

#####the length of the confidence interval
L<-log(a2/a1)/log(m)
L*log(m)

## [1] 4.722806

########## get The plot of age versus mean
m<-1.01
n.hat<-numeric(20)
cilow<- numeric(20)
cihigh<- numeric(20)
for (i in 1:20) {
  b<-m*(1-c)^2
  pi<-(1-b-c)/(c*(1-c))
  n.hat[i]<-log((1-pi)*875+pi)/log(m)
  a1<- -log(1-alpha1)
  a2<- -log(alpha2)
  #get the confidence interval
  cilow[i]<-n.hat[i]-log(a2)/log(m)
  cihigh[i]<-n.hat[i]-log(a1)/log(m)
  m<-1.01+0.01*i
}
mm<-seq(1.01,1.20,0.01)
plot(mm,n.hat)

library(ggplot2)
df<-data.frame(Mean=mm,Age=n.hat,l=cilow,h=cihigh)
sp<-ggplot()+geom_line(data=df,mapping=aes(df$Mean,df$Age))+geom_line(data = df,mapping = aes(Mean,l),lty="dashed")+
  geom_line(data = df,mapping = aes(Mean,h),lty="dashed")+xlim(c(1.01,1.2))
sp + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  xlab("Mean")+ theme(axis.title.x = element_text(size = 15))+
  ylab("Age")+ theme(axis.title.y = element_text(size = 15))

############################few k=4 in the situaion of young variant
c<-0.4
alpha1<- 0.0415
alpha2<- 0.05-alpha1
############m=1.02
m<-1.02
b<-m*(1-c)^2
pi<-(1-b-c)/(c*(1-c))
n.hat<-log((1-pi)*4+pi)/log(m)
a1<- -log(1-alpha1)
a2<- -log(alpha2)
#get the confidence interval
cilow<-n.hat-log(a2)/log(m)
cihigh<-n.hat-log(a1)/log(m)
###########m=1.05
m<-1.05
b<-m*(1-c)^2
pi<-(1-b-c)/(c*(1-c))
n.hat<-log((1-pi)*4+pi)/log(m)
a1<- -log(1-alpha1)
a2<- -log(alpha2)
#get the confidence interval
cilow<-n.hat-log(a2)/log(m)
cihigh<-n.hat-log(a1)/log(m)
#############m=1.1
m<-1.1
b<-m*(1-c)^2
pi<-(1-b-c)/(c*(1-c))
n.hat<-log((1-pi)*4+pi)/log(m)
a1<- -log(1-alpha1)
a2<- -log(alpha2)
#get the confidence interval
cilow<-n.hat-log(a2)/log(m)
cihigh<-n.hat-log(a1)/log(m)