Skip to content

Commit

Permalink
issue #35 fixed; caused by pigamma not being initialized to 0 so that…
Browse files Browse the repository at this point in the history 
… sampling ended early
  • Loading branch information
@merliseclyde
merliseclyde committed Aug 14, 2018
1 parent bd1a17f commit 03d80b1
Show file tree
Hide file tree
Showing 3 changed files with 71 additions and 60 deletions.
2 changes: 1 addition & 1 deletion NEWS.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,11 +18,11 @@

## Bugs

* fixed [issue #35](https://github.com/merliseclyde/BAS/issues/35) for `method="MCMC+BAS"` in `bas.glm` in `glm_mcmcbas.c` when no values are provided for `MCMC.iterations` or `n.models` and defaults are used. Added unit test in `test-bas-glm.R`

* fixed [issue #32](https://github.com/merliseclyde/BAS/issues/32)
to allow vectorization for `phi1` function in R/cch.R
and added unit test to "tests/testthat/test-special-functions.R"
* fixed [issue #35](https://github.com/merliseclyde/BAS/issues/35) for `method="MCMC+BAS"` in `bas.glm` in `glm_mcmcbas.c` when no values are provided for `MCMC.iterations` or `n.models` and defaults are used. Added unit test in `test-bas-glm.R`

* fixed [issue #30](https://github.com/merliseclyde/BAS/issues/30) added n as hyperparameter if NULL and coerced to be a REAL for `intrinsic` prior

Expand Down
71 changes: 40 additions & 31 deletions src/glm_mcmcbas.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -197,11 +197,14 @@ SEXP glm_mcmcbas(SEXP Y, SEXP X, SEXP Roffset, SEXP Rweights,
if (nUnique < k) {
int *modelwork= ivecalloc(p);
double *pigamma = vecalloc(p);
memset(pigamma, 0.0, p*sizeof(double));

update_probs(probs, vars, mcurrent, k, p);
update_tree(modelspace, tree, modeldim, vars, k,p,n,mcurrent, modelwork);

// now sample
for (m = nUnique; m < k && pigamma[0] < 1.0; m++) {

for (m = nUnique; (m < k) && (pigamma[0] < 1.0); m++) {
for (i = n; i < p; i++) {
INTEGER(modeldim)[m] += model[vars[i].index];
}
Expand All @@ -211,46 +214,52 @@ SEXP glm_mcmcbas(SEXP Y, SEXP X, SEXP Roffset, SEXP Rweights,
bestmodel, Rparents);

/* Now subtract off the visited probability mass. */
branch=tree;
Substract_visited_probability_mass(branch, vars, model, n, m, pigamma,eps);
branch=tree;
Substract_visited_probability_mass(branch, vars, model, n, m, pigamma,eps);

/* Now get model specific calculations */
pmodel = INTEGER(modeldim)[m];
PROTECT(Rmodel_m = allocVector(INTSXP,pmodel));
GetModel_m(Rmodel_m, model, p);
pmodel = INTEGER(modeldim)[m];
PROTECT(Rmodel_m = allocVector(INTSXP,pmodel));
GetModel_m(Rmodel_m, model, p);

glm_fit = PROTECT(glm_FitModel(X, Y, Rmodel_m, Roffset, Rweights,
glm_fit = PROTECT(glm_FitModel(X, Y, Rmodel_m, Roffset, Rweights,
glmfamily, Rcontrol, Rlaplace,
betapriorfamily));
prior_m = compute_prior_probs(model,pmodel,p, modelprior);
logmargy = REAL(getListElement(getListElement(glm_fit, "lpy"),"lpY"))[0];
shrinkage_m = REAL(getListElement(getListElement(glm_fit, "lpy"),
"shrinkage"))[0];
SetModel2(logmargy, shrinkage_m, prior_m, sampleprobs, logmarg, shrinkage, priorprobs, m);
SetModel1(glm_fit, Rmodel_m, beta, se, modelspace, deviance, R2, Q, Rintercept, m);
UNPROTECT(2);
prior_m = compute_prior_probs(model,pmodel,p, modelprior);
logmargy = REAL(getListElement(getListElement(glm_fit, "lpy"),"lpY"))[0];
shrinkage_m = REAL(getListElement(getListElement(glm_fit, "lpy"),
"shrinkage"))[0];
SetModel2(logmargy, shrinkage_m, prior_m, sampleprobs, logmarg, shrinkage, priorprobs, m);
SetModel1(glm_fit, Rmodel_m, beta, se, modelspace, deviance, R2, Q, Rintercept, m);
UNPROTECT(2);

REAL(sampleprobs)[m] = pigamma[0];
REAL(sampleprobs)[m] = pigamma[0];


//update marginal inclusion probs
if (m > 1) {
double mod;
double rem = modf((double) m/(double) update, &mod);
if (rem == 0.0) {
int mcurrent = m;
compute_modelprobs(modelprobs, logmarg, priorprobs,mcurrent);
compute_margprobs(modelspace, modeldim, modelprobs, probs, mcurrent, p);
if (update_probs(probs, vars, mcurrent, k, p) == 1) {
if (m > 1) {
double mod;
double rem = modf((double) m/(double) update, &mod);
if (rem == 0.0) {
mcurrent = m;
compute_modelprobs(modelprobs, logmarg, priorprobs,mcurrent);
compute_margprobs(modelspace, modeldim, modelprobs, probs, mcurrent, p);
if (update_probs(probs, vars, mcurrent, k, p) == 1) {
// ("Updating Model Tree %d \n", m);
update_tree(modelspace, tree, modeldim, vars, k,p,n,mcurrent, modelwork);
update_tree(modelspace, tree, modeldim, vars, k,p,n,mcurrent, modelwork);
}
}
}
}
}
}
//Rprintf("Nunique = %d, m = %d, k = %d, mcurrent = %d %lf\n",
// nUnique, m, k, mcurrent, pigamma[0]);
if (m < k) {
mcurrent = m;
}
else {mcurrent = k;}
}

/* if (mcurrent < k) {
if (mcurrent < k) { // truncate vectors
SETLENGTH(modelspace, mcurrent);
SETLENGTH(logmarg, mcurrent);
SETLENGTH(modelprobs, mcurrent);
Expand All @@ -267,12 +276,12 @@ SEXP glm_mcmcbas(SEXP Y, SEXP X, SEXP Roffset, SEXP Rweights,
SETLENGTH(R2, mcurrent);
SETLENGTH(Rintercept, mcurrent);
}
*/

compute_modelprobs(modelprobs, logmarg, priorprobs, k);
compute_margprobs(modelspace, modeldim, modelprobs, probs, k, p);

INTEGER(NumUnique)[0] = k;
compute_modelprobs(modelprobs, logmarg, priorprobs, mcurrent);
compute_margprobs(modelspace, modeldim, modelprobs, probs, mcurrent, p);

INTEGER(NumUnique)[0] = mcurrent;
SET_VECTOR_ELT(ANS, 0, Rprobs);
SET_STRING_ELT(ANS_names, 0, mkChar("probne0"));

Expand Down
58 changes: 30 additions & 28 deletions tests/testthat/test-bas-glm.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -111,33 +111,6 @@ test_that("diagnostic plot for glm MCMC", {
expect_error(diagnostics(pima_MCMC, type = "pip"))
})

# FIXME issue #35
test_that("MCMC+BAS: missing MCMC.iterations and n.models arg", {
data(Pima.tr, package = "MASS")
set.seed(1)
pima_BAS <- bas.glm(type ~ .,
data = Pima.tr, method = "BAS",
betaprior = bic.prior(),
family = binomial(),
modelprior = uniform())
set.seed(1)
pima_1 <- bas.glm(type ~ ., data=Pima.tr,
method = "MCMC+BAS",
betaprior = bic.prior(),
family = binomial(),
modelprior = uniform())
set.seed(1)
pima_2 <- bas.glm(type ~ ., data = Pima.tr,
method = "MCMC+BAS",
betaprior = bic.prior(),
n.models=2^7,
MCMC.iterations=10000, #default
family = binomial(),
modelprior = uniform())
expect_equal(pima_1$probne0, pima_2$probne0)
expect_equal(pima_BAS$probne0, pima_2$probne0)
expect_equal(pima_BAS$n.models, pima_1$n.models)
})

test_that("missing data arg", {
data(Pima.tr, package = "MASS")
Expand Down Expand Up @@ -360,6 +333,7 @@ test_that("include always", {
# object 'prob' not found
})

# issue #33
test_that("Jeffreys & MCMC", {
data(Pima.tr, package="MASS")
pima_BAS <- bas.glm(type ~ .,
Expand All @@ -368,5 +342,33 @@ test_that("Jeffreys & MCMC", {
family = binomial(),
modelprior = tr.beta.binomial(1, 1, 4))
# expect_equal(0, sum(pima_BAS$probne0 > 1))
# issue #33
})

# FIXME issue #35
test_that("MCMC+BAS: missing MCMC.iterations and n.models arg", {
data(Pima.tr, package = "MASS")
set.seed(1)
pima_BAS <- bas.glm(type ~ .,
data = Pima.tr, method = "BAS",
betaprior = bic.prior(),
family = binomial(),
modelprior = uniform())
set.seed(1)
pima_1 <- bas.glm(type ~ ., data=Pima.tr,
method = "MCMC+BAS",
betaprior = bic.prior(),
family = binomial(),
modelprior = uniform())
set.seed(1)
pima_2 <- bas.glm(type ~ ., data = Pima.tr,
method = "MCMC+BAS",
betaprior = bic.prior(),
n.models=2^7,
MCMC.iterations=10000, #default
family = binomial(),
modelprior = uniform())
expect_equal(pima_1$probne0, pima_2$probne0)
expect_equal(pima_BAS$probne0, pima_2$probne0)
expect_equal(pima_BAS$n.models, pima_1$n.models)
expect_equal(pima_BAS$n.models, pima_2$n.models)
})

0 comments on commit 03d80b1

Please sign in to comment.