skillTABLE <- function(n11,n01,n10,n00,x,theta=0.5,CI=FALSE,t=1,u=0) {
    # No error checking yet
    #
    # n_ij is the 2x2 table of counts
    # theta in (0,1) is the asymmetric loss criteria
    # CI = TRUE gives 95% confidnece intervals
    # t and u are measurement error parameters; see Briggs Pocernich and Ruppert
    # Either n_ij can be a vector or theta can be;
    # If n_ij a vector, it gives skill scores at each value;
    # If theta is a vector, it gives skill scores for fixed n11 etc. at each theta
    # Both n_ij and theta can NOT be vectors.
    # briggs mattstat@gmail.com

    n<-n11+n01+n10+n00;

    # to determine which of RARE of COMM(ON) is used for each set
    z<-((n10+n11)/n);
    z<-(z<=theta)*1;

        q11<-n11/(n11+n01);
        p11<-(q11-u)/(t-u);
        q00<-n00/(n10+n00);
        p00<-(q00-(1-t))/(t-u);
        px<-(n11+n01)/n
        py<-(n11+n10-n*u)/(n*(t-u))
    k_RARE <- (p11-theta)*px/((1-theta)*py);
    G_RARE <- 2*n11*log(p11/theta) + 2*n01*log((1-p11)/(1-theta));
    k_COMM <- (p00-1+theta)*(1-px)/(theta*(1-py));
    G_COMM <- 2*n00*log(p00/(1-theta))  +2*n10*log((1-p00)/(theta));
    k<-k_RARE*z+k_COMM*(1-z);
    G<-G_RARE*z+G_COMM*(1-z);
    # if k<0 then G=0 and p=.5
    z<-(k>0)*1;
    G<-G*z;
    p<-pchisq(G,1,lower.tail=FALSE)/2;
    p<-p*z+0.5*(1-z);
    bigP<-p11*z+p00*(1-z)
    if(CI){
        ciLO<-0;ciHI<-0;
        rootLO<-0;rootHI<-0;
        f_RARE<-function(p11,n11,n01,theta)
            2*n11*log(p11/theta) + 2*n01*log((1-p11)/(1-theta))- (2*n11*log((n11/(n11+n01))/theta) + 2*n01*log((n01/(n11+n01))/(1-theta))-3.84);
                # extremely rotten way to test for skill; this is a chi^2-like sig value
                # 3.84 for climate
                # 4.25 for markov
        f_COMM<-function(p00,n00,n10,theta)
            2*n00*log(p00/(1-theta)) + 2*n10*log((1-p00)/theta)- (2*n00*log((n00/(n00+n10))/(1-theta)) + 2*n10*log((n10/(n00+n10))/theta)-3.84);
        if(length(theta)>1){
            # program can act funny here because of finding roots is hard!
            for (i in 1:length(theta)){
                    zz<-((n10[i]+n11[i])/n[i]);
                    zz<-(zz<=theta[i]);            
                    tol<-0.0001
                    if(zz){
                        hi<-p11[i]-tol
                        lo<-p11[i]+tol; # print(c(0,hi,lo,n11[i],n01[i],theta[i]))
                        if (is.na(hi) | is.na(lo)){
                            ciLO[i] <- NA;
                            ciHI[i] <- NA;
                        } else {
                            rootLO[i]<-uniroot(f_RARE,lower=0.001,upper=hi,n11=(n11[i]+tol),n01=(n01[i]+tol),theta=theta[i])$root
                            ciLO[i] <- (rootLO[i]-theta[i])*px[i]/((1-theta[i])*py[i]);
                            if(n01[i]==0) {
                                rootHI[i]<-1
                            } else {
                                rootHI[i]<-uniroot(f_RARE,lower=lo,upper=(1-tol),n11=(n11[i]+tol),n01=(n01[i]+tol),theta=theta[i])$root
                            }
                            ciHI[i] <- (rootHI[i]-theta[i])*px[i]/((1-theta[i])*py[i]); 
                        }
                    } else {
                        hi<-p00[i]-tol
                        lo<-p00[i]+tol;  #print(c(1,hi,lo,n00[i],n10[i],theta[i]))
                        if (is.na(hi) | is.na(lo)){
                            ciLO[i] <- NA;
                            ciHI[i] <- NA;
                        } else {
                            rootLO[i]<-uniroot(f_COMM,lower=0.001,upper=hi,n00=(n00[i]+tol),n10=(n10[i]+tol),theta=theta[i])$root    
                            ciLO[i] <- (rootLO[i]-1+theta[i])*(1-px[i])/(theta[i]*(1-py[i]));
                            rootHI[i]<-uniroot(f_COMM,lower=lo,upper=(1-tol),n00=(n00[i]+tol),n10=(n10[i]+tol),theta=theta[i])$root    
                            ciHI[i] <- (rootHI[i]-1+theta[i])*(1-px[i])/(theta[i]*(1-py[i]));
                        }
                    }
            }
        } else {
            for (i in 1:length(n11)){
                    zz<-((n10[i]+n11[i])/n[i]);
                    zz<-(zz<=theta);            
                    tol<-0.0001
                    if(zz){
                        hi<-p11[i]-tol
                        lo<-p11[i]+tol; # print(c(0,hi,lo,n11[i],n01[i],theta))
                        if ((is.na(hi) | is.na(lo)) | (hi<0 | lo<0)){
                            ciLO[i] <- NA;
                            ciHI[i] <- NA;
                        } else {
                            #print(hi);
                            rootLO[i]<-uniroot(f_RARE,lower=0.001,upper=hi,n11=(n11[i]+tol),n01=(n01[i]+tol),theta=theta)$root
                            ciLO[i] <- (rootLO[i]-theta)*px[i]/((1-theta)*py[i]);
                            if(n01[i]==0) {
                                rootHI[i]<-1
                            } else {
                                rootHI[i]<-uniroot(f_RARE,lower=lo,upper=(1-tol),n11=(n11[i]+tol),n01=(n01[i]+tol),theta=theta)$root
                            }
                            ciHI[i] <- (rootHI[i]-theta)*px[i]/((1-theta)*py[i]);                        
                        }
                    } else {
                        hi<-p00[i]-tol
                        lo<-p00[i]+tol; #print(c(1,hi,lo,n00[i],n10[i],p00[i],theta))
                        if ((is.na(hi) | is.na(lo)) | (hi<0 | lo<0)){
                            ciLO[i] <- NA;                        
                            ciHI[i] <- NA;                        
                        } else {
                            rootLO[i]<-uniroot(f_COMM,lower=0.001,upper=hi,n00=(n00[i]+tol),n10=(n10[i]+tol),theta=theta)$root    
                            ciLO[i] <- (rootLO[i]-1+theta)*(1-px[i])/(theta*(1-py[i]));                        
                            rootHI[i]<-uniroot(f_COMM,lower=lo,upper=(1-tol),n00=(n00[i]+tol),n10=(n10[i]+tol),theta=theta)$root    
                            ciHI[i] <- (rootHI[i]-1+theta)*(1-px[i])/(theta*(1-py[i]));                                                
                        }
                    }
            }
        
        
        }
    } else {
        ciLO <- CI
        ciHI <- CI
    }
    imx = which(k==max(k,na.rm=T))
    #print(imx)
    cat("Maximum skill value(s) = ", k[imx],"\n")
    cat("Corresponding x value(s) = ", x[imx],"\n")

    answer<-list(z=z,n=n,k=k,x=x,G=G,p=p,theta=theta,ciLO=ciLO,ciHI=ciHI);
}

