code: def estimateHeritabilities(self, K, verbose=False):
          
          
          # Fit single trait model
          varg  = SP.zeros(self.P)
          varn  = SP.zeros(self.P)
          fixed = SP.zeros((1,self.P))
          
          for p in range(self.P):
              y = self.Y[:,p:p+1]
              lmm = limix.CLMM()
              lmm.setK(K)
              lmm.setSNPs(SP.ones((K.shape[0],1)))
              lmm.setPheno(y)
              lmm.setCovs(SP.zeros((K.shape[0],1)))
              lmm.setVarcompApprox0(-20, 20, 1000)
              lmm.process()
              delta = SP.exp(lmm.getLdelta0()[0,0])
              Vtot  = SP.exp(lmm.getLSigma()[0,0])
      
              varg[p] = Vtot
              varn[p] = delta*Vtot
              fixed[:,p] = lmm.getBetaSNP()
      
              if verbose: print(p)
      
          sth = {}
          sth['varg']  = varg
          sth['varn']  = varn
          sth['fixed'] = fixed
  
          return sth
docs: estimate variance components and fixed effects
          from a single trait model having only two terms