#This function takes lambda, data set, pos (market positions), #and locus (number of possible positions of x_s) #Generate a sequence, lsupp, between .001 and .999 of size nlocus. #To get simple estimators, use lambda= lambda hat obtained in simple.est.null #Fix omegahat = max (LI/I) #For fixed x_s, maximize log-likelihood function over delta, using nlminb. #Store log-likelihood into a vector (loglik) at each x_s #Return omegahat,lambdahat, deltahat, lsupp, and loglik simple.est.alt<-function(lambda,data,pos,nlocus){ LI<-apply(data,2,function(x){sum(x==1)}) I<-apply(data,2,function(x){sum(x!=-1)}) I[I==0]<-Inf #to prevent from dividing by 0 omegahat<- max(LI/I) lsupp<-seq(.001, .999, length=nlocus) loglik<-c(NA,nlocus) deltahat<-c(NA,nlocus) message<-c(NA,nlocus) for(i in 1:nlocus) { # print(i) if(.001<(omegahat/2)) { profile<-nlminb(start=(omegahat/2),mloglik.alt,omega=omegahat, leta=log(1/lambda),locus=lsupp[i], data=data,pos=pos,lower=0.001,upper=omegahat) deltahat[i]<-profile$parameters loglik[i]<-(-profile$objective) message[i]<-profile$message } else { deltahat[i]<-NA loglik[i]<-NA message[i]<-NA } } loglik.alt<-max(loglik) index<-(loglik==loglik.alt) deltahat<-deltahat[index] locushat<-lsupp[index] return(list(omegahat=omegahat,lambdahat=lambda,deltahat=deltahat, locushat=locushat,loglik.alt=loglik.alt, loglik=loglik,lsupp=lsupp,message=message)) } ######## Example ###################################### pro<-dget("pro.data") thy<-dget("thy.data") fjp<-dget("fjp.data") pro.simple.est.alt<-simple.est.alt(lambda= pro.simple.est.null$lambdahat, data=pro$data, pos=pro$pos, nlocus=25) thy.simple.est.alt<-simple.est.alt(lambda= thy.simple.est.null$lambdahat, data=thy$data, pos=thy$pos, nlocus=25) fjp.simple.est.alt<-simple.est.alt(lambda= fjp.simple.est.null$lambdahat, data=pro$data, pos=fjp$pos, nlocus=25) ########################################################