skill <- function(x,y,theta=0.5,CI=FALSE,t=1,u=0) {
    # briggs mattstat@gmail.com 
    #if((any(x>1)) | (any(x<0))){
        n<-length(x);
        # x needs to be sorted
        w<-sort(x,index.return=TRUE)
        x<-w$x
        y<-as.numeric(y[w$ix])

        n11<-0;n01<-0;n10<-0;n00<-0;
        for (i in 1:n){
            n11[i]<-0;n01[i]<-0;n10[i]<-0;n00[i]<-0;
            n11[i]<-sum(((x>x[i])*1==1 & y==1),na.rm=T);
            n01[i]<-sum(((x>x[i])*1==1 & y==0),na.rm=T);
            n10[i]<-sum(((x>x[i])*1==0 & y==1),na.rm=T);
            n00[i]<-sum(((x>x[i])*1==0 & y==0),na.rm=T);
        }
        answer<-skillTABLE(n11,n01,n10,n00,x,theta,CI,t=1,u=0);    
    # old test for drawing curves for fixed tables; keep for fun
    #} else {
    #    n<-length(theta);
    #    n11<-0;n01<-0;n10<-0;n00<-0;
    #    for (i in 1:n){
    #        fit<-table((x>theta[i])*1,y)
    #        n11[i]<-0;n01[i]<-0;n10[i]<-0;n00[i]<-0;
    #        if(fit[2,2]) n11[i]<-fit[2,2];
    #        if(fit[2,1]) n01[i]<-fit[2,1];
    #        if(fit[1,2]) n10[i]<-fit[1,2];
    #        if(fit[1,1]) n00[i]<-fit[1,1];
    #    }
    #    answer<-skillTABLE(n11,n01,n10,n00,theta,CI,t=1,u=0);
    #}
}


