tlnise_function(Y,V,W="NA",V0="NA",Vs="NA",P=1,N=1000, alpha=.05,
 Astart="NA",ds="NA",nuo="NA",Tol=1e-6,maxiter=1000,intercept=T, 
 labelW="NA",labelY="NA",digits=3,srt=F,brief=0,prnt=T){
if(!is.loaded(Fortran.symbol("mode2dfb0u"))||
   !is.loaded(Fortran.symbol("postu")))
# 			!!!IMPORTANT!!! 
# you must change the path in the following line to specify the directory in
# which you stored the Fortran object file tlnise.o.
 dyn.load("/home/everson/Frtrn/Nise/TLNise/tlnise.o")
#          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -> change to agree with your path
#
# TLNise: Two-Level Normal independent sampling estimation 
#  Provides inference for 2-level Normal hierarchical 
#  models with UNIVARIATE (p = 1) outcomes. 
# Calls Fortran subroutine modefB0u.f to find posterior 
#  mode for sampling variable B0 = V0/(V0+A).
# Calls Fortran subroutine postu.f to process draws from
#  approximating truncated Gamma distribution and provide
#  posterior estimates of the r-dimensional level-two 
#  regression coefficient Gamma, the level-two variance
#  component A, and of the level-one mean parameters 
#  theta1-thetaJ. 
#
# inputs:
#        Y: J-vector of Normal outcomes
#        V: J-vector of Level-1 variances
#        W: Jxr covariate matrix (adds column of 1's if not
#	      included and intercept = T)
#       V0: "typical" Vj (default is average of Vj's)
#       Vs: V* for defining prior on B* = V*/(V*+A)
#        P: prior specification; 
#         P = 0: uniform on B0; 
#	  P=0.5: Jeffreys' prior;
#	  P = 1: uniform on A;
#         P = 2: uniform on B* = V*/(V*+A) (default V* = V0)
#	  P = 3: p(B*) = B*^(ds/2)         (default ds = -2)
#        N: number of truncated Gamma draws for inference
#    alpha: specifies (1-alpha)100% confidence level (default=0.05, for 95%)
#   Astart: starting value for locating mode of B0
#	ds: required for prior option P=3
#      nu0: alternate value of df - default is J-r (nu0 <= J-r)
#      Tol: tolorance for determining modal convergence
#  maxiter: maximum number of EM iterations for finding mode
#intercept: if True, an intercept term is included in the regression
#   labelW: optional names vector for covariates
#   labelY: optional names vector for observations
#   digits: number of digits beyond decimal for reporting results
#      srt: if True, sorts observations on V
#    brief: level of output, from 0 (minimum) to 2 (maximum)
#     prnt: Controls printing.
#
#
J_length(Y)
if(missing(W))
 W_rep(1,J)
W_as.matrix(W)
if(intercept&&max(abs(W[,1]-1))!=0) W_cbind(rep(1,J),W)
r_dim(W)[[2]]
Y_matrix(Y,ncol=1)
if(missing(V0)) V0_mean(V)
if(missing(Astart)) Astart_V0
#
# determine prior df nu*
#CASE P = 0: flat prior distribution assumed for B0 (find MLE of A)
if(P == 0) {
 if(prnt) print("flat prior distribution assumed for B0", quote = F)
 ds <- 0
 Vs <- V0
 }
if(P==0.5) {
 if(prnt) print("Jeffreys' prior distribution assumed for B0", quote = F)
 ds_-2
 Vs_V0
 }
#  CASE P = 1: flat prior distribution assumed for A
if(P == 1) {
 if(prnt) print("flat prior distribution assumed for A", quote = F)
 ds <- -4 
 Vs <- V0
 }
#  Case P = 2: flat prior distribution assumed for Bs
if(P == 2) {
 if(prnt) print("flat prior distribution assumed for B*", quote = F)
 if(missing(Vs)) {
  print("Vs is MISSING; set equal to V0", quote = F)
  Vs <- V0
  }
 ds <- 0
 }
# Case P = 3: p(Bs) = |Bs|^(ds/2)
if(P == 3) {
 if(prnt) print("prior on Bs defined by input ds: p(Bs)=Bs^(ds/2)", quote = F)
 if(missing(ds)) {
  print("missing ds; set equal to -2 (flat prior distribution for B*)", 
	quote = F)
  ds <- -2
  }
 if(missing(Vs)) {
  if(prnt) print("Vs IS MISSING; set equal to V0", quote = F)
  Vs <- V0
   }
  }
#
# determine df nu for truncated Gamma(nu/2,Sig/2;1)
if(missing(nuo)) nuo_J-r
nus_ds+2
nu_nuo+nus
# check for "nearly equal variance" case:
if((max(V)/min(V))<2&missing(Astart)){
 mu_W%*%solve(t(W)%*%W)%*%t(W)%*%Y
 Astart_max(Tol,sum((Y-mu)^2)/(J-r+ds) - mean(V))}
#
# find posterior mode of B0=V0/(V0+A) via EM algorithm
#
if(prnt) print("********* Locating posterior mode of B0 ************",quote=F)
#
out.mode_niseEMu.v1.f(Astart,Y,V,V0,Vs,W,t(W),ds,J,r,Tol,maxiter)
modeA_out.mode$A
modeB0_out.mode$modeB0
lf1mode_out.mode$lf
d1_out.mode$d1
if(d1>.5){
 if(prnt) print("matching derivatives at modal value", quote=F)
# boundary mode - match derivatives at mode
 S0_(nu-2)/modeB0-2*d1
 }
else{
 S0_(nu-2)/modeB0
 }
lf0mode_log(dtgam(modeB0,nu/2,S0/2,1))
#}
adj_lf0mode-lf1mode
cat("\n")
if(prnt) print("******** Drawing truncated Gammas *********** ",quote=F)
cat("\n")
B0draws_-sort(-rtgam(N,nu/2,S0/2,1))
Adraws_V0*(1/B0draws -1)  # ordered smallest to largest
if(prnt) print("************ Processing draws *************** ",quote=F)
f0_dtgam(B0draws,nu/2,S0/2,1)
quant_c(0.05,0.10,0.25,0.5,0.75,0.90,0.95)
out.post_postu.f(Adraws,N,ds,f0,adj,Y,V,V0,Vs,W,t(W),J,r,quant,alpha)
#
ratio_out.post$rat/out.post$maxrat
if(prnt) print(paste("Average Importance Weight =", round(1000*mean(ratio))/1000),quote=F)
mlt_10^digits
if(prnt) print("Estimated quantiles of f(A|Y):",quote=F)
Aqnt_round(mlt*out.post$Aqnt)/mlt
B0qnt_round(mlt*(V0/(V0+Aqnt))[7:1])/mlt
names(Aqnt)_names(B0qnt)_c("5%","10%","25%","Median","75%","90%","95%")
if(prnt){
 print(Aqnt)
 print("sqrt(A):",quote=F)
 print(round(mlt*sqrt(Aqnt))/mlt)
 print("Estimated quantiles of f(B0|Y):",quote=F)
 print(paste("B0 = ",signif(V0,4),"/(",signif(V0,4),"+ A ) "),quote=F)
 print(B0qnt)
 }
gammat_cbind(out.post$gammap,
#out.post$gamlp,out.post$gamup,
	sqrt(diag(out.post$Dgp)),out.post$gammap/sqrt(diag(out.post$Dgp)))
if(missing(labelW)) labelW_c(0:(r-1))
if(length(labelW)==(r-1)) labelW_c("intrcpt",labelW)
dimnames(gammat)_list(labelW,c("est",
#"lower","upper",
"sd","ratio"))
gammat_apply(mlt*gammat,2,round)/mlt
if(r==1){
 theta.mat_cbind(Y,sqrt(V),W*gammat[1],V/(V+Aqnt[4]),out.post$thetap,
 ,out.post$thetlp,out.post$thetup,sqrt(out.post$Vp))}
else{
 theta.mat_cbind(Y,sqrt(V),W%*%gammat[,1],V/(V+Aqnt[4]),out.post$thetap,
# out.post$thetlp,out.post$thetup,
sqrt(out.post$Vp))}
theta.mat_apply(mlt*theta.mat,2,round)/mlt
if(missing(labelY)) labelY_1:J
dimnames(theta.mat)_list(labelY,c("Y","rt(V)","mu","B","theta*",
#"lower","upper",
"rt(Vtheta)"))
if(srt)
 theta.mat_theta.mat[order(-V),]
if(brief==0)
 out_list(gamma=gammat,theta=theta.mat)
if(brief==1)
 out_list(gamma=gammat,theta=theta.mat,A=Aqnt,modeB0=modeB0,param=c(nu=nu,S0=S0),
  Adraws=out.post$Adraws,ratio=ratio)
if(brief==2)
 out_list(gamma=gammat, theta=theta.mat, A=Aqnt, modeB0=modeB0, 
	param=c(nu=nu,S0=S0), Adraws=out.post$Adraws, ratio=ratio, 
	thetam=out.post$thetam,Vm=out.post$Vm, gammam=out.post$gammam, Dga=out.post$Dga, f1=out.post$f1, f0=f0)
if(brief==3)
 out_list(out.mode=out.mode,out2=c(nu,S0),out.post=out.post,lf=lf1mode,adj=adj)
out}

