R/mixcomp_poisson.R

In airbornemint/pogit: Bayesian Variable Selection for a Poisson-Logistic Model

## Approximation of the negative log-gamma distribution by Gaussian mixtures.
## Compute parameters of Gaussian mixture, given degrees of freedom v: 
## with R(v) components, weights w(v), means m(v) and variances v(v) 

## This file is based on the bayesf MATLAB package found on S. Fruehwirth-Schnatters'
## website 
## <http://statmath.wu.ac.at/~fruehwirth/monographie/>.

#------ individual mixtures ---------
# R(v)= 10 components for 1<= v <= 4          
# R(v)= 9 components for 5 <= v <= 19         
#----- parameterized mixtures -------
# R(v)= 4 components for 20 <= v <= 49        ... range = 3
# R(v)= 3 componnets for 50 <= v <= 439       ... range = 4
# R(v)= 2 components for 440 <= v <= 1599     ... range = 5
# R(v)= 2 components for 1600 <= v <= 10000   ... range = 6
# R(v)= 2 components for 10000 <= v <= 30000  ... range = 7
#----------------------------------------------------------------------------------------------

#----------------------------------------------------------------------------------------------
# M I X T U R E   C O M P O N E N T S
#----------------------------------------------------------------------------------------------

mixcomp_poisson <- function(){
  w <- matrix(0, 30000, 10)
  m <- matrix(0, 30000, 10)
  v <- matrix(0, 30000, 10)

  #---------------   I N D I V I D U A L    M I X T U R E S  -------------------------------------
  # weights
  w1  <- c(0.00396984425,0.0396244597,0.16776747,0.147036501,0.125306271,0.1014852,0.103758531,0.115972617,0.107065659,0.0880134482)
  w2  <- c(0.00396860865,0.039627049,0.16777003,0.147037759,0.125304523,0.101481714,0.103759705,0.115973128,0.107065554,0.0880119305)
  w3  <- c(0.00396904734,0.0396280642,0.167770514,0.14703607,0.12530043,0.10148242,0.103759287,0.115974323,0.107066971,0.0880128738)
  w4  <- c(0.00396883344,0.039627504,0.167771274,0.147036659,0.125301189,0.101481755,0.103760036,0.115974339,0.107065718,0.0880126919)
  w5  <- c(0.04358203,0.16779435,0.14704072,0.12531065,0.10147112,0.10376347,0.11597388,0.10705620,0.08800758)
  w6  <- c(0.04358170,0.16779599,0.14704260,0.12531102,0.10147067,0.10376308,0.11597286,0.10705544,0.08800665)
  w7  <- c(0.04357986,0.16779709,0.14704207,0.12531329,0.10147098,0.10376185,0.11597323,0.10705435,0.08800728) 
  w8  <- c(0.04357889,0.16779743,0.14704199,0.12531387,0.10147092,0.10376158,0.11597314,0.10705400,0.08800818)
  w9  <- c(0.04357867,0.16779774,0.14704280,0.12531355,0.10147095,0.10376139,0.11597319,0.10705364,0.08800807)
  w10 <- c(0.04357793,0.16779779,0.14704273,0.12531407,0.10147145,0.10376142,0.11597319,0.10705347,0.08800795)
  w11 <- c(0.04357676,0.16780138,0.14704262,0.12531407,0.10147055,0.10376107,0.11597323,0.10705197,0.08800836)
  w12 <- c(0.04357718,0.16780110,0.14704244,0.12531386,0.10147030,0.10376152,0.11597319,0.10705242,0.08800798)
  w13 <- c(0.04357785,0.16780049,0.14704195,0.12531391,0.10147008,0.10376171,0.11597361,0.10705276,0.08800765)
  w14 <- c(0.04357864,0.16779926,0.14704212,0.12531339,0.10147055,0.10376238,0.11597379,0.10705254,0.08800733)
  w15 <- c(0.04358150,0.16779787,0.14704361,0.12531252,0.10146908,0.10376217,0.11597341,0.10705436,0.08800546)
  w16 <- c(0.04358114,0.16779895,0.14704369,0.12531262,0.10146892,0.10376205,0.11597342,0.10705397,0.08800525)
  w17 <- c(0.04358126,0.16779887,0.14704452,0.12531232,0.10146888,0.10376173,0.11597269,0.10705405,0.08800568)
  w18 <- c(0.04358087,0.16779900,0.14704453,0.12531268,0.10146895,0.10376147,0.11597264,0.10705388,0.08800597)
  w19 <- c(0.04358073,0.16779923,0.14704464,0.12531292,0.10146881,0.10376128,0.11597263,0.10705366,0.08800610)
  
  # means
  m1  <- c(5.14164960,3.28872998,1.81888375,1.23821155,0.76378235,0.39080161,0.04282805,-0.30646483,-0.67322634,-1.05559858)
  m2  <- c(1.81625904,1.06274371,0.33888820,0.04113215,-0.24333721,-0.48957362,-0.73371533,-0.98183872,-1.24066737,-1.49981491)
  m3  <- c(0.6071746,0.1653763,-0.3306968,-0.5444110,-0.7646286,-0.9633027,-1.1614760,-1.3606682,-1.5696376,-1.7756108)
  m4  <- c(-0.02821078,-0.34388039,-0.74845607,-0.93586649,-1.12166807,-1.28756741,-1.45622175,-1.62854829,-1.80876206,-1.98114148)
  m5  <- c(-0.8892953,-1.0526376,-1.1959978,-1.3928605,-1.5573236,-1.6247121,-1.8141066,-1.9960747,-2.1547473)
  m6  <- c(-1.169928,-1.298515,-1.422049,-1.598956,-1.759445,-1.808373,-1.982529,-2.148634,-2.291363)
  m7  <- c(-1.396288,-1.501396,-1.609664,-1.768774,-1.913815,-1.969604,-2.128834,-2.280606,-2.409574)
  m8  <- c(-1.584982,-1.672755,-1.768398,-1.908550,-2.056099,-2.109670,-2.254861,-2.395355,-2.513576)
  m9  <- c(-1.746432,-1.820395,-1.905389,-2.033826,-2.179360,-2.234645,-2.364767,-2.497469,-2.606721)
  m10 <- c(-1.885814,-1.949444,-2.027657,-2.148290,-2.287289,-2.350002,-2.459770,-2.588560,-2.690438)
  m11 <- c(-2.011021,-2.066343,-2.138910,-2.256248,-2.389999,-2.425558,-2.554919,-2.672373,-2.767415)
  m12 <- c(-2.122447,-2.172019,-2.237011,-2.345187,-2.474035,-2.528341,-2.632660,-2.749537,-2.838774)
  m13 <- c(-2.224250,-2.267041,-2.327989,-2.432455,-2.561510,-2.595176,-2.713155,-2.821338,-2.904757)
  m14 <- c(-2.315246,-2.353100,-2.412718,-2.514480,-2.644595,-2.664389,-2.782449,-2.885121,-2.964992)
  m15 <- c(-2.406138,-2.439038,-2.489970,-2.583989,-2.705827,-2.741717,-2.842697,-2.946913,-3.023404)
  m16 <- c(-2.480686,-2.511804,-2.563433,-2.655585,-2.771971,-2.802656,-2.907052,-3.004767,-3.077292)
  m17 <- c(-2.558039,-2.586419,-2.632158,-2.718344,-2.835261,-2.859474,-2.958664,-3.058637,-3.129204)
  m18 <- c(-2.624498,-2.649278,-2.694619,-2.780596,-2.892908,-2.921123,-3.015017,-3.109458,-3.176608)
  m19 <- c(-2.685950,-2.710506,-2.755594,-2.839381,-2.946400,-2.973889,-3.069450,-3.157098,-3.221553)
  
  # variances
  v1  <- c(4.31964673,1.99470699,1.09530437,0.42217332,0.19750911,0.10690578,0.07788990,0.07671881,0.09483729,0.14615846)
  v2  <- c(1.54538942,0.75454554,0.44427581,0.19804942,0.10038033,0.05800283,0.04565027,0.04848328,0.06389837,0.09978319)
  v3  <- c(0.85390901,0.43862561,0.26946132,0.12825508,0.06840449,0.04267804,0.03621836,0.03959501,0.05188316,0.07899662)
  v4  <- c(0.54760673,0.28572995,0.18135815,0.09163303,0.05235144,0.03468812,0.03019878,0.03318938,0.04374090,0.06627978)
  v5  <- c(0.34820369,0.16490241,0.07536732,0.04566448,0.04580879,0.02944360,0.02743986,0.03439750,0.05481320)
  v6  <- c(0.27660867,0.13512482,0.06304644,0.03865726,0.04005947,0.02674554,0.02529060,0.03131686,0.04890635)
  v7  <- c(0.22873691,0.11502189,0.05498601,0.03404040,0.03323795,0.02392815,0.02279867,0.02846645,0.04425618)
  v8  <- c(0.19451997,0.10017776,0.04862145,0.02930406,0.02651962,0.02232816,0.02098391,0.02623811,0.04053028)
  v9  <- c(0.16894062,0.08862642,0.04318970,0.02535628,0.02176994,0.02116887,0.01943677,0.02431963,0.03743835)
  v10 <- c(0.14877694,0.07908391,0.03873882,0.02263928,0.01873837,0.02088987,0.01859239,0.02294265,0.03493172)
  v11 <- c(0.13309193,0.07189149,0.03588324,0.02182041,0.02042703,0.01759560,0.01695284,0.02153022,0.03270945)
  v12 <- c(0.12002140,0.06596064,0.03300613,0.01934954,0.01611542,0.01791754,0.01605716,0.02014741,0.03073844)
  v13 <- c(0.10948010,0.06070346,0.03051874,0.01804678,0.01596139,0.01564856,0.01460856,0.01893691,0.02902857)
  v14 <- c(0.10024957,0.05580046,0.02823760,0.01710759,0.01542400,0.01556254,0.01443473,0.01845353,0.02767434)
  v15 <- c(0.09351666,0.05340978,0.02760381,0.01644813,0.01375095,0.01587233,0.01373580,0.01734491,0.02620031)
  v16 <- c(0.08612932,0.04925710,0.02558762,0.01560139,0.01368409,0.01412096,0.01302018,0.01662699,0.02505744)
  v17 <- c(0.08093024,0.04736369,0.02531342,0.01555457,0.01279792,0.01543469,0.01264462,0.01580596,0.02386073)
  v18 <- c(0.07556599,0.04404920,0.02322975,0.01425995,0.01205152,0.01425304,0.01224597,0.01542578,0.02300035)
  v19 <- c(0.07059393,0.04139277,0.02204407,0.01378145,0.01222891,0.01280162,0.01194975,0.01513357,0.02222269)
  
  #---------------   P A R A M E T E R I Z E D    M I X T U R E S  -------------------------------
  #-------------------------------- R = 4 components ---------------------------------------------
  # 20 <= yt <= 49 (range=3)
  nc3 <- 4
  coeff.w3 <- matrix(c(-5.644536495326424e-009,-1.266992312620729e-006,1.000000000000000e+000,4.730022618640170e+000,
                       7.299190941771616e-009,1.387196986613049e-006,1.000000000000000e+000,3.672627139063664e+000,-1.788056445701375e-008,    
                       -2.391966642312227e-006,1.000000000000000e+000,4.871566292199141e+000,9.259794020097245e-009,1.224613603301172e-006,
                       1.000000000000000e+000,3.215154075255811e+000),ncol=nc3)
  
  coeff.m3 <- matrix(c(-4.552797222245886e-005,4.638009105860874e-002,1.000000000000000e+000,9.627160143020319e-002,
                       1.143772956135763e+000,2.284729919321826e-005,-1.095058888699896e-002,1.000000000000000e+000,-1.690501643196285e-002,
                       1.944583776810219e+000,-3.900177124793736e-005,4.731686443506078e-002,1.000000000000000e+000,-5.610095109268514e-001,
                       -6.033854619021440e+000,-2.486737015927609e-005,2.978371498898014e-002,1.000000000000000e+000,-4.643825308039595e-002,
                       -1.105498133466881e+000),ncol=nc3)
  
  coeff.v3 <- matrix(c(-2.191015160635246e-005,9.939739739229481e-002,1.000000000000000e+000,9.208564449363621e-002,
                       7.148740127686324e-001,7.060864706964765e-005,1.143203813437594e-001,1.000000000000000e+000,1.548617268517774e-001,
                       2.428636911969079e+000,1.823003483480967e-004,1.675101633325085e-001,1.000000000000000e+000,2.735383307281154e-001,
                       4.861423133312201e+000,1.613752763707232e-004,1.943336591436946e-001,1.000000000000000e+000,2.797653349939537e-001,
                       3.840341872064528e+000),ncol=nc3)
  
  # calculation and transformation of (unstandardized) means and variances for range=3
  # Range 3: 20 <= yt <= 49
  nu  <- 20
  wr3 <- matrix(0, (49-20+1), nc3)
  mr3 <- matrix(0, (49-20+1), nc3)
  vr3 <- matrix(0, (49-20+1), nc3)
  
  for(i in 1:(49-20+1)){
    wr3[i,] <- (coeff.w3[3,]+coeff.w3[1,]*nu^2+coeff.w3[2,]*nu)/coeff.w3[4,]
    mr3[i,] <- ((coeff.m3[3,]+coeff.m3[1,]*nu^2+coeff.m3[2,]*nu)/(coeff.m3[4,]*nu+coeff.m3[5,]))*sqrt(trigamma(nu))-digamma(nu)
    vr3[i,] <- ((coeff.v3[3,]+coeff.v3[1,]*nu^2+coeff.v3[2,]*nu)/(coeff.v3[4,]*nu+coeff.v3[5,]))*trigamma(nu)
    nu <- nu + 1
  }
  
  #-------------------------------- R = 3 components ---------------------------------------------
  # 50 <= yt <= 439 (range=4)
  nc4 <- 3
  coeff.w4 <- matrix(c(-5.639545796991280e-010,4.698743002874532e-007,1.000000000000000e+000,4.730482920811330e+000,2.651836392450035e-010,
                       -1.280380156002802e-007,1.000000000000000e+000,2.093982718501769e+000,-2.384482520627535e-011,-1.227680572544847e-007,
                       1.000000000000000e+000,3.214956149674574e+000),ncol=nc4)
  coeff.m4 <- matrix(c(-1.653173201148335e-006,1.036578627632170e-002,1.000000000000000e+000,2.349390607953272e-002,1.432904956685477e+000,
                       -8.298537364426537e-007,5.017456263052972e-003,1.000000000000000e+000,5.123168011502721e-002,6.453910704667408e+000,-1.431525987300163e-006,
                       8.386323466104712e-003,1.000000000000000e+000,-1.841057020139425e-002,-1.410602407670769e+000),ncol=nc4)
  coeff.v4 <- matrix(c(-2.726183914412441e-007,2.788507874891710e-002,1.000000000000000e+000,2.777086294607445e-002,8.369406298984288e-001,
                       1.118379212729684e-006,2.433214514397419e-002,1.000000000000000e+000,2.778340896223197e-002,1.489387981224663e+000,2.197737873275589e-006,
                       3.186581505796005e-002,1.000000000000000e+000,3.808382220884354e-002,1.958805931276004e+000),ncol=nc4)
  
  # calculation and transformation of (unstandardized) means and variances for range=4
  # Range 4: 50 <= yt <= 439
  nu  <- 50
  wr4 <- matrix(0, (439-50+1), nc4)
  mr4 <- matrix(0, (439-50+1), nc4)
  vr4 <- matrix(0, (439-50+1), nc4)
  
  for(i in 1:(439-50+1)){
    wr4[i,] <- (coeff.w4[3,]+coeff.w4[1,]*nu^2+coeff.w4[2,]*nu)/coeff.w4[4,]
    mr4[i,] <- ((coeff.m4[3,]+coeff.m4[1,]*nu^2+coeff.m4[2,]*nu)/(coeff.m4[4,]*nu+coeff.m4[5,]))*sqrt(trigamma(nu))-digamma(nu)
    vr4[i,] <- ((coeff.v4[3,]+coeff.v4[1,]*nu^2+coeff.v4[2,]*nu)/(coeff.v4[4,]*nu+coeff.v4[5,]))*trigamma(nu)
    nu <- nu + 1
  }
  
  #-------------------------------- R = 2 components ---------------------------------------------
  # 440 <= yt <= 1599 (range=5)
  nc5 <- 2
  coeff.w5 <- matrix(c(1.034981541036597e-010,-2.445177000398938e-007,1.000000000000000e+000,1.451229377475864e+000,-2.291586556531707e-010,
                       5.414543692806514e-007,1.000000000000000e+000,3.216167113242079e+000),ncol=nc5)
  coeff.m5 <- matrix(c(-6.578325435644067e-008,1.648723149067166e-003,1.000000000000000e+000,1.594968525045459e-002,5.566082591106806e+000,
                       -6.292364160498604e-008,1.618047470065775e-003,1.000000000000000e+000,-7.091699113800587e-003,-2.516741952410371e+000),ncol=nc5)
  coeff.v5 <- matrix(c(-2.802162650788337e-009,4.051503597935380e-003,1.000000000000000e+000,4.018981069179972e-003,9.654061278849895e-001,
                       3.776558110733883e-008,5.022619018941299e-003,1.000000000000000e+000,5.525253413878772e-003,1.450507513327352e+000),ncol=nc5)
  
  # calculation and transformation of (unstandardized) means and variances for range=5
  # Range 4: 440 <= yt <= 1599
  nu  <- 440
  wr5 <- matrix(0, (1599-440+1), nc5)
  mr5 <- matrix(0, (1599-440+1), nc5)
  vr5 <- matrix(0, (1599-440+1), nc5)
  
  for(i in 1:(1599-440+1)){
    wr5[i,] <- (coeff.w5[3,]+coeff.w5[1,]*nu^2+coeff.w5[2,]*nu)/coeff.w5[4,]
    mr5[i,] <- ((coeff.m5[3,]+coeff.m5[1,]*nu^2+coeff.m5[2,]*nu)/(coeff.m5[4,]*nu+coeff.m5[5,]))*sqrt(trigamma(nu))-digamma(nu)
    vr5[i,] <- ((coeff.v5[3,]+coeff.v5[1,]*nu^2+coeff.v5[2,]*nu)/(coeff.v5[4,]*nu+coeff.v5[5,]))*trigamma(nu)
    nu <- nu + 1
  }
  
  #-------------------------------- R = 2 components ---------------------------------------------
  # 1600 <= yt <= 9999 (range=6)
  nc6 <- 2
  coeff.w6 <- matrix(c(-1.586037487404490e-013,3.575996226727867e-009,1.000000000000000e+000,2.228310599179340e+000,1.291237745205579e-013,
                       -2.911316152726367e-009,1.000000000000000e+000,1.814126328168031e+000),ncol=nc6)
  coeff.m6 <- matrix(c(-2.419956255676409e-009,3.245753451748892e-004,1.000000000000000e+000,1.895335618211674e-003,3.388553853864067e+000,
                       -2.419411092563945e-009,3.245669014250788e-004,1.000000000000000e+000,-2.327930564510444e-003,-4.162274939236667e+000),ncol=nc6)
  coeff.v6 <- matrix(c(-6.024563976875348e-011,-6.540694956580495e-004,1.000000000000000e+000,-6.582951415419203e-004,1.006399508694657e+000,
                       5.024317053887777e-010,8.898044793516080e-004,1.000000000000000e+000,9.246987493760628e-004,1.149073788967684e+000),ncol=nc6)
  
  # calculation and transformation of (unstandardized) means and variances for range=6
  # Range 4: 1600 <= yt <= 9999
  nu  <- 1600
  wr6 <- matrix(0, (9999-1600+1), nc6)
  mr6 <- matrix(0, (9999-1600+1), nc6)
  vr6 <- matrix(0, (9999-1600+1), nc6)
  
  for(i in 1:(9999-1600+1)){
    wr6[i,] <- (coeff.w6[3,]+coeff.w6[1,]*nu^2+coeff.w6[2,]*nu)/coeff.w6[4,]
    mr6[i,] <- ((coeff.m6[3,]+coeff.m6[1,]*nu^2+coeff.m6[2,]*nu)/(coeff.m6[4,]*nu+coeff.m6[5,]))*sqrt(trigamma(nu))-digamma(nu)
    vr6[i,] <- ((coeff.v6[3,]+coeff.v6[1,]*nu^2+coeff.v6[2,]*nu)/(coeff.v6[4,]*nu+coeff.v6[5,]))*trigamma(nu)
    nu <- nu + 1
  }
  
  #-------------------------------- R = 2 components ---------------------------------------------
  # 10000 <= yt <= 30000 (range=7)
  nc7 <- 2
  coeff.w7 <- matrix(c(-1.663426552872397e-014,1.141056828884990e-009,1.000000000000000e+000,2.228285989630589e+000,1.354267905471566e-014,
                       -9.289835742028532e-010,1.000000000000000e+000,1.814142639751731e+000),ncol=nc7)
  coeff.m7 <- matrix(c(-8.929405559006038e-011,6.319814700961324e-005,1.000000000000000e+000,4.785131212048377e-004,4.271922830906078e+000,
                       -8.931137480031157e-011,6.320244393309693e-005,1.000000000000000e+000,-5.877524860249395e-004,-5.247218808668549e+000),ncol=nc7)
  coeff.v7 <- matrix(c(-1.418731402291282e-012,-5.512224505288543e-006,1.000000000000000e+000,-5.638714069888806e-006,1.006201804172733e+000,
                       1.576782807097003e-011,1.914006058179041e-004,1.000000000000000e+000,1.959753272178233e-004,1.087101027065273e+000),ncol=nc7)
  
  # calculation and transformation of (unstandardized) means and variances for range=7
  # Range 4: 10000 <= yt <= 30000
  nu  <- 10000
  wr7 <- matrix(0, (30000-10000+1), nc7)
  mr7 <- matrix(0, (30000-10000+1), nc7)
  vr7 <- matrix(0, (30000-10000+1), nc7)
  
  for(i in 1:(30000-10000+1)){
    wr7[i,] <- (coeff.w7[3,]+coeff.w7[1,]*nu^2+coeff.w7[2,]*nu)/coeff.w7[4,]
    mr7[i,] <- ((coeff.m7[3,]+coeff.m7[1,]*nu^2+coeff.m7[2,]*nu)/(coeff.m7[4,]*nu+coeff.m7[5,]))*sqrt(trigamma(nu))-digamma(nu)
    vr7[i,] <- ((coeff.v7[3,]+coeff.v7[1,]*nu^2+coeff.v7[2,]*nu)/(coeff.v7[4,]*nu+coeff.v7[5,]))*trigamma(nu)
    nu <- nu + 1
  }
  
  #------------------------------------------------------------------------------------------------
  # weights
  w[1:4,]            <- rbind(w1,w2,w3,w4)
  w[5:19,1:9]        <- rbind(w5,w6,w7,w8,w9,w10,w11,w12,w13,w14,w15,w16,w17,w18,w19)
  w[20:49,1:4]       <- wr3
  w[50:439,1:3]      <- wr4
  w[440:1599,1:2]    <- wr5
  w[1600:9999,1:2]   <- wr6
  w[10000:30000,1:2] <- wr7
  
  # means
  m[1:4,]            <- rbind(m1,m2,m3,m4)
  m[5:19,1:9]        <- rbind(m5,m6,m7,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18,m19)
  m[20:49,1:4]       <- mr3
  m[50:439,1:3]      <- mr4
  m[440:1599,1:2]    <- mr5
  m[1600:9999,1:2]   <- mr6
  m[10000:30000,1:2] <- mr7
  
  # variances
  v[1:4,]            <- rbind(v1,v2,v3,v4)
  v[5:19,1:9]        <- rbind(v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15,v16,v17,v18,v19)
  v[20:49,1:4]       <- vr3
  v[50:439,1:3]      <- vr4
  v[440:1599,1:2]    <- vr5
  v[1600:9999,1:2]   <- vr6
  v[10000:30000,1:2] <- vr7
  
  
  return(list(m = m, v = v, w = w)
  )
}