skillFIT <- function(x,y,theta=0.5,dist="normal"){
    p<-mean(y,na.rm=TRUE);
    Ip<-I(p<=theta);
    if(Ip){
        fity<-fitdistr(x[y==1],dist)
        fitA<-fitdistr(x,dist)
        K<-0*x;
        if(dist=="normal"){
            K<-(1/(p*(1-theta)))*(p*pnorm(x,mean=fity$estimate[1],sd=fity$estimate[2],lower.tail=FALSE) 
             - theta*pnorm(x,mean=fitA$estimate[1],sd=fitA$estimate[2],lower.tail=FALSE)) 
        }
        if(dist=="t"){
            K<-(1/(p*(1-theta)))*(p*pt(x,df=fity$estimate[3],lower.tail=FALSE)
                - theta*pt(x,df=fitA$estimate[3],lower.tail=FALSE))
        }
        if (dist=="gamma") {
            K<-(1/(p*(1-theta)))*(p*pgamma(x,shape=fity$estimate[1],rate=fity$estimate[2],lower.tail=FALSE)
                - theta*pgamma(x,shape=fitA$estimate[1],rate=fitA$estimate[2],lower.tail=FALSE))
        }
        if (dist=="weibull") {
            K<-(1/(p*(1-theta)))*(p*pweibull(x,shape=fity$estimate[1],scale=fity$estimate[2],lower.tail=FALSE)
                - theta*pweibull(x,shape=fitA$estimate[1],scale=fitA$estimate[2],lower.tail=FALSE)) 
        }
        if (dist=="cauchy") {
            K<-(1/(p*(1-theta)))*(p*pcauchy(x,location=fity$estimate[1],scale=fity$estimate[2],lower.tail=FALSE)
                - theta*pcauchy(x,location=fitA$estimate[1],scale=fitA$estimate[2],lower.tail=FALSE))
        }
    } else {
        fity<-fitdistr(x[y==0],dist)
        fitA<-fitdistr(x,dist)
        K<-0*x;
        if(dist=="normal"){
            K<-(1/((1-p)*theta))*((1-p)*pnorm(x,mean=fity$estimate[1],sd=fity$estimate[2])
               - (1-theta)*pnorm(x,mean=fitA$estimate[1],sd=fitA$estimate[2])) 
        }
        if(dist=="t"){
            K<-(1/((1-p)*theta))*((1-p)*pt(x,df=fity$estimate[3])
               - (1-theta)*pt(x,df=fitA$estimate[3]))
        }
        if (dist=="gamma") {
            K<-(1/((1-p)*theta))*((1-p)*pgamma(x,shape=fity$estimate[1],rate=fity$estimate[2])
               - (1-theta)*pgamma(x,shape=fitA$estimate[1],rate=fitA$estimate[2]))
        }
        if (dist=="weibull") {
            K<-(1/((1-p)*theta))*((1-p)*pweibull(x,shape=fity$estimate[1],scale=fity$estimate[2])
               - (1-theta)*pweibull(x,shape=fitA$estimate[1],scale=fitA$estimate[2]))
        }
        if (dist=="cauchy") {
            K<-(1/((1-p)*theta))*((1-p)*pcauchy(x,location=fity$estimate[1],scale=fity$estimate[2])
                - (1-theta)*pcauchy(x,location=fitA$estimate[1],scale=fitA$estimate[2]))
        }
    }
    answer<-list(K=K)
}