niseEMu.v1.f_function(Astart,Y,V,V0,Vs,W,Wt,ds,J,r,tol,maxit){
#
# calls Fortran subroutine "mode2DfB0u" to find the posterior 
# mode of B0 = V0/(V0+A) assuming p(Bs) = Bs^{ds/2}
#
newA_lf_lA_iter_0
gam_rep(0,r)
rt_rep(0,3)
d1_0
ratevec_rep(0,500)
Dgam_matrix(0,ncol=r,nrow=r)
out_.Fortran("mode2dfb0u",
	as.double(Astart),
	as.double(Y),
	as.double(V),
	as.double(V0),
	as.double(Vs),#5
	as.double(W), 
	as.double(Wt),
	as.integer(ds),
	as.double(ds),
	as.integer(J), #10
	as.integer(r), 
	as.double(tol),
	as.integer(maxit),
	newA=newA,
	lf=lf,
	lA=lA, #16
	as.integer(0),  # iter
	d1=d1,
	as.double(V),    #deltinv 19
	as.double(V),    #B
	as.double(Dgam), #Dgam
	as.double(gam),  #WDY 
	as.double(gam),  #gam 
	as.double(V),    #mu  24
	as.double(V0),   #H
	as.double(V),    #us
	as.double(V),    #Dus 
	as.double(lf),   #ttol 
	as.double(lf),   #ldDg 29
	as.double(Dgam), #temprr
	as.double(gam),  #tempr
	as.double(W),    #tempJr
	as.double(V),    #tempJ
	as.double(rt))
   iter_out[[17]]
   if(iter == maxit)
     print(paste("did not converge within tolerance in", maxit,"iterations"), quote = F)
   else{
    print(paste("converged to within tolerance in", iter-1,"EM iterations"), 
	quote = F)
	}
   A_out$newA
   modeB0_V0/(V0+A)
   print(paste("posterior mode for B0 =",round(1000*modeB0)/1000,"( A =",round(100*A)/100,")"),quote=F)
   lA_out$lA
   lf_out$lf
   d1_out$d1
   d2_out$d2
   ratevec_out$ratevec
   list(modeB0=modeB0,A=A,lf=lf,lA=lA,d1=d1)
}

