R/compmix.R

In binomlogit: Efficient MCMC for Binomial Logit Models

# Approximation of the type III generalized logistic distribution by Gaussian mixtures
# Computes mixture parameters, given degrees of freedom

compmix=function(k){
  
      # mixture parameters for degrees of freedom <=60
      A <-
      structure(c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 
                16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 
                32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 
                48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 5, 5, 4, 
                4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 
                2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 
                2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 0.0443337, 0.0635862, 
                0.206397, 0.204072, 0.204101, 0.204665, 0.38065, 0.376986, 0.379651, 
                0.379554, 0.401391, 0.401546, 0.382557, 0.384717, 0.381464, 0.380527, 
                0.380767, 0.379757, 0.37807, 0.380199, 0.537661, 0.537773, 0.538551, 
                0.540635, 0.542054, 0.54409, 0.542922, 0.541641, 0.528778, 0.528906, 
                0.53011, 0.528813, 0.527609, 0.52759, 0.52762, 0.526492, 0.526431, 
                0.525753, 0.52538, 0.525902, 0.524433, 0.523196, 0.524131, 0.525114, 
                0.525819, 0.526319, 0.526935, 0.526284, 0.526893, 0.52706, 0.527385, 
                0.526799, 0.525431, 0.521378, 0.522114, 0.521635, 0.519731, 0.514642, 
                0.514756, 0.514866, 0.294977, 0.450637, 0.518711, 0.522425, 0.52318, 
                0.522745, 0.550368, 0.554703, 0.547164, 0.546832, 0.543566, 0.535463, 
                0.537188, 0.541301, 0.547412, 0.548631, 0.550416, 0.549739, 0.547679, 
                0.549532, 0.462339, 0.462227, 0.461449, 0.459365, 0.457946, 0.45591, 
                0.457078, 0.458359, 0.471222, 0.471094, 0.46989, 0.471187, 0.472391, 
                0.47241, 0.47238, 0.473508, 0.473569, 0.474247, 0.47462, 0.474098, 
                0.475567, 0.476804, 0.475869, 0.474886, 0.474181, 0.473681, 0.473065, 
                0.473716, 0.473107, 0.47294, 0.472615, 0.473201, 0.474569, 0.478622, 
                0.477886, 0.478365, 0.480269, 0.485358, 0.485244, 0.485134, 0.429806, 
                0.405544, 0.251858, 0.251401, 0.250535, 0.250057, 0.0689821, 
                0.0683109, 0.0731848, 0.0736144, 0.055043, 0.0629911, 0.0802551, 
                0.0739824, 0.0711245, 0.0708426, 0.068817, 0.0705042, 0.0742508, 
                0.0702688, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                0, 0, 0.207597, 0.0769788, 0.0230344, 0.0221022, 0.022184, 0.0225329, 
                0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0232861, 0.00325404, 0, 
                0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.241146, 0.393076, 
                0.567173, 0.614533, 0.649475, 0.676405, 0.754877, 0.768741, 0.782233, 
                0.792839, 0.806967, 0.815231, 0.818203, 0.824808, 0.829623, 0.834588, 
                0.83931, 0.843504, 0.847291, 0.851328, 0.882606, 0.885332, 0.888014, 
                0.890789, 0.893276, 0.895743, 0.897479, 0.899096, 0.898494, 0.90023, 
                0.902069, 0.903384, 0.90465, 0.90606, 0.907417, 0.908524, 0.909759, 
                0.910843, 0.911936, 0.913126, 0.913956, 0.914785, 0.915936, 0.91705, 
                0.918084, 0.919055, 0.920013, 0.920756, 0.921658, 0.92247, 0.923282, 
                0.923941, 0.924469, 0.924611, 0.925396, 0.925999, 0.926381, 0.926298, 
                0.926942, 0.92757, 0.470351, 0.687823, 0.911262, 0.930135, 0.942488, 
                0.950939, 1.09877, 1.09254, 1.08786, 1.08421, 1.09668, 1.0905, 
                1.07496, 1.07543, 1.07309, 1.07094, 1.06982, 1.0674, 1.0644, 
                1.06465, 1.13647, 1.13337, 1.13066, 1.12851, 1.12631, 1.12441, 
                1.12177, 1.11923, 1.11387, 1.11198, 1.11046, 1.10841, 1.10647, 
                1.10489, 1.10339, 1.10169, 1.1003, 1.09882, 1.09747, 1.09635, 
                1.09487, 1.09349, 1.09258, 1.09171, 1.09083, 1.08993, 1.08909, 
                1.08803, 1.08724, 1.0864, 1.0856, 1.08467, 1.08362, 1.08211, 
                1.0815, 1.08069, 1.07966, 1.07814, 1.07749, 1.07686, 0.91554, 
                1.2185, 1.41674, 1.36085, 1.32143, 1.29187, 1.56324, 1.52393, 
                1.47202, 1.44192, 1.45214, 1.40792, 1.3644, 1.35874, 1.35096, 
                1.33882, 1.33042, 1.31721, 1.30237, 1.29866, 0, 0, 0, 0, 0, 0, 
                0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.80022, 2.05516, 2.31566, 
                2.10397, 1.94996, 1.83732, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                3.57791, 3.95253, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                0), .Dim = c(60L, 12L), .Dimnames = list(NULL, c("V1", "V2",
                "V3", "V4", "V5", "V6", "V7", "V8", "V9", "V10", "V11", "V12"
                )))
      
      coeff_v <-
      structure(c(-4.97e-08, 0.024036101, 1, 0.024403982, 1.165357272, 
                -4.28e-07, 0.027883529, 1, 0.027266081, 0.84360355), .Dim = c(5L,2L),
                .Dimnames = list(NULL, c("V1", "V2")))
      
      K=A[,1]             # degrees of freedom
      NC=A[,2]            # number of components    
      
      # internal function polyval
      polyval=function(a,x) sum(a*x^((length(a)-1):0))
      
      # internal function for rational approximation
      ratval=function(p,x){
        m=which(p==1)
        a=p[1:m]
        b=p[(m+1):length(p)]
        y=polyval(a,x)/polyval(b,x)
        return(y)
      }
  
      # ----------------------- start computing mixture parameters -----------------------
      
      if(k<1) print("Argument must be at least 1")    # error
      
      if(k>=1 & k<=60){                               # look up table
        nc=NC[k]
        mixcomp=data.frame(probs=A[k,3:(3+nc-1)],var=A[k,8:(8+nc-1)])
      }
      
      if(k>60 & k<=500){                              # rational approximation to the parameters
        nc=2
        v=numeric(2)
        p=numeric(2)
        for(i in 1:nc){
          v[i]=ratval(coeff_v[,i],k)
        }
        p[2]=(1-v[1])/(v[2]-v[1])
        p[1]=1-p[2]
        mixcomp=data.frame(probs=p,var=v)
      }
      
      if(k>500){                                      # approximation by standard normal distribution
        nc=1
        mixcomp=data.frame(probs=1,var=1)
      }
      
      # compute non-standardized components
      mixcomp$var=mixcomp$var*2*psigamma(k,1)
      
      # return data.frame
      return(mixcomp)
}