skillFIT2 <- function(x,y,theta=0.5,dist="normal"){
    p<-mean(y,na.rm=TRUE);
    Ip<-I(p<=theta);
    if(Ip){
        fity<-fitdistr(x[y==1],dist)
        fitA<-fitdistr(x[y==0],dist)
        K<-0*x;
        if(dist=="normal"){
            K<-(1/(p*(1-theta)))*((1-theta)*pnorm(x,mean=fity$estimate[1],sd=fity$estimate[2],lower.tail=FALSE)*p-
			theta*pnorm(x,mean=fitA$estimate[1],sd=fitA$estimate[2],lower.tail=FALSE)*(1-p))
        }
        if(dist=="t"){
            K<-(1/(p*(1-theta)))*((1-theta)*pt(x,df=fity$estimate[3],lower.tail=FALSE)*p-
			theta*pt(x,df=fitA$estimate[3],lower.tail=FALSE)*(1-p))
        }
        if (dist=="gamma") {
            K<-(1/(p*(1-theta)))*((1-theta)*pgamma(x,shape=fity$estimate[1],rate=fity$estimate[2],lower.tail=FALSE)*p-
			theta*pgamma(x,shape=fitA$estimate[1],rate=fitA$estimate[2],lower.tail=FALSE)*(1-p))
        }
        if (dist=="weibull") {
            K<-(1/(p*(1-theta)))*((1-theta)*pweibull(x,shape=fity$estimate[1],scale=fity$estimate[2],lower.tail=FALSE)*p-
			theta*pweibull(x,shape=fitA$estimate[1],scale=fitA$estimate[2],lower.tail=FALSE)*(1-p))
        }
        if (dist=="cauchy") {
            K<-(1/(p*(1-theta)))*((1-theta)*pcauchy(x,location=fity$estimate[1],scale=fity$estimate[2],lower.tail=FALSE)*p-
			theta*pcauchy(x,location=fitA$estimate[1],scale=fitA$estimate[2],lower.tail=FALSE)*(1-p))
        }
    } else {
        fity<-fitdistr(x[y==1],dist)
        fitA<-fitdistr(x[y==0],dist)
        K<-0*x;
        if(dist=="normal"){
            K<-(1/((1-p)*(1-theta)))*(-(1-theta)*pnorm(x,mean=fity$estimate[1],sd=fity$estimate[2],lower.tail=FALSE)*p+
			theta*pnorm(x,mean=fitA$estimate[1],sd=fitA$estimate[2],lower.tail=FALSE)*(1-p))
        }
        if(dist=="t"){
            K<-(1/((1-p)*(1-theta)))*(-(1-theta)*pt(x,df=fity$estimate[3])*p+
			theta*pt(x,df=fitA$estimate[3])*(1-p))
        }
        if (dist=="gamma") {
            K<-(1/((1-p)*(1-theta)))*(-(1-theta)*pgamma(x,shape=fity$estimate[1],rate=fity$estimate[2])*p+
			theta*pgamma(x,shape=fitA$estimate[1],rate=fitA$estimate[2])*(1-p))
        }
        if (dist=="weibull") {
            K<-(1/((1-p)*(1-theta)))*(-(1-theta)*pweibull(x,shape=fity$estimate[1],scale=fity$estimate[2])*p+
			theta*pweibull(x,shape=fitA$estimate[1],scale=fitA$estimate[2])*(1-p))
        }
        if (dist=="cauchy") {
            K<-(1/((1-p)*(1-theta)))*(-(1-theta)*pcauchy(x,location=fity$estimate[1],scale=fity$estimate[2])*p+
			theta*pcauchy(x,location=fitA$estimate[1],scale=fitA$estimate[2])*(1-p))
        }
    }
    answer<-list(K=K)
}