postu.f_function(Adraws,N,ds,f0,adj,Y,V,V0,Vs,W,Wt,J,r,quant,alpha){
#
# calls Fortran subroutine "postu" to processes UNIVARIATE B0 draws
#from an approximating density f0, calculating the true posterior
#density f(B0|Y) for a prior distribution p(B*) = B*^(d*/2); B* =
#V*/(V*+A), the normailzed importance weights f(B0|Y)/f0(B0), and the 
#conditional means and variances of the thetas and of gamma given each
#draw B0.
#
zstar_qnorm(1-alpha/2)
f1_rat_rep(0,N)
thetam_Vm_matrix(0,ncol=N,nrow=J)
gammam_matrix(0,ncol=N,nrow=r)
Dga_array(0,c(r,r,N))
gammap_gamup_gamlp_rep(0,r)
thetap_Vp_thetlp_thetup_rep(0,J)
Dgp_matrix(0,ncol=r,nrow=r)
Ap_0
Aqnt_quant
maxrat_sumrat_1
Dgam_matrix(0,ncol=r,nrow=r)
   out_.Fortran("postu",
	as.double(Adraws),
	as.integer(N),
	as.integer(ds),
	as.double(ds),
	as.double(f0),
	as.double(adj), 	# 6
	as.double(Y),
	as.double(V),   
	as.double(V0),
	as.double(Vs),
	as.double(W),		# 11
	as.double(Wt),   
	as.integer(J),
	as.integer(r),
	as.double(quant),
	as.double(zstar),
	thetap=thetap,		# 17
	thetup=thetup,
	thetlp=thetlp,
	Vp=Vp,
	gammap=gammap,
	gamlp=gamlp,
	gamup=gamup,
	Dgp=Dgp,	
	Ap=Ap,
	Aqnt=Aqnt,		# 25
	f1=f1,    
	rat=rat, 		   
	thetam=thetam, 
	Vm=Vm,
	gammam=gammam,		# 30
	Dga=Dga, 
	maxrat=maxrat,		
	sumrat=sumrat,		
	as.double(V),     	# B 	
	as.double(V),		# deltinv 32
	as.double(gammap),	# WDY
	as.double(gammap),	# gam    
	as.double(Dgp),		# Dgam  		   
	as.double(V),		# mu		
	as.double(gammap),	# tempr   37
	as.double(Dgp),		# temprr
	as.double(W),		# tempJr
	as.double(V)		# tempJ  
	)
	out
}