R/ComputeMixture.R
In airbornemint/pogit: Bayesian Variable Selection for a Poisson-Logistic Model
## Approximation of the negative log-gamma distribution by Gaussian mixtures
## Compute mixture parameters, given degrees of freedom
## Author: R. Fruhwirth, HEPHY Vienna
## Version 1.0, 15/09/2006
## Adapted to R by Jesse Windle, UT Austin, 10/9/2012.
## This file was original found in R. Fruhwirth's public bitbucket repository.
## Version 2.0 of this file can be found in the bayesf MATLAB package found on
## S. Fruhwirth-Schnatter's website
## <http://statmath.wu.ac.at/~fruehwirth/monographie/>.
range <- matrix(ncol=2, byrow=TRUE, data=
c(1, 4, ## 10 components
5, 19, ## 9 components
20, 49, ## 4 components
50, 439, ## 3 components
440, 1599, ## 2 components
1600, 10000, ## 2 components
10000, 30000) ## 2 components
)
nrow.range = nrow(range);
..P = list()
..M = list()
..V = list()
Coeff.p = list()
Coeff.m = list()
Coeff.v = list()
## Explicit parameters for n up to 19
..P[[1]]=c(0.00396984425, 0.0396244597, 0.16776747, 0.147036501, 0.125306271, 0.1014852, 0.103758531, 0.115972617, 0.107065659, 0.0880134482);
..M[[1]]=c(3.55887454, 2.11415904, 0.968124631, 0.51537638, 0.145465449, -0.145346445, -0.416660312, -0.689002855, -0.974965634, -1.27310004);
..V[[1]]=c(2.62603032, 1.21263644, 0.66586521, 0.256650604, 0.120071142, 0.0649909219, 0.0473513798, 0.046639443, 0.0576541602, 0.0888536903);
..P[[2]]=c(0.00396860865, 0.039627049, 0.16777003, 0.147037759, 0.125304523, 0.101481714, 0.103759705, 0.115973128, 0.107065554, 0.0880119305);
..M[[2]]=c(2.78807754, 1.84979328, 0.94844169, 0.577673108, 0.223449219, -0.0831666379, -0.387174155, -0.69613969, -1.01843553, -1.34112844);
..V[[2]]=c(2.39619753, 1.16995764, 0.688870128, 0.307084756, 0.155644328, 0.0899360571, 0.0707828448, 0.0751755614, 0.0990773728, 0.15471843);
..P[[3]]=c(0.00396904734, 0.0396280642, 0.167770514, 0.14703607, 0.12530043, 0.10148242, 0.103759287, 0.115974323, 0.107066971, 0.0880128738);
..M[[3]]=c(2.43454312, 1.7315327, 0.942157786, 0.60208557, 0.251664821, -0.0644746918, -0.379817508, -0.696781518, -1.0293035, -1.35705784);
..V[[3]]=c(2.16215586, 1.11062998, 0.682294453, 0.324750601, 0.173204837, 0.108063698, 0.0917073596, 0.100257256, 0.131371692, 0.200024832);
..P[[4]]=c(0.00396883344, 0.039627504, 0.167771274, 0.147036659, 0.125301189, 0.101481755, 0.103760036, 0.115974339, 0.107065718, 0.0880126919);
..M[[4]]=c(2.30484474, 1.71231656, 0.952907078, 0.601128034, 0.252368847, -0.059032783, -0.375605704, -0.699071542, -1.03734211, -1.3609072);
..V[[4]]=c(1.92939547, 1.00671896, 0.638983371, 0.322852776, 0.18445103, 0.122217472, 0.106400052, 0.116936918, 0.154113316, 0.233525098);
..P[[5]]=c(0.0435820277, 0.167794347, 0.147040722, 0.125310654, 0.10147112, 0.10376347, 0.115973878, 0.107056197, 0.0880075845);
..M[[5]]=c(1.31113348, 0.963928895, 0.659198795, 0.240742429, -0.108844644, -0.252087404, -0.6546691, -1.04146524, -1.37874376);
..V[[5]]=c(1.5732832, 0.745075965, 0.340530976, 0.206325108, 0.206977107, 0.133034557, 0.123981078, 0.155417698, 0.247661591);
..P[[6]]=c(0.0435817033, 0.167795985, 0.1470426, 0.125311016, 0.101470666, 0.103763084, 0.115972864, 0.107055437, 0.0880066471);
..M[[6]]=c(1.25919247, 0.957217299, 0.66710982, 0.251658342, -0.125234491, -0.240137829, -0.64912733, -1.03921002, -1.37439461);
..V[[6]]=c(1.52550277, 0.745216293, 0.347702459, 0.213195645, 0.220928839, 0.147502243, 0.139478204, 0.17271313, 0.269719569);
..P[[7]]=c(0.0435798639, 0.167797087, 0.147042073, 0.125313291, 0.101470979, 0.103761847, 0.115973234, 0.107054351, 0.0880072753);
..M[[7]]=c(1.21602216, 0.94778507, 0.671484869, 0.265435387, -0.104709908, -0.24708343, -0.653441223, -1.04076324, -1.36988994);
..V[[7]]=c(1.48970429, 0.74910777, 0.35810967, 0.221696291, 0.216470192, 0.155837875, 0.148481868, 0.185394632, 0.28822907);
..P[[8]]=c(0.043578895, 0.167797426, 0.147041988, 0.125313875, 0.101470922, 0.103761581, 0.115973137, 0.107054001, 0.0880081751);
..M[[8]]=c(1.18027937, 0.939725546, 0.67760436, 0.293497817, -0.110879079, -0.257696481, -0.655613756, -1.0406543, -1.36465528);
..V[[8]]=c(1.46105103, 0.752441091, 0.365198621, 0.220104509, 0.199190433, 0.167708126, 0.15761138, 0.197076001, 0.304425302);
..P[[9]]=c(0.0435786725, 0.167797743, 0.1470428, 0.125313553, 0.101470946, 0.103761391, 0.115973188, 0.10705364, 0.0880080663);
..M[[9]]=c(1.14996911, 0.934206664, 0.686267712, 0.311595579, -0.112948479, -0.274222612, -0.653808807, -1.04092104, -1.35962481);
..V[[9]]=c(1.43764551, 0.754190306, 0.367534375, 0.215776065, 0.185257157, 0.180142183, 0.165402413, 0.206954388, 0.318591695);
..P[[10]]=c(0.0435779307, 0.167797788, 0.147042734, 0.125314068, 0.101471449, 0.10376142, 0.115973187, 0.107053473, 0.0880079505);
..M[[10]]=c(1.12841748, 0.932206841, 0.69102714, 0.319038554, -0.109581301, -0.302963892, -0.641448217, -1.03858769, -1.35274157);
..V[[10]]=c(1.41468216, 0.75198881, 0.368357589, 0.215271168, 0.178178434, 0.198636491, 0.176790288, 0.218155881, 0.332156859);
..P[[11]]=c(0.043576761, 0.167801375, 0.147042624, 0.125314075, 0.101470546, 0.103761069, 0.115973226, 0.107051966, 0.0880083593);
..M[[11]]=c(1.10451126, 0.925180162, 0.689947194, 0.309587296, -0.123979787, -0.239246368, -0.658582798, -1.03932069, -1.347407);
..V[[11]]=c(1.39851898, 0.755429842, 0.377058085, 0.229287048, 0.214645547, 0.18489307, 0.178139004, 0.226237823, 0.343708183);
..P[[12]]=c(0.0435771819, 0.167801103, 0.147042441, 0.125313864, 0.101470305, 0.103761519, 0.11597319, 0.107052417, 0.0880079809);
..M[[12]]=c(1.08624068, 0.918081034, 0.697616213, 0.330655882, -0.106424319, -0.290644969, -0.644517493, -1.04099153, -1.34370607);
..V[[12]]=c(1.38111403, 0.759024378, 0.379809227, 0.222659694, 0.185443843, 0.206181273, 0.184773494, 0.231840962, 0.353714302);
..P[[13]]=c(0.0435778469, 0.167800486, 0.147041951, 0.125313914, 0.101470076, 0.103761707, 0.115973611, 0.107052756, 0.0880076518);
..M[[13]]=c(1.0671125, 0.915784215, 0.70024231, 0.330800476, -0.125598534, -0.244656951, -0.661886313, -1.04447342, -1.33948264);
..V[[13]]=c(1.36922993, 0.759197249, 0.381687395, 0.225704876, 0.199623554, 0.195711194, 0.18270427, 0.236837387, 0.363050264);
..P[[14]]=c(0.0435786417, 0.16779926, 0.147042119, 0.125313391, 0.101470554, 0.103762378, 0.115973792, 0.107052537, 0.0880073289);
..M[[14]]=c(1.05721516, 0.918099637, 0.698999193, 0.325014717, -0.153165358, -0.225909041, -0.659788653, -1.03711782, -1.33064663);
..V[[14]]=c(1.35398708, 0.753650144, 0.381381699, 0.231057971, 0.208319112, 0.210190241, 0.194957855, 0.249236388, 0.373774124);
..P[[15]]=c(0.043581505, 0.167797871, 0.147043608, 0.125312521, 0.101469081, 0.103762173, 0.115973414, 0.107054363, 0.0880054639);
..M[[15]]=c(1.02150943, 0.896206397, 0.702224917, 0.344137939, -0.119895501, -0.256590721, -0.641185291, -1.03810889, -1.32943558);
..V[[15]]=c(1.35652837, 0.774748407, 0.400413698, 0.238592235, 0.199467639, 0.230239828, 0.19924794, 0.251600772, 0.380054821);
..P[[16]]=c(0.0435811435, 0.167798952, 0.147043687, 0.125312616, 0.101468918, 0.103762052, 0.115973417, 0.107053968, 0.0880052462);
..M[[16]]=c(1.02508782, 0.902555642, 0.699256309, 0.336391119, -0.121902141, -0.242730179, -0.6538063, -1.0385784, -1.32415888);
..V[[16]]=c(1.33546695, 0.763749521, 0.396745563, 0.241905327, 0.212176877, 0.218950701, 0.201882762, 0.257807637, 0.388524892);
..P[[17]]=c(0.0435812603, 0.167798873, 0.147044518, 0.125312321, 0.101468879, 0.103761729, 0.115972692, 0.107054049, 0.0880056789);
..M[[17]]=c(0.997274184, 0.88197491, 0.696155279, 0.3460138, -0.128981232, -0.227346713, -0.630322077, -1.03647508, -1.32316505);
..V[[17]]=c(1.33575722, 0.781739895, 0.417799104, 0.256728889, 0.211230256, 0.254750255, 0.208700024, 0.26087813, 0.393822372);
..P[[18]]=c(0.0435808733, 0.167799002, 0.147044529, 0.125312675, 0.101468951, 0.103761472, 0.115972643, 0.107053883, 0.0880059719);
..M[[18]]=c(0.995086849, 0.891409674, 0.70171109, 0.341992158, -0.127906113, -0.245952673, -0.638792902, -1.03392281, -1.31486719);
..V[[18]]=c(1.3227643, 0.771070524, 0.406631212, 0.249617029, 0.210958958, 0.249496089, 0.214362668, 0.270024593, 0.402615529);
..P[[19]]=c(0.0435807283, 0.167799231, 0.14704464, 0.12531292, 0.101468814, 0.103761275, 0.115972628, 0.107053662, 0.088006103);
..M[[19]]=c(0.997741814, 0.892112396, 0.698155553, 0.337731787, -0.122630195, -0.240878604, -0.651951415, -1.02898878, -1.3062535);
..V[[19]]=c(1.30630549, 0.765952536, 0.407914566, 0.255018833, 0.226289944, 0.236887588, 0.221124118, 0.280039124, 0.411219814);
## n from 20 to 49
Coeff.p[[3]] <- matrix(ncol=4, byrow=TRUE, data=
c(-5.644536495326424e-009, 7.299190941771616e-009, -1.788056445701375e-008, 9.259794020097245e-009,
-1.266992312620729e-006, 1.387196986613049e-006, -2.391966642312227e-006, 1.224613603301172e-006,
1.000000000000000e+000, 1.000000000000000e+000, 1.000000000000000e+000, 1.000000000000000e+000,
0, 0, 0, 0,
4.730022618640170e+000, 3.672627139063664e+000, 4.871566292199141e+000, 3.215154075255811e+000))
Coeff.m[[3]] <- matrix(ncol=4, byrow=TRUE, data=
c(-4.552797222245886e-005, 2.284729919321826e-005, -3.900177124793736e-005, -2.486737015927609e-005,
4.638009105860874e-002, -1.095058888699896e-002, 4.731686443506078e-002, 2.978371498898014e-002,
1.000000000000000e+000, 1.000000000000000e+000, 1.000000000000000e+000, 1.000000000000000e+000,
9.627160143020319e-002, -1.690501643196285e-002, -5.610095109268514e-001, -4.643825308039595e-002,
1.143772956135763e+000, 1.944583776810219e+000, -6.033854619021440e+000, -1.105498133466881e+000))
Coeff.v[[3]] <- matrix(ncol=4, byrow=TRUE, data=
c(-2.191015160635246e-005, 7.060864706964765e-005, 1.823003483480967e-004, 1.613752763707232e-004,
9.939739739229481e-002, 1.143203813437594e-001, 1.675101633325085e-001, 1.943336591436946e-001,
1.000000000000000e+000, 1.000000000000000e+000, 1.000000000000000e+000, 1.000000000000000e+000,
9.208564449363621e-002, 1.548617268517774e-001, 2.735383307281154e-001, 2.797653349939537e-001,
7.148740127686324e-001, 2.428636911969079e+000, 4.861423133312201e+000, 3.840341872064528e+000));
## n from 50 to 439
Coeff.p[[4]] <- matrix(ncol=3, byrow=TRUE, data=
c(-5.639545796991280e-010, 2.651836392450035e-010, -2.384482520627535e-011,
4.698743002874532e-007, -1.280380156002802e-007, -1.227680572544847e-007,
1.000000000000000e+000, 1.000000000000000e+000, 1.000000000000000e+000,
0, 0, 0,
4.730482920811330e+000, 2.093982718501769e+000, 3.214956149674574e+000))
Coeff.m[[4]] <- matrix(ncol=3, byrow=TRUE, data=
c(-1.653173201148335e-006, -8.298537364426537e-007, -1.431525987300163e-006,
1.036578627632170e-002, 5.017456263052972e-003, 8.386323466104712e-003,
1.000000000000000e+000, 1.000000000000000e+000, 1.000000000000000e+000,
2.349390607953272e-002, 5.123168011502721e-002, -1.841057020139425e-002,
1.432904956685477e+000, 6.453910704667408e+000, -1.410602407670769e+000))
Coeff.v[[4]] <- matrix(ncol=3, byrow=TRUE, data=
c(-2.726183914412441e-007, 1.118379212729684e-006, 2.197737873275589e-006,
2.788507874891710e-002, 2.433214514397419e-002, 3.186581505796005e-002,
1.000000000000000e+000, 1.000000000000000e+000, 1.000000000000000e+000,
2.777086294607445e-002, 2.778340896223197e-002, 3.808382220884354e-002,
8.369406298984288e-001, 1.489387981224663e+000, 1.958805931276004e+000))
## n from 440 to 1599
Coeff.p[[5]] <- matrix(ncol=2, byrow=TRUE, data=
c(1.034981541036597e-010, -2.291586556531707e-010,
-2.445177000398938e-007, 5.414543692806514e-007,
1.000000000000000e+000, 1.000000000000000e+000,
0, 0,
1.451229377475864e+000, 3.216167113242079e+000))
Coeff.m[[5]] <- matrix(ncol=2, byrow=TRUE, data=
c(-6.578325435644067e-008, -6.292364160498604e-008,
1.648723149067166e-003, 1.618047470065775e-003,
1.000000000000000e+000, 1.000000000000000e+000,
1.594968525045459e-002, -7.091699113800587e-003,
5.566082591106806e+000, -2.516741952410371e+000))
Coeff.v[[5]] <- matrix(ncol=2, byrow=TRUE, data=
c(-2.802162650788337e-009, 3.776558110733883e-008,
4.051503597935380e-003, 5.022619018941299e-003,
1.000000000000000e+000, 1.000000000000000e+000,
4.018981069179972e-003, 5.525253413878772e-003,
9.654061278849895e-001, 1.450507513327352e+000))
## n from 1600 to 10000
Coeff.p[[6]] <- matrix(ncol=2, byrow=TRUE, data=
c(-1.586037487404490e-013, 1.291237745205579e-013,
3.575996226727867e-009, -2.911316152726367e-009,
1.000000000000000e+000, 1.000000000000000e+000,
0, 0,
2.228310599179340e+000, 1.814126328168031e+000))
Coeff.m[[6]] <- matrix(ncol=2, byrow=TRUE, data=
c(-2.419956255676409e-009, -2.419411092563945e-009,
3.245753451748892e-004, 3.245669014250788e-004,
1.000000000000000e+000, 1.000000000000000e+000,
1.895335618211674e-003, -2.327930564510444e-003,
3.388553853864067e+000, -4.162274939236667e+000))
Coeff.v[[6]] <- matrix(ncol=2, byrow=TRUE, data=
c(-6.024563976875348e-011, 5.024317053887777e-010,
-6.540694956580495e-004, 8.898044793516080e-004,
1.000000000000000e+000, 1.000000000000000e+000,
-6.582951415419203e-004, 9.246987493760628e-004,
1.006399508694657e+000, 1.149073788967684e+000))
## n from 10000 to 30000
Coeff.p[[7]] <- matrix(ncol=2, byrow=TRUE, data=
c(-1.663426552872397e-014, 1.354267905471566e-014,
1.141056828884990e-009, -9.289835742028532e-010,
1.000000000000000e+000, 1.000000000000000e+000,
0, 0,
2.228285989630589e+000, 1.814142639751731e+000))
Coeff.m[[7]] <- matrix(ncol=2, byrow=TRUE, data=
c(-8.929405559006038e-011, -8.931137480031157e-011,
6.319814700961324e-005, 6.320244393309693e-005,
1.000000000000000e+000, 1.000000000000000e+000,
4.785131212048377e-004, -5.877524860249395e-004,
4.271922830906078e+000, -5.247218808668549e+000))
Coeff.v[[7]] <- matrix(ncol=2, byrow=TRUE, data=
c(-1.418731402291282e-012, 1.576782807097003e-011,
-5.512224505288543e-006, 1.914006058179041e-004,
1.000000000000000e+000, 1.000000000000000e+000,
-5.638714069888806e-006, 1.959753272178233e-004,
1.006201804172733e+000, 1.087101027065273e+000))
################################################################################
polyval <- function(coeff, x, n=length(coeff))
{
## columns are monomials evalued at x.
X = t(outer(x, (n-1):0, FUN="^"))
y = colSums(X * coeff);
## y = sum( coeff * x^((n-1):0) );
y
}
ratval <- function(p, x)
{
## Evaluate rational function at x from coefficients p
n = length(p)
m=(1:n)[p==1];
a=p[1:m];
b=p[(m+1):n];
y=polyval(a, x)/polyval(b, x);
y
}
compute.mixture.lg <- function(n)
{
## n: degrees of freedom.
p = rep(0, 10)
m = rep(0, 10)
v = rep(0, 10)
## Check input
if (n<=0) {
print('n must be positive');
return(NA);
}
## if n<=range(3, 1) & n~=round(n)
if (n<=range[3,1] && n!=round(n)) {
print(paste('If n <= ', 19, ' it must be integer'));
return(NA);
}
## Find appropriate range
if (n>range[length(range)]) {
p=1;
m=0;
v=1;
} else {
r = max( (1:nrow.range)[rowSums(range<=n)>0] );
if (r<=2) { ## Individual mixtures
p=..P[[n]];
m=..M[[n]];
v=..V[[n]];
} else { ## Parameterized mixtures
coeff_p=Coeff.p[[r]];
coeff_m=Coeff.m[[r]];
coeff_v=Coeff.v[[r]];
nc=ncol(coeff_p);
p = rep(0, nc)
m = rep(0, nc)
v = rep(0, nc)
for (ic in 1:nc) {
p[ic]=ratval(coeff_p[, ic], n);
m[ic]=ratval(coeff_m[, ic], n);
v[ic]=ratval(coeff_v[, ic], n);
}
}
}
## Scale with mean and standard deviation:
## Mean and variance of - log Ga(n,1). The normal mixture is for the
## normalized version of - log Ga(n,1) so we have to transform the means and
## variances. You can actually compute these moments by hand via the MGF
## (transform to y = ln(x) simplify and then back). Evidently, psi is called
## the polygamma function. In R it is called psigamma.
mu = -1 * digamma(n); ## psigamma(n, 0)
sig2 = trigamma(n); ## psigamma(n, 1)
sig = sqrt(sig2);
m = sig*m + mu;
v = sig2 * v;
nc=length(p);
out = list("p"=p, "m"=m, "v"=v, "nc"=nc);
out
}
sample.mixture <- function(nsamp, nmix)
{
r = sample.int(nmix$nc, size=nsamp, replace=TRUE, prob=nmix$p);
y = rnorm(nsamp, nmix$m[r], sqrt(nmix$v[r]));
y
}
################################################################################
test.nmix <- function(df, nsamp)
{
n = length(df)
report = matrix(nrow=n, ncol=4);
colnames(report) = c("R.mean", "R.sd", "nmix.mean", "nmix.sd");
for (i in 1:n) {
samp.R = -1 * log(rgamma(nsamp, df[i]))
nmix = compute.mixture(df[i]);
samp.N = sample.mixture(nsamp, nmix);
report[i,1:2] = c(mean(samp.R), sd(samp.R));
report[i,3:4] = c(mean(samp.N), sd(samp.N));
}
report
}
if (FALSE) {
the.test = test.nmix(1:5000, 1000)
diff.mean = the.test[,1] - the.test[,3]
diff.sd = the.test[,2] - the.test[,4]
sd(diff.mean)
sd(diff.sd)
}
################################################################################
## COMPUTE MIXTURE GENERAL ##
################################################################################
compute.mixture <- function(shape, type=c("log.gamma", "logistic.iii"))
{
nm = NA
if (type[1]=="log.gamma") {
nm = compute.mixture.lg(shape)
} else
if (type[1]=="logistic.iii") {
nm = compmix(shape)
nm$m = rep(0, length(nm$v))
} else {
cat("compute.mixture: unrecognized type:", type[1], "\n");
}
nm
}
airbornemint/pogit documentation built on May 31, 2019, 1:49 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
## Approximation of the negative log-gamma distribution by Gaussian mixtures
## Compute mixture parameters, given degrees of freedom
## Author: R. Fruhwirth, HEPHY Vienna
## Version 1.0, 15/09/2006
## Adapted to R by Jesse Windle, UT Austin, 10/9/2012.
## This file was original found in R. Fruhwirth's public bitbucket repository.
## Version 2.0 of this file can be found in the bayesf MATLAB package found on
## S. Fruhwirth-Schnatter's website
## <http://statmath.wu.ac.at/~fruehwirth/monographie/>.
range <- matrix(ncol=2, byrow=TRUE, data=
c(1, 4, ## 10 components
5, 19, ## 9 components
20, 49, ## 4 components
50, 439, ## 3 components
440, 1599, ## 2 components
1600, 10000, ## 2 components
10000, 30000) ## 2 components
)
nrow.range = nrow(range);
..P = list()
..M = list()
..V = list()
Coeff.p = list()
Coeff.m = list()
Coeff.v = list()
## Explicit parameters for n up to 19
..P[[1]]=c(0.00396984425, 0.0396244597, 0.16776747, 0.147036501, 0.125306271, 0.1014852, 0.103758531, 0.115972617, 0.107065659, 0.0880134482);
..M[[1]]=c(3.55887454, 2.11415904, 0.968124631, 0.51537638, 0.145465449, -0.145346445, -0.416660312, -0.689002855, -0.974965634, -1.27310004);
..V[[1]]=c(2.62603032, 1.21263644, 0.66586521, 0.256650604, 0.120071142, 0.0649909219, 0.0473513798, 0.046639443, 0.0576541602, 0.0888536903);
..P[[2]]=c(0.00396860865, 0.039627049, 0.16777003, 0.147037759, 0.125304523, 0.101481714, 0.103759705, 0.115973128, 0.107065554, 0.0880119305);
..M[[2]]=c(2.78807754, 1.84979328, 0.94844169, 0.577673108, 0.223449219, -0.0831666379, -0.387174155, -0.69613969, -1.01843553, -1.34112844);
..V[[2]]=c(2.39619753, 1.16995764, 0.688870128, 0.307084756, 0.155644328, 0.0899360571, 0.0707828448, 0.0751755614, 0.0990773728, 0.15471843);
..P[[3]]=c(0.00396904734, 0.0396280642, 0.167770514, 0.14703607, 0.12530043, 0.10148242, 0.103759287, 0.115974323, 0.107066971, 0.0880128738);
..M[[3]]=c(2.43454312, 1.7315327, 0.942157786, 0.60208557, 0.251664821, -0.0644746918, -0.379817508, -0.696781518, -1.0293035, -1.35705784);
..V[[3]]=c(2.16215586, 1.11062998, 0.682294453, 0.324750601, 0.173204837, 0.108063698, 0.0917073596, 0.100257256, 0.131371692, 0.200024832);
..P[[4]]=c(0.00396883344, 0.039627504, 0.167771274, 0.147036659, 0.125301189, 0.101481755, 0.103760036, 0.115974339, 0.107065718, 0.0880126919);
..M[[4]]=c(2.30484474, 1.71231656, 0.952907078, 0.601128034, 0.252368847, -0.059032783, -0.375605704, -0.699071542, -1.03734211, -1.3609072);
..V[[4]]=c(1.92939547, 1.00671896, 0.638983371, 0.322852776, 0.18445103, 0.122217472, 0.106400052, 0.116936918, 0.154113316, 0.233525098);
..P[[5]]=c(0.0435820277, 0.167794347, 0.147040722, 0.125310654, 0.10147112, 0.10376347, 0.115973878, 0.107056197, 0.0880075845);
..M[[5]]=c(1.31113348, 0.963928895, 0.659198795, 0.240742429, -0.108844644, -0.252087404, -0.6546691, -1.04146524, -1.37874376);
..V[[5]]=c(1.5732832, 0.745075965, 0.340530976, 0.206325108, 0.206977107, 0.133034557, 0.123981078, 0.155417698, 0.247661591);
..P[[6]]=c(0.0435817033, 0.167795985, 0.1470426, 0.125311016, 0.101470666, 0.103763084, 0.115972864, 0.107055437, 0.0880066471);
..M[[6]]=c(1.25919247, 0.957217299, 0.66710982, 0.251658342, -0.125234491, -0.240137829, -0.64912733, -1.03921002, -1.37439461);
..V[[6]]=c(1.52550277, 0.745216293, 0.347702459, 0.213195645, 0.220928839, 0.147502243, 0.139478204, 0.17271313, 0.269719569);
..P[[7]]=c(0.0435798639, 0.167797087, 0.147042073, 0.125313291, 0.101470979, 0.103761847, 0.115973234, 0.107054351, 0.0880072753);
..M[[7]]=c(1.21602216, 0.94778507, 0.671484869, 0.265435387, -0.104709908, -0.24708343, -0.653441223, -1.04076324, -1.36988994);
..V[[7]]=c(1.48970429, 0.74910777, 0.35810967, 0.221696291, 0.216470192, 0.155837875, 0.148481868, 0.185394632, 0.28822907);
..P[[8]]=c(0.043578895, 0.167797426, 0.147041988, 0.125313875, 0.101470922, 0.103761581, 0.115973137, 0.107054001, 0.0880081751);
..M[[8]]=c(1.18027937, 0.939725546, 0.67760436, 0.293497817, -0.110879079, -0.257696481, -0.655613756, -1.0406543, -1.36465528);
..V[[8]]=c(1.46105103, 0.752441091, 0.365198621, 0.220104509, 0.199190433, 0.167708126, 0.15761138, 0.197076001, 0.304425302);
..P[[9]]=c(0.0435786725, 0.167797743, 0.1470428, 0.125313553, 0.101470946, 0.103761391, 0.115973188, 0.10705364, 0.0880080663);
..M[[9]]=c(1.14996911, 0.934206664, 0.686267712, 0.311595579, -0.112948479, -0.274222612, -0.653808807, -1.04092104, -1.35962481);
..V[[9]]=c(1.43764551, 0.754190306, 0.367534375, 0.215776065, 0.185257157, 0.180142183, 0.165402413, 0.206954388, 0.318591695);
..P[[10]]=c(0.0435779307, 0.167797788, 0.147042734, 0.125314068, 0.101471449, 0.10376142, 0.115973187, 0.107053473, 0.0880079505);
..M[[10]]=c(1.12841748, 0.932206841, 0.69102714, 0.319038554, -0.109581301, -0.302963892, -0.641448217, -1.03858769, -1.35274157);
..V[[10]]=c(1.41468216, 0.75198881, 0.368357589, 0.215271168, 0.178178434, 0.198636491, 0.176790288, 0.218155881, 0.332156859);
..P[[11]]=c(0.043576761, 0.167801375, 0.147042624, 0.125314075, 0.101470546, 0.103761069, 0.115973226, 0.107051966, 0.0880083593);
..M[[11]]=c(1.10451126, 0.925180162, 0.689947194, 0.309587296, -0.123979787, -0.239246368, -0.658582798, -1.03932069, -1.347407);
..V[[11]]=c(1.39851898, 0.755429842, 0.377058085, 0.229287048, 0.214645547, 0.18489307, 0.178139004, 0.226237823, 0.343708183);
..P[[12]]=c(0.0435771819, 0.167801103, 0.147042441, 0.125313864, 0.101470305, 0.103761519, 0.11597319, 0.107052417, 0.0880079809);
..M[[12]]=c(1.08624068, 0.918081034, 0.697616213, 0.330655882, -0.106424319, -0.290644969, -0.644517493, -1.04099153, -1.34370607);
..V[[12]]=c(1.38111403, 0.759024378, 0.379809227, 0.222659694, 0.185443843, 0.206181273, 0.184773494, 0.231840962, 0.353714302);
..P[[13]]=c(0.0435778469, 0.167800486, 0.147041951, 0.125313914, 0.101470076, 0.103761707, 0.115973611, 0.107052756, 0.0880076518);
..M[[13]]=c(1.0671125, 0.915784215, 0.70024231, 0.330800476, -0.125598534, -0.244656951, -0.661886313, -1.04447342, -1.33948264);
..V[[13]]=c(1.36922993, 0.759197249, 0.381687395, 0.225704876, 0.199623554, 0.195711194, 0.18270427, 0.236837387, 0.363050264);
..P[[14]]=c(0.0435786417, 0.16779926, 0.147042119, 0.125313391, 0.101470554, 0.103762378, 0.115973792, 0.107052537, 0.0880073289);
..M[[14]]=c(1.05721516, 0.918099637, 0.698999193, 0.325014717, -0.153165358, -0.225909041, -0.659788653, -1.03711782, -1.33064663);
..V[[14]]=c(1.35398708, 0.753650144, 0.381381699, 0.231057971, 0.208319112, 0.210190241, 0.194957855, 0.249236388, 0.373774124);
..P[[15]]=c(0.043581505, 0.167797871, 0.147043608, 0.125312521, 0.101469081, 0.103762173, 0.115973414, 0.107054363, 0.0880054639);
..M[[15]]=c(1.02150943, 0.896206397, 0.702224917, 0.344137939, -0.119895501, -0.256590721, -0.641185291, -1.03810889, -1.32943558);
..V[[15]]=c(1.35652837, 0.774748407, 0.400413698, 0.238592235, 0.199467639, 0.230239828, 0.19924794, 0.251600772, 0.380054821);
..P[[16]]=c(0.0435811435, 0.167798952, 0.147043687, 0.125312616, 0.101468918, 0.103762052, 0.115973417, 0.107053968, 0.0880052462);
..M[[16]]=c(1.02508782, 0.902555642, 0.699256309, 0.336391119, -0.121902141, -0.242730179, -0.6538063, -1.0385784, -1.32415888);
..V[[16]]=c(1.33546695, 0.763749521, 0.396745563, 0.241905327, 0.212176877, 0.218950701, 0.201882762, 0.257807637, 0.388524892);
..P[[17]]=c(0.0435812603, 0.167798873, 0.147044518, 0.125312321, 0.101468879, 0.103761729, 0.115972692, 0.107054049, 0.0880056789);
..M[[17]]=c(0.997274184, 0.88197491, 0.696155279, 0.3460138, -0.128981232, -0.227346713, -0.630322077, -1.03647508, -1.32316505);
..V[[17]]=c(1.33575722, 0.781739895, 0.417799104, 0.256728889, 0.211230256, 0.254750255, 0.208700024, 0.26087813, 0.393822372);
..P[[18]]=c(0.0435808733, 0.167799002, 0.147044529, 0.125312675, 0.101468951, 0.103761472, 0.115972643, 0.107053883, 0.0880059719);
..M[[18]]=c(0.995086849, 0.891409674, 0.70171109, 0.341992158, -0.127906113, -0.245952673, -0.638792902, -1.03392281, -1.31486719);
..V[[18]]=c(1.3227643, 0.771070524, 0.406631212, 0.249617029, 0.210958958, 0.249496089, 0.214362668, 0.270024593, 0.402615529);
..P[[19]]=c(0.0435807283, 0.167799231, 0.14704464, 0.12531292, 0.101468814, 0.103761275, 0.115972628, 0.107053662, 0.088006103);
..M[[19]]=c(0.997741814, 0.892112396, 0.698155553, 0.337731787, -0.122630195, -0.240878604, -0.651951415, -1.02898878, -1.3062535);
..V[[19]]=c(1.30630549, 0.765952536, 0.407914566, 0.255018833, 0.226289944, 0.236887588, 0.221124118, 0.280039124, 0.411219814);
## n from 20 to 49
Coeff.p[[3]] <- matrix(ncol=4, byrow=TRUE, data=
c(-5.644536495326424e-009, 7.299190941771616e-009, -1.788056445701375e-008, 9.259794020097245e-009,
-1.266992312620729e-006, 1.387196986613049e-006, -2.391966642312227e-006, 1.224613603301172e-006,
1.000000000000000e+000, 1.000000000000000e+000, 1.000000000000000e+000, 1.000000000000000e+000,
0, 0, 0, 0,
4.730022618640170e+000, 3.672627139063664e+000, 4.871566292199141e+000, 3.215154075255811e+000))
Coeff.m[[3]] <- matrix(ncol=4, byrow=TRUE, data=
c(-4.552797222245886e-005, 2.284729919321826e-005, -3.900177124793736e-005, -2.486737015927609e-005,
4.638009105860874e-002, -1.095058888699896e-002, 4.731686443506078e-002, 2.978371498898014e-002,
1.000000000000000e+000, 1.000000000000000e+000, 1.000000000000000e+000, 1.000000000000000e+000,
9.627160143020319e-002, -1.690501643196285e-002, -5.610095109268514e-001, -4.643825308039595e-002,
1.143772956135763e+000, 1.944583776810219e+000, -6.033854619021440e+000, -1.105498133466881e+000))
Coeff.v[[3]] <- matrix(ncol=4, byrow=TRUE, data=
c(-2.191015160635246e-005, 7.060864706964765e-005, 1.823003483480967e-004, 1.613752763707232e-004,
9.939739739229481e-002, 1.143203813437594e-001, 1.675101633325085e-001, 1.943336591436946e-001,
1.000000000000000e+000, 1.000000000000000e+000, 1.000000000000000e+000, 1.000000000000000e+000,
9.208564449363621e-002, 1.548617268517774e-001, 2.735383307281154e-001, 2.797653349939537e-001,
7.148740127686324e-001, 2.428636911969079e+000, 4.861423133312201e+000, 3.840341872064528e+000));
## n from 50 to 439
Coeff.p[[4]] <- matrix(ncol=3, byrow=TRUE, data=
c(-5.639545796991280e-010, 2.651836392450035e-010, -2.384482520627535e-011,
4.698743002874532e-007, -1.280380156002802e-007, -1.227680572544847e-007,
1.000000000000000e+000, 1.000000000000000e+000, 1.000000000000000e+000,
0, 0, 0,
4.730482920811330e+000, 2.093982718501769e+000, 3.214956149674574e+000))
Coeff.m[[4]] <- matrix(ncol=3, byrow=TRUE, data=
c(-1.653173201148335e-006, -8.298537364426537e-007, -1.431525987300163e-006,
1.036578627632170e-002, 5.017456263052972e-003, 8.386323466104712e-003,
1.000000000000000e+000, 1.000000000000000e+000, 1.000000000000000e+000,
2.349390607953272e-002, 5.123168011502721e-002, -1.841057020139425e-002,
1.432904956685477e+000, 6.453910704667408e+000, -1.410602407670769e+000))
Coeff.v[[4]] <- matrix(ncol=3, byrow=TRUE, data=
c(-2.726183914412441e-007, 1.118379212729684e-006, 2.197737873275589e-006,
2.788507874891710e-002, 2.433214514397419e-002, 3.186581505796005e-002,
1.000000000000000e+000, 1.000000000000000e+000, 1.000000000000000e+000,
2.777086294607445e-002, 2.778340896223197e-002, 3.808382220884354e-002,
8.369406298984288e-001, 1.489387981224663e+000, 1.958805931276004e+000))
## n from 440 to 1599
Coeff.p[[5]] <- matrix(ncol=2, byrow=TRUE, data=
c(1.034981541036597e-010, -2.291586556531707e-010,
-2.445177000398938e-007, 5.414543692806514e-007,
1.000000000000000e+000, 1.000000000000000e+000,
0, 0,
1.451229377475864e+000, 3.216167113242079e+000))
Coeff.m[[5]] <- matrix(ncol=2, byrow=TRUE, data=
c(-6.578325435644067e-008, -6.292364160498604e-008,
1.648723149067166e-003, 1.618047470065775e-003,
1.000000000000000e+000, 1.000000000000000e+000,
1.594968525045459e-002, -7.091699113800587e-003,
5.566082591106806e+000, -2.516741952410371e+000))
Coeff.v[[5]] <- matrix(ncol=2, byrow=TRUE, data=
c(-2.802162650788337e-009, 3.776558110733883e-008,
4.051503597935380e-003, 5.022619018941299e-003,
1.000000000000000e+000, 1.000000000000000e+000,
4.018981069179972e-003, 5.525253413878772e-003,
9.654061278849895e-001, 1.450507513327352e+000))
## n from 1600 to 10000
Coeff.p[[6]] <- matrix(ncol=2, byrow=TRUE, data=
c(-1.586037487404490e-013, 1.291237745205579e-013,
3.575996226727867e-009, -2.911316152726367e-009,
1.000000000000000e+000, 1.000000000000000e+000,
0, 0,
2.228310599179340e+000, 1.814126328168031e+000))
Coeff.m[[6]] <- matrix(ncol=2, byrow=TRUE, data=
c(-2.419956255676409e-009, -2.419411092563945e-009,
3.245753451748892e-004, 3.245669014250788e-004,
1.000000000000000e+000, 1.000000000000000e+000,
1.895335618211674e-003, -2.327930564510444e-003,
3.388553853864067e+000, -4.162274939236667e+000))
Coeff.v[[6]] <- matrix(ncol=2, byrow=TRUE, data=
c(-6.024563976875348e-011, 5.024317053887777e-010,
-6.540694956580495e-004, 8.898044793516080e-004,
1.000000000000000e+000, 1.000000000000000e+000,
-6.582951415419203e-004, 9.246987493760628e-004,
1.006399508694657e+000, 1.149073788967684e+000))
## n from 10000 to 30000
Coeff.p[[7]] <- matrix(ncol=2, byrow=TRUE, data=
c(-1.663426552872397e-014, 1.354267905471566e-014,
1.141056828884990e-009, -9.289835742028532e-010,
1.000000000000000e+000, 1.000000000000000e+000,
0, 0,
2.228285989630589e+000, 1.814142639751731e+000))
Coeff.m[[7]] <- matrix(ncol=2, byrow=TRUE, data=
c(-8.929405559006038e-011, -8.931137480031157e-011,
6.319814700961324e-005, 6.320244393309693e-005,
1.000000000000000e+000, 1.000000000000000e+000,
4.785131212048377e-004, -5.877524860249395e-004,
4.271922830906078e+000, -5.247218808668549e+000))
Coeff.v[[7]] <- matrix(ncol=2, byrow=TRUE, data=
c(-1.418731402291282e-012, 1.576782807097003e-011,
-5.512224505288543e-006, 1.914006058179041e-004,
1.000000000000000e+000, 1.000000000000000e+000,
-5.638714069888806e-006, 1.959753272178233e-004,
1.006201804172733e+000, 1.087101027065273e+000))
################################################################################
polyval <- function(coeff, x, n=length(coeff))
{
## columns are monomials evalued at x.
X = t(outer(x, (n-1):0, FUN="^"))
y = colSums(X * coeff);
## y = sum( coeff * x^((n-1):0) );
y
}
ratval <- function(p, x)
{
## Evaluate rational function at x from coefficients p
n = length(p)
m=(1:n)[p==1];
a=p[1:m];
b=p[(m+1):n];
y=polyval(a, x)/polyval(b, x);
y
}
compute.mixture.lg <- function(n)
{
## n: degrees of freedom.
p = rep(0, 10)
m = rep(0, 10)
v = rep(0, 10)
## Check input
if (n<=0) {
print('n must be positive');
return(NA);
}
## if n<=range(3, 1) & n~=round(n)
if (n<=range[3,1] && n!=round(n)) {
print(paste('If n <= ', 19, ' it must be integer'));
return(NA);
}
## Find appropriate range
if (n>range[length(range)]) {
p=1;
m=0;
v=1;
} else {
r = max( (1:nrow.range)[rowSums(range<=n)>0] );
if (r<=2) { ## Individual mixtures
p=..P[[n]];
m=..M[[n]];
v=..V[[n]];
} else { ## Parameterized mixtures
coeff_p=Coeff.p[[r]];
coeff_m=Coeff.m[[r]];
coeff_v=Coeff.v[[r]];
nc=ncol(coeff_p);
p = rep(0, nc)
m = rep(0, nc)
v = rep(0, nc)
for (ic in 1:nc) {
p[ic]=ratval(coeff_p[, ic], n);
m[ic]=ratval(coeff_m[, ic], n);
v[ic]=ratval(coeff_v[, ic], n);
}
}
}
## Scale with mean and standard deviation:
## Mean and variance of - log Ga(n,1). The normal mixture is for the
## normalized version of - log Ga(n,1) so we have to transform the means and
## variances. You can actually compute these moments by hand via the MGF
## (transform to y = ln(x) simplify and then back). Evidently, psi is called
## the polygamma function. In R it is called psigamma.
mu = -1 * digamma(n); ## psigamma(n, 0)
sig2 = trigamma(n); ## psigamma(n, 1)
sig = sqrt(sig2);
m = sig*m + mu;
v = sig2 * v;
nc=length(p);
out = list("p"=p, "m"=m, "v"=v, "nc"=nc);
out
}
sample.mixture <- function(nsamp, nmix)
{
r = sample.int(nmix$nc, size=nsamp, replace=TRUE, prob=nmix$p);
y = rnorm(nsamp, nmix$m[r], sqrt(nmix$v[r]));
y
}
################################################################################
test.nmix <- function(df, nsamp)
{
n = length(df)
report = matrix(nrow=n, ncol=4);
colnames(report) = c("R.mean", "R.sd", "nmix.mean", "nmix.sd");
for (i in 1:n) {
samp.R = -1 * log(rgamma(nsamp, df[i]))
nmix = compute.mixture(df[i]);
samp.N = sample.mixture(nsamp, nmix);
report[i,1:2] = c(mean(samp.R), sd(samp.R));
report[i,3:4] = c(mean(samp.N), sd(samp.N));
}
report
}
if (FALSE) {
the.test = test.nmix(1:5000, 1000)
diff.mean = the.test[,1] - the.test[,3]
diff.sd = the.test[,2] - the.test[,4]
sd(diff.mean)
sd(diff.sd)
}
################################################################################
## COMPUTE MIXTURE GENERAL ##
################################################################################
compute.mixture <- function(shape, type=c("log.gamma", "logistic.iii"))
{
nm = NA
if (type[1]=="log.gamma") {
nm = compute.mixture.lg(shape)
} else
if (type[1]=="logistic.iii") {
nm = compmix(shape)
nm$m = rep(0, length(nm$v))
} else {
cat("compute.mixture: unrecognized type:", type[1], "\n");
}
nm
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.