calibration <- function(x,y,binned=TRUE){
    if(binned){
        x<-round(x*10)/10;
    }
    u<-sort(unique(x))
    p<-0;
        for (i in 1:length(u)){
            j<-(x==u[i])
            p[i]<-sum(y[j]==1,na.rm=T)/length(is.na(y[j]))
        }
    answer<-list(x=u,y=p)
}


#### nonparametric fit functions
rvec<-function(x){cbind(x, x^2)} # transformation of vector x used by Qin and Zhang.

qzfit<-function(t,resp){
   #Perform the semi-parametric fit for each value of t.
   #t corresponds to the covariate x in a traditional logistic regression.
   #v2 is the response y with 0 the non-diseased and 1 the diseased.
   #transform so that x are non-diseased and y are diseased.
   x<-t[resp==0];
   y<-t[resp==1];

   n0<-length(x);
   n1<-length(y);
   rho<-n1/n0;
   #t<-c(x,y);
   #resp<-v2#c(rep(1,n0),rep(1,n1));
   t2<-t^2; t3<-t^3;


   #parameter estimates from Qin and Zhang biometrika 2003.
   stdlogit<-glm(resp~t+t2,family=binomial);
   semi.parm<-optim(stdlogit$coef,betaroot,betagrad,method="Nelder-Mead",t=t,y=y);
   semi.parm2<-optim(semi.parm$par,betaroot,betagrad,method="BFGS",t=t,y=y);

   tillogit<-exp(semi.parm2$par[1]+rvec(t)%*%semi.parm2$par[-1]);

   gtilde<-1/(n0*(1+rho*tillogit));
   ftilde<- tillogit*gtilde;

   return<-cbind(cumsum(ftilde),cumsum(gtilde));
}


npskill<-function(resp,t,theta=.5){
   p<-mean(resp,na.rm=TRUE);
   Ip<-I(p<=theta);
   if(Ip){
      FG<-qzfit(t,resp);
      K<-(1/(p*(1-theta)))*((1-theta)*(1-FG[,1])*p-theta*(1-FG[,2])*(1-p));
   }
   else {
      FG<-qzfit(t,resp);
      K<-(1/((1-p)*(1-theta)))*(-(1-theta)*(1-FG[,1])*p+theta*(1-FG[,2])*(1-p))
   }
   answer<-K;
}

betaroot<-function(parm,t,y){
   n1<-length(y);
   n<-length(t);
   rho=n1/(n-n1);
   lgit<-rho*exp(parm[1]+rvec(t)%*%parm[-1]);
   alpha<-n1-sum(lgit/(1+lgit));
   tmat<-rvec(t);
   sz<-dim(tmat);
   lgitmat<-(lgit/(1+lgit))%*%matrix(1,1,sz[2]);
   beta<-apply(rvec(y),2,sum) - apply(lgitmat*tmat,2,sum);
   answer<-.5*(alpha^2+t(beta)%*%beta);
}


betagrad<-function(parm,t,y){
   n1<-length(y);
   n<-length(t);
   rho=n1/(n-n1);
   lgit<-rho*exp(parm[1]+rvec(t)%*%parm[-1]);
   tmat<-rvec(t);
   sz<-dim(tmat);
   lgitmat<-(lgit/(1+lgit))%*%matrix(1,1,sz[2]);
   lgitmat2<--apply(((lgit/(1+lgit)^2)%*%matrix(1,1,sz[2]))*tmat,2,sum);

   alpha<-n1-sum(lgit/(1+lgit));
   beta<-apply(rvec(y),2,sum) - apply(lgitmat*tmat,2,sum);

   alphaderiv<--alpha*sum(lgit/(1+lgit)^2) + t(beta)%*%lgitmat2;
   part1<-t(tmat%*%beta)%*%(((lgit/(1+lgit)^2)%*%matrix(1,1,sz[2]))*tmat);
   betaderiv<-alpha*lgitmat2-part1;
   answer<-cbind(alphaderiv,betaderiv);
}