Aug. 2012


The original code was written in 1992, and runs on Splus on an HP unix workstation.

We have since stopped paying for Splus lisence and using R nowadays.

I just dig up the code and did minimal changes and then 
try the code on R version 2.14.1. And it runs. 
The original code has two versions of the isot( ) function; one runs in C code 
and thus a lot faster, but making the C to work in R will take more effort. 

So this is just the R version below (and a bit slower).  :-(




Example run:   The regression model is y= x1 + x2 + epsilon

 set.seed(1234)
 x <- matrix(c(rnorm(250), rnorm(250, mean=1)), ncol=2, byrow=F)
 z <- x[,1] + x[,2] + rnorm(250)
 delta <- rep(1,250)
 delta[z<0] <- -1 

> bibj(x=x, delta=delta, error=0.000001)
[1] 1.0000000 0.8107448 9.0000000

So, in 9 interations the estimator is 1.0 and 0.8107 (true value is 1)
because of the standardization, the beta1 is always taking as 1


Aug.20, 2012

I made a small change so that it also output the estimated F( ) distribution of epsilon.


=================================================
Now the code itself

bibj <- function(x, delta, maxiter = 20, error = 0.0001) {  
### x is a design matrix of N rows, delta is a vector of length N.
### delta =1 mean the actual obs (y) is > 0. delta = -1 means y < 0. 
### Output: beta is a vector of length = no. of column(x)
 
dimofX <- dim(x)
if( dimofX[1] != length(delta) ) stop("check dim of x and /or delta")

      newtemp <- matrix(NA, ncol=dimofX[2], nrow=4)

      newtemp[1,]<-lsfit(x,delta,intercept=F)$coef 
      newtemp[1,]<-newtemp[1,]/newtemp[1,1]


      for (i in 2:4) { 
      newtemp[i,]<-iter(x,delta,newtemp[i-1,])
      newtemp[i,]<-newtemp[i,]/newtemp[i,1]
      } 

      num<-3 
      while(num <= maxiter && error <sum((newtemp[3,]-newtemp[4,])^2) ) {
             newtemp[3,] <- newtemp[4,]
             newtemp[4,]<-iter(x,delta,newtemp[3,] )
             newtemp[4,]<-newtemp[4,]/newtemp[4,1]
             num <- num + 1 
      } 
      delta[delta<0.5] <- 0
      u <- - x %*% newtemp[4,]
      tempF <- isot(x=u, y=1-delta, a=2, b=2)
      list( est=newtemp[4,], iterN=num,  distFx=tempF$x, distFy=tempF$y )
}




iter <- function(x,delta,beta) {  
### x is a matrix of N rows, delta is a vector of length N.
### delta =1 mean the actual obs (y) is > 0.
### beta is a vector of length = no. of column(x)

N <- length(delta)
k <- length(as.vector(beta))

if (dim(x)[1] != N) stop("check dim of x")
if (dim(x)[2] != k) stop("check length of beta and dim of x") 
u <- - x %*% beta 
ystar <- - u      # should I just let  ystar <- delta ?

dd <- delta
dd[delta<0.5] <- 0 
temp<-isot(u,1-dd,2,2)
tx<-temp$x
ty<-temp$y
TN<-length(tx)

    for (i in 1:N) {

               m<-length(tx[tx <= u[i]])
               if(delta[i] > 0.5 && ty[m] < 1 ) {
                             ystar[i]<- - u[i] + 0.5*sum((tx[(m+1):TN]+tx[m:(TN-1)])*diff(ty[m:TN])) / 
                                     (1- ty[m] ) 
                                }
              else if( ty[m] > 0 ) {
                    ystar[i]<- - u[i] + 0.5*sum((tx[2:m]+tx[1:(m-1)])*diff(ty[1:m])) / ty[m]
                    }

    }

return( lsfit(x,ystar,intercept=FALSE)$coef ) 

}




isot <- function(x,y,a,b) {  
### the parameter a, b are to argument output F
### so that F(x[1]-a)=0 and F(x[n]+b) =1 always 
### y are 0's and 1's.
if ( length(x) != length(y) ) stop("check length of x and/or y")
 
xsort <- unique(sort(x))
m<-length(xsort)
ysort <- rep(NA,m)
nsort <- rep(NA,m)
for ( i in 1:m ) {
ysort[i] <- sum(y[x == xsort[i]])
nsort[i] <- sum(x == xsort[i])
}
y <- ysort
n <- nsort

ms <- 1
ys <- y[1]
ns <- n[1]

for ( i in 2:length(y) ) {

ys <- c(ys,y[i])
ns <- c(ns,n[i])
ms <- c(ms,1)
ps <- ys/ns

j <- length(ys)
while ( length(ps) > 1 && ps[j-1] > ps[j] ) {
ys[j-1] <- ys[j-1] + ys[j] ; ys <- ys[-j]
ns[j-1] <- ns[j-1] + ns[j] ; ns <- ns[-j]
ms[j-1] <- ms[j-1] + ms[j] ; ms <- ms[-j]
ps <- ys/ns
j <- j - 1
}
}

return( list(x=c(xsort[1]-a,xsort,xsort[m]+b), y=c(0,rep(ps,ms),1)) )

}


A final note: The main bottle neck of the running time is the isot( ) function. 
Call a C function equivalent inside R to speed it up if you can. 
I do have a C function from 1992. Just need to search/find it, if there is sufficient interest.
(minor changes needed to make it call-able from within R).

Mai Zhou