Nothing
## 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)
)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.