View source: R/simulationpipeline.R
| EnvelopeSort | R Documentation |
Sorts Enveloping function for simulation. how frequently each component of the resulting grid should be sampled during simulation.
EnvelopeSort(
l1,
l2,
GIndex,
G3,
cbars,
logU,
logrt,
loglt,
logP,
LLconst,
PLSD,
a1,
E_draws,
lg_prob_factor = NULL,
UB2min = NULL
)
l1 |
dimension for model (number of independent variables in X matrix) |
l2 |
dimension for Envelope (number of components) |
GIndex |
matrix containing information on how each dimension should be sampled (1 means left tail of a restricted normal, 2 center, 3 right tail, and 4 the entire line) |
G3 |
A matrix containing the points of tangencies associated with each component of the grid |
cbars |
A matrix containing the gradients for the negative log-likelihood at each tangency |
logU |
A matrix containing the log of the cummulative probability associated with each dimension |
logrt |
A matrix containing the log of the probability associated with the right tail (i.e. that to the right of the lower bound) |
loglt |
A matrix containing the log of the probability associated with the left tail (i.e., that to the left of the upper bound) |
logP |
A matrix containing log-probabilities related to the components of the grid |
LLconst |
A vector containing constant for each component of the grid used during the accept-reject procedure |
PLSD |
A vector containing the probability of each component in the Grid |
a1 |
A vector containing the diagonal of the standardized precision matrix |
E_draws |
Bound on Expected number of candidates per accepted draw |
lg_prob_factor |
vector of lg_prob_factors used for the Envelope connected to the independent normal gamma prior |
UB2min |
Vector containing min for UB2 for each component (relevant for EnvelopeDispersionBuild) |
This function sorts the envelope in descending order based on the
probability associated with each component in the Grid. Sorting helps
speed up simulation once the envelope is constructed. If memory allocation
fails (e.g. for very large grids), the function returns the unsorted envelope
with sort_ok = FALSE; the sampler remains valid but may have poorer acceptance.
Used after EnvelopeBuild and (for Normal–Gamma models)
EnvelopeDispersionBuild; see \insertCiteNygren2006,glmbayesChapterA08glmbayes.
The function returns a list consisting of the following components (the first six of which are matrics with number of rows equal to the number of components in the Grid and columns equal to the number of parameters):
GridIndex |
A matrix containing information on how each dimension should be sampled (1 means left tail of a restricted normal, 2 center, 3 right tail, and 4 the entire line) |
thetabars |
A matrix containing the points of tangencies associated with each component of the grid |
cbars |
A matrix containing the gradients for the negative log-likelihood at each tangency |
logU |
A matrix containing the log of the cummulative probability associated with each dimension |
logrt |
A matrix containing the log of the probability associated with the right tail (i.e. that to the right of the lower bound) |
loglt |
A matrix containing the log of the probability associated with the left tail (i.e., that to the left of the upper bound) |
LLconst |
A vector containing constant for each component of the grid used during the accept-reject procedure |
logP |
A matrix containing log-probabilities related to the components of the grid |
PLSD |
A vector containing the probability of each component in the Grid |
E_draws |
A containing a computed theoretical bound on the expected number of draws |
sort_ok |
Logical; |
EnvelopeBuild, EnvelopeOrchestrator,
EnvelopeDispersionBuild, rNormal_reg, rglmb.
data(menarche,package="MASS")
Age2=menarche$Age-13
summary(menarche)
plot(Menarche/Total ~ Age, data=menarche)
x<-matrix(as.numeric(1.0),nrow=length(Age2),ncol=2)
x[,2]=Age2
y=menarche$Menarche/menarche$Total
wt=menarche$Total
mu<-matrix(as.numeric(0.0),nrow=2,ncol=1)
mu[2,1]=(log(0.9/0.1)-log(0.5/0.5))/3
V1<-1*diag(as.numeric(2.0))
# 2 standard deviations for prior estimate at age 13 between 0.1 and 0.9
## Specifies uncertainty around the point estimates
V1[1,1]<-((log(0.9/0.1)-log(0.5/0.5))/2)^2
V1[2,2]=(3*mu[2,1]/2)^2 # Allows slope to be up to 1 times as large as point estimate
famfunc<-glmbfamfunc(binomial(logit))
f1<-famfunc$f1
f2<-famfunc$f2
f3<-famfunc$f3
f5<-famfunc$f5
f6<-famfunc$f6
dispersion2<-as.numeric(1.0)
start <- mu
offset2=rep(as.numeric(0.0),length(y))
P=solve(V1)
n=1000
## Appears that the type for some of these arguments are important/problematic
### This example constructs and sorts an envelope without calling the
### lower-level sampler directly.
###### Adjust weight for dispersion
wt2=wt/dispersion2
######################### Shift mean vector to offset so that adjusted model has 0 mean
alpha=x%*%as.vector(mu)+offset2
mu2=0*as.vector(mu)
P2=P
x2=x
##### Optimization step to find posterior mode and associated Precision
parin=start-mu
opt_out=optim(parin,f2,f3,y=as.vector(y),x=as.matrix(x),mu=as.vector(mu2),
P=as.matrix(P),alpha=as.vector(alpha),wt=as.vector(wt2),
method="BFGS",hessian=TRUE
)
bstar=opt_out$par ## Posterior mode for adjusted model
bstar
bstar+as.vector(mu) # mode for actual model
A1=opt_out$hessian # Approximate Precision at mode
## Standardize Model
Standard_Mod=glmb_Standardize_Model(y=as.vector(y), x=as.matrix(x),
P=as.matrix(P),bstar=as.matrix(bstar,ncol=1), A1=as.matrix(A1))
bstar2=Standard_Mod$bstar2
A=Standard_Mod$A
x2=Standard_Mod$x2
mu2=Standard_Mod$mu2
P2=Standard_Mod$P2
L2Inv=Standard_Mod$L2Inv
L3Inv=Standard_Mod$L3Inv
Env2=EnvelopeBuild(as.vector(bstar2), as.matrix(A),y, as.matrix(x2),
as.matrix(mu2,ncol=1),as.matrix(P2),as.vector(alpha),as.vector(wt2),
family="binomial",link="logit",Gridtype=as.integer(3),
n=as.integer(n),sortgrid=FALSE)
## Extract l1, l2 from envelope (as done in C++) and call EnvelopeSort
l1 <- ncol(Env2$cbars)
l2 <- nrow(Env2$cbars)
logP_mat <- matrix(Env2$logP, ncol = 1)
Env_sorted <- EnvelopeSort(l1, l2,
GIndex = Env2$GridIndex,
G3 = Env2$thetabars,
cbars = Env2$cbars,
logU = Env2$logU,
logrt = Env2$logrt,
loglt = Env2$loglt,
logP = logP_mat,
LLconst = Env2$LLconst,
PLSD = Env2$PLSD,
a1 = Env2$a1,
E_draws = Env2$E_draws
)
Env_sorted
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.