Nothing
EMVS.probit=function(y,x,epsilon=.0005,v0s=.025,nu.1=100,nu.gam=1,a=1,b=ncol(x),beta.initial=NULL,
sigma.initial=1,theta.inital=.5,temp=1,p=ncol(x),n=nrow(x)){ #set reasonable defauly for v0
if(length(beta.initial)==0){
beta.initial=rep(NaN,times=ncol(x))
while(sum(is.nan(beta.initial))>0)
{
beta.initial=CSDCD.logistic(p,n,x,ifelse(y==1,1,-1),rep(.0001,p),75)
}
}
scrap=0
cat("\n")
L=length(v0s)
cat("\n","Running Probit across v0's","\n")
cat(rep("",times=(L+1)),sep="|")
cat("\n")
intersects=numeric(L) # intersection points between posterior weighted spike and slab
log_post=numeric(L) # logarithm of the g-function models associated with v0s
sigma.Vec=numeric(L)
theta.Vec=numeric(L)
log_post=numeric(L)
index.Vec=numeric(L)
beta.Vec=matrix(0,L,p) # L x p matrix of MAP beta estimates for each spike
p.Star.Vec=matrix(0,L,p) # L x p matrix of conditional posterior inclusion probabilities
for (i in (1:L)){
bail.count=0
nu.0=v0s[i]
beta.Current=beta.initial
beta.new=beta.initial
sigma.EM=sigma.initial
theta.EM=theta.inital
eps=epsilon+1
iter.index=1
while(eps>epsilon && iter.index<20){ #
if(bail.count>=3)
{
cat("\n","Iteration scrapped!","\n",sep="")
list=list(betas=NULL,intersects=NULL,sigmas=NULL,
niters=NULL,posts=NULL,thetas=NULL,v0s=NULL,scrap=1)
return(list)
}
#print(iter.index)
############### E STEP #######################
d.Star=rep(NA,p)
p.Star=rep(NA,p)
for(j in 1:p){
gam.one=dnorm(beta.Current[j],0,sqrt(nu.1))**temp*theta.EM**temp
gam.zero=dnorm(beta.Current[j],0,sqrt(nu.0))**temp*(1-theta.EM)**temp
p.Star[j]=gam.one/(gam.one+gam.zero)
d.Star[j]=((1-p.Star[j])/nu.0)+(p.Star[j]/nu.1)
}
############# LATENT VARIABLE FOR PROBIT ##################
M=rep(NA,n)
y.Star=rep(NA,n)
for(j in 1:n){
M[j]=ifelse(y[j]==0,-dnorm(-x[j,]%*%beta.Current,0,1)/pnorm(-x[j,]%*%beta.Current,0,1),
dnorm(-x[j,]%*%beta.Current,0,1)/(1-pnorm(-x[j,]%*%beta.Current,0,1)))
y.Star[j]=x[j,]%*%beta.Current+M[j]
}
############### M STEP #######################
d.Star.Mat=diag(d.Star,p)
beta.old=beta.new
######## Sherman-Morrison-Woodbury Formula ####################
#beta.EM[i,]=(solve(v.Star.Mat)-(solve(v.Star.Mat)%*%t(x))%*%solve(diag(1,n)+x%*%solve(v.Star.Mat)%*%t(x))%*%(x%*%solve(v.Star.Mat)))%*%(t(x)%*%y.Star)
#beta.Current=solve(t(x)%*%x+d.Star.Mat)%*%t(x)%*%y.Star
beta.Current=CSDCD.random(p,n,x,y.Star,d.Star)
#print(beta.Current[1:10])
if(sum(is.nan(beta.Current))==0)
{
sigma.EM=1
theta.EM=(sum(p.Star)+a-1)/(a+b+p-2)
eps=max(abs(beta.new-beta.Current))
beta.new=beta.Current
iter.index=iter.index+1
}
if(sum(is.nan(beta.Current))>0){
beta.Current=beta.old
bail.count=bail.count+1
}
}
p.Star.Vec[i,]=p.Star
beta.Vec[i,]=beta.new
sigma.Vec[i]=sigma.EM
theta.Vec[i]=theta.EM
index.Vec[i]=iter.index
index=p.Star>0.5
c=sqrt(nu.1/v0s[i])
w=(1-theta.Vec[i])/theta.Vec[i]
if (w>0){
intersects[i]=sigma.Vec[i]*sqrt(v0s[i])*sqrt(2*log(w*c)*c^2/(c^2-1))}else{
intersects[i]=0}
cat("|",sep="")
}
list=list(betas=beta.Vec,intersects=intersects,sigmas=sigma.Vec,
niters=index.Vec,posts=p.Star.Vec,thetas=theta.Vec,v0s=v0s,scrap=0)
return(list)
}
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.