code: def lpointbiserialr(x,y):
      
      TINY = 1e-30
      if len(x) != len(y):
          raise ValueError('INPUT VALUES NOT PAIRED IN pointbiserialr.  ABORTING.')
      data = pstat.abut(x,y)
      categories = pstat.unique(x)
      if len(categories) != 2:
          raise ValueError("Exactly 2 categories required for pointbiserialr().")
      else:   # there are 2 categories, continue
          codemap = pstat.abut(categories,range(2))
          recoded = pstat.recode(data,codemap,0)
          x = pstat.linexand(data,0,categories[0])
          y = pstat.linexand(data,0,categories[1])
          xmean = mean(pstat.colex(x,1))
          ymean = mean(pstat.colex(y,1))
          n = len(data)
          adjust = math.sqrt((len(x)/float(n))*(len(y)/float(n)))
          rpb = (ymean - xmean)/samplestdev(pstat.colex(data,1))*adjust
          df = n-2
          t = rpb*math.sqrt(df/((1.0-rpb+TINY)*(1.0+rpb+TINY)))
          prob = betai(0.5*df,0.5,df/(df+t*t))  # t already a float
          return rpb, prob
docs: Calculates a point-biserial correlation coefficient and the associated
  probability value.  Taken from Heiman's Basic Statistics for the Behav.
  Sci (1st), p.194.
  
  Usage:   lpointbiserialr(x,y)      where x,y are equal-length lists
  Returns: Point-biserial r, two-tailed p-value