########input argument#########
###x is a vector containing the expression intensities of a gene. The noncancer samples 
are listed before the cancer samples (or the group to be tested).
###nn is the sample size of the noncancer group.
###nc is the sample size of the cancer group.
###m0 is the low-end of the range of x to be tested.
###m1 is the up-end of the range of x to be tested.

########output argument########
###lrtsq is the likeihood ratio test statistic.
###pvalue is the p-value of the test.

getcp=function(x,nn,nc,m0,m1){
    
  xsort=c(x[1:nn], sort(x[(nn+1):(nn+nc)]))
   
  nqtest=nn+nc      
  qaimidx=m0:m1
  
  sk=cumsum(xsort)[qaimidx]
  lrtsqi=(qaimidx*sum(xsort)/nqtest-sk)/(qaimidx*(1-qaimidx/nqtest))^0.5  

  lrtsq=max(lrtsqi)   

 if (lrtsq==-Inf) lrtsq==0
 
 if (lrtsq<=0){
    lrtsq=0;xqchangepoint=0; xqcpprop=0;pvalue=1}
  
 if (lrtsq>0){

  xqcpidx=qaimidx[lrtsqi==lrtsq]   
  xqcpprop=(m1-xqcpidx+1)/nc 
  b=lrtsq 
 
  gamma=b/nqtest^0.5
  
  intf1=function(x){
     output=(1-x^2)^((nqtest-4)/2)
     return(output)
  }

 if (gamma<1){
  int1=integrate(intf1,lower=gamma,upper=1)$value 

  intf2=function(x){
      xnew=x+b^2/(nqtest*(1-gamma^2)*x)
      nuxnew=exp(-0.583*xnew)
      output=1/x*nuxnew
      return(output)
  }

  lower2=b*((1/(qaimidx[length(qaimidx)])-1/nqtest)/(1-gamma^2))^0.5  
  upper2=b*((1/(qaimidx[1])-1/nqtest)/(1-gamma^2))^0.5  

  int2=integrate(intf2,lower=lower2,upper=upper2)$value
  pvalue=(nqtest/(2*pi))^0.5*int1+(2*pi)^(-0.5)*b*(1-b^2/nqtest)^((nqtest-4)/2)*int2
 }

 if (gamma>=1) pvalue=NA
 }

 return(list(lrtsq=lrtsq,pvalue=pvalue))
}