Revision c04c770bc32df0fa9b80e8783414392020eba35a authored by Wesley Tansey on 08 July 2016, 19:53:59 UTC, committed by Wesley Tansey on 08 July 2016, 19:53:59 UTC
1 parent 8598e9e
braindt.cpp
#include"braindt.h"
//#include <iostream>
//#include<vector>
#include <cfloat>//for IsFiniteNumber and IsNumber function
#include <algorithm>
#include<cstdio>
#include<cstdlib>
#include<string>
//#include<cstring>
#define _USE_MATH_DEFINES
#include<math.h>
using namespace std;
const double braindt::TINY=1e-20;
const double braindt::mu0=0;
const double braindt::sigma0_sq=1;
//constructor function of class braindt
braindt::braindt(int &LX, int &LY, int &LZ,int &Num,char *argv ,vector<int> &Loca,int &PN,vector<double>& est,int &L_hat,int sweep_b, int sweep_r,int sweep_lis,int iter_max):\
Lx(LX),Ly(LY),Lz(LZ),N(Num),pN(PN),y(PN),L(L_hat),U(2),invI(3),I(3),swp_b(sweep_b),swp_r(sweep_r),max_iter(iter_max),initial(PN ),delta(2),\
H_mean0(2),H0(sweep_r,vector<double>(2)),H_cond_mean(2),mu1_hat_pre(L_hat ),sigma1_sq_hat_pre(L_hat ),pEll_hat_pre(L_hat ),LIS(PN),swp_lis(sweep_lis)
{
vector<double> y_ori;//observed data on the whole image grid including the brain part and the non-brain part
ifstream ifs(argv);//
double temp_y;
while(ifs>>temp_y)
y_ori.push_back(temp_y);// read in y_ori
int i;
loca=Loca.begin();//the location of the brain region of interest (ROI) on the Lx*Ly*Lz image grid
for(int pi=0;pi<pN;pi++)
{
i=loca[pi];
y[pi]=y_ori[i];//the observed data of the ROI
}
YbyX=Ly*Lx;
beta_hat=est[0];
h_hat=est[1];
int i3;
for(int i=0;i<L;i++)
{
i3=i*3;
mu1_hat.push_back(est[2+i3]);
sigma1_sq_hat.push_back(est[3+i3]);
pEll_hat.push_back(est[4+i3]);
}
}
inline double braindt::normpdf(double &x, double &mu, double &sigma_sq)// The probability density function of the normal distribution with mean=mu, variance=sigma_sq
{
return exp(-1*(x-mu)*(x-mu)/(2*sigma_sq))/(sqrt(2*M_PI*sigma_sq));
}
inline bool braindt::IsFiniteNumber(double &x)// Judge whether x is a finite numebr
{
return (x<=DBL_MAX && x>=-DBL_MAX);
}
inline bool braindt::IsNumber(double &x)//Judge whether x is a numebr
{
return (x==x);
}
double braindt::sigmasum(vector<double> &p, int begin, int end)//sum of elements of vector p
{
double sum=0;
for(int i=begin; i<end; i++)
{
sum+=p[i];
}
return sum;
}
double braindt::sigmasum(vector<double> &p1, vector<double> &p2, int begin, int end)//sum of elements of the entrywise-product of vectors p1 and p2
{
double sum=0;
for(int i=begin; i<end; i++)
{
sum+=p1[i]*p2[i];
}
return sum;
}
double braindt::sigmasum(vector<double> &p1, vector<double> &p2, double p2minus, int begin, int end)//sum of elements of the entrywise-product of vectors p1 and p3 squared, where p3=p2-p2minus
{
double sum=0;
double minus;
for(int i=begin; i<end; i++)
{
minus=p2[i]-p2minus;
sum+=p1[i]*minus*minus;
}
return sum;
}
void braindt::two_d_var_matrix_inv(ofstream &fileout, int splnumber)// compute the inverse matrix of a 2 by 2 matrix
{
double s1,s2;
vector<double> var(3);
for(int i=0;i<splnumber;i++)
{
s1=H[i][0]-H_mean[0];
var[0]+=s1*s1;
s2=H[i][1]-H_mean[1];
var[2]+=s2*s2;
var[1]+=s1*s2;//non-diagonal entry
}
splnumber-=1;
var[0]/=splnumber;
var[2]/=splnumber;
var[1]/=splnumber;
for(int t=0;t<3;t++)
{fileout<<"I["<<t<<"]="<<var[t]<<endl;}
double denom=var[0]*var[2]-var[1]*var[1];
if(denom>0)
{
invI[0]=var[2]/denom;
invI[2]=var[0]/denom;
invI[1]=-var[1]/denom;
if(!(IsFiniteNumber(invI[0])&&IsFiniteNumber(invI[1])&&IsFiniteNumber(invI[2])&&IsNumber(invI[0])&&IsNumber(invI[1])&&IsNumber(invI[2])))
{
fileout<<"invI is almost non-invertible matrix"<<endl;
var[0]+=TINY;
var[2]+=TINY;
denom=var[0]*var[2]-var[1]*var[1];
invI[0]=var[2]/denom;
invI[2]=var[0]/denom;
invI[1]=-var[1]/denom;
if(!(IsFiniteNumber(invI[0])&&IsFiniteNumber(invI[1])&&IsFiniteNumber(invI[2])&&IsNumber(invI[0])&&IsNumber(invI[1])&&IsNumber(invI[2])))
{
fileout<<"we choose invI=I"<<endl;
invI[0]=1;
invI[1]=0;
invI[2]=1;
}
}
}
else {
fileout<<"invI is non-invertible matrix"<<endl;
var[0]+=TINY;
var[2]+=TINY;
denom=var[0]*var[2]-var[1]*var[1];
invI[0]=var[2]/denom;
invI[2]=var[0]/denom;
invI[1]=-var[1]/denom;
if(!(IsFiniteNumber(invI[0])&&IsFiniteNumber(invI[1])&&IsFiniteNumber(invI[2])&&IsNumber(invI[0])&&IsNumber(invI[1])&&IsNumber(invI[2])))
{
fileout<<"we choose invI=I"<<endl;
invI[0]=1;
invI[1]=0;
invI[2]=1;
}
}
I=var;
}
void braindt::two_d_var_matrix_inv( int splnumber)// compute the inverse matrix of a 2by2 matrix
{
double s1,s2;
vector<double> var(3);
for(int i=0;i<splnumber;i++)
{
s1=H[i][0]-H_mean[0];
var[0]+=s1*s1;
s2=H[i][1]-H_mean[1];
var[2]+=s2*s2;
var[1]+=s1*s2;//non-diagonal entry
}
splnumber-=1;
var[0]/=splnumber;
var[2]/=splnumber;
var[1]/=splnumber;
double denom=var[0]*var[2]-var[1]*var[1];
if(denom>0)
{
invI[0]=var[2]/denom;
invI[2]=var[0]/denom;
invI[1]=-var[1]/denom;
if(!(IsFiniteNumber(invI[0])&&IsFiniteNumber(invI[1])&&IsFiniteNumber(invI[2])&&IsNumber(invI[0])&&IsNumber(invI[1])&&IsNumber(invI[2])))
{
var[0]+=TINY;
var[2]+=TINY;
denom=var[0]*var[2]-var[1]*var[1];
invI[0]=var[2]/denom;
invI[2]=var[0]/denom;
invI[1]=-var[1]/denom;
if(!(IsFiniteNumber(invI[0])&&IsFiniteNumber(invI[1])&&IsFiniteNumber(invI[2])&&IsNumber(invI[0])&&IsNumber(invI[1])&&IsNumber(invI[2])))
{
invI[0]=1;
invI[1]=0;
invI[2]=1;
}
}
}
else {
var[0]+=TINY;
var[2]+=TINY;
denom=var[0]*var[2]-var[1]*var[1];
invI[0]=var[2]/denom;
invI[2]=var[0]/denom;
invI[1]=-var[1]/denom;
if(!(IsFiniteNumber(invI[0])&&IsFiniteNumber(invI[1])&&IsFiniteNumber(invI[2])&&IsNumber(invI[0])&&IsNumber(invI[1])&&IsNumber(invI[2])))
{
invI[0]=1;
invI[1]=0;
invI[2]=1;
}
}
I=var;
}
void braindt::gem(ofstream &fileout,char* filename)//The Generalized EM Algorithm
{
int i,j,sum,sum_cond,sweep,ell,pi;
int t;///////
fileout<<"The beginning of GEM with initial values:"<<endl;
fileout<<"beta_hat="<<beta_hat<<endl;
fileout<<"h_hat="<<h_hat<<endl;
for(ell=0;ell<L;ell++)
{
fileout<<"mu1_hat["<<ell<<"]="<<mu1_hat[ell]<<endl;
fileout<<"sigma1_sq_hat["<<ell<<"]="<<sigma1_sq_hat[ell]<<endl;
fileout<<"pEll_hat["<<ell<<"]="<<pEll_hat[ell]<<endl;
}
fileout<<"End the beginning of GEM"<<endl<<endl;
double TINY_L=TINY/L;
int error_time=0;
double temp_like;
int seed=0;
CRandomMersenne rnd1(seed);
StochasticLib1 rnd2(seed);
double a,b;//for penalized MLE (PMLE) when L>=2
a=1;
b=2;
double del=1e-3;
int succ_t=0;// to see the stopping rule successive time in the t iteration
double exp_x[7];
vector<double> exp_y(pN );
double beta_hat_pre,h_hat_pre;//previous values for t iteration, i.e., the interations of GEM algorithm
double d1=1;//error in t iteration of GEM
vector<double> diff(2+3*L);
double p_1;
double temp2;
vector<double> constant(L );
vector<int> x(N );//x is the unobservable state 0 or 1
vector<int> x_cond(N );// conditional x on y
double constant_0=sqrt(2*M_PI*sigma0_sq);
vector<double> temp1(pN);
for(pi=0;pi<pN ;pi++)
{
initial[pi]=rnd1.Random()<0.5?1:0;
temp1[pi]=constant_0*exp(pow(y[pi]-mu0,2)/(2*sigma0_sq));
}
vector<double> log_q2_sub(swp_r); // -phi^T*H(x)
alam=0;
int n_tout=1;
for(t=0;n_tout && t<max_iter;t++)
{
fileout<<"The GEM iteration t="<<t<<endl;
beta_hat_pre=beta_hat;
h_hat_pre=h_hat;
vector<double> gamma_1(pN);
for(i=0;i<L;i++)
{
constant[i]=pEll_hat[i]/sqrt(2*M_PI*sigma1_sq_hat[i]);
}
for(pi=0;pi<pN;pi++)
{
temp2=0;
i=loca[pi];
x[i]=x_cond[i]=initial[pi];
for(j=0;j<L;j++)
{
temp2+=constant[j]*exp(-pow(y[pi]-mu1_hat[j],2)/(2*sigma1_sq_hat[j]));
}
exp_y[pi]=temp1[pi]*temp2;
}
for(i=0; i<7; i++)
{
exp_x[i]=exp(beta_hat*i+h_hat);
}
H_cond_mean=H_mean0;//set to be zeros
for(sweep=1;sweep<swp_b;sweep++)
{
for(pi=0; pi<pN;pi++)
{
i=loca[pi];
sum=x[i+1]+x[i-1]+x[i+Lx]+x[i-Lx]+x[i+YbyX]+x[i-YbyX];
sum_cond=x_cond[i+1]+x_cond[i-1]+x_cond[i+Lx]+x_cond[i-Lx]+x_cond[i+YbyX]+x_cond[i-YbyX];
p_1=1-1/(1+exp_x[sum]);
x[i]=rnd2.Bernoulli(p_1);
p_1=1-1/(1+exp_x[sum_cond]*exp_y[pi]);
x_cond[i]=rnd2.Bernoulli(p_1);
}
}
H_mean=H_mean0;//set to be zeros
H=H0;
for(sweep=0;sweep<swp_r;sweep++)
{
for(pi=0; pi<pN;pi++)
{
i=loca[pi];
sum=x[i+1]+x[i-1]+x[i+Lx]+x[i-Lx]+x[i+YbyX]+x[i-YbyX];
sum_cond=x_cond[i+1]+x_cond[i-1]+x_cond[i+Lx]+x_cond[i-Lx]+x_cond[i+YbyX]+x_cond[i-YbyX];
p_1=1-1/(1+exp_x[sum]);
x[i]=rnd2.Bernoulli(p_1);
p_1=1-1/(1+exp_x[sum_cond]*exp_y[pi]);
x_cond[i]=rnd2.Bernoulli(p_1);
}
for(pi=0; pi<pN;pi++)
{
i=loca[pi];
sum=x[i+1]+x[i-1]+x[i+Lx]+x[i-Lx]+x[i+YbyX]+x[i-YbyX];
sum_cond=x_cond[i+1]+x_cond[i-1]+x_cond[i+Lx]+x_cond[i-Lx]+x_cond[i+YbyX]+x_cond[i-YbyX];
H_cond_mean[0]+=x_cond[i]*sum_cond;//warning:1/2 in the later, for compute twice//Don't forgot to divide swp_r at the end
H_cond_mean[1]+=x_cond[i];
H[sweep][0]+=x[i]*sum*0.5;
H[sweep][1]+=x[i];
gamma_1[pi]+=x_cond[i];//r_[i][1]
//////////////////////////////////////////////////////////////
}
log_q2_sub[sweep]=-H[sweep][0]*beta_hat-H[sweep][1]*h_hat;
//////////////////////////////////////////////////////////////////
H_mean[0]+=H[sweep][0];//don't forgot to divide swp_r at the end
H_mean[1]+=H[sweep][1];
}
///////////////////////////////////////////////////////////////////////
temp_like=*max_element(log_q2_sub.begin(), log_q2_sub.end());
log_plikelihood=0;
for(sweep=0;sweep<swp_r;sweep++)
{
log_plikelihood+=exp(log_q2_sub[sweep]-temp_like);
}
q2=temp_like+log(log_plikelihood);//-logZ
//fileout<<"q2="<<q2<<endl;
///////////////// U
H_cond_mean[0]/=2;
fileout<<"The score function U:"<<endl;
for(i=0;i<2;i++)
{
H_cond_mean[i]/=swp_r;
H_mean[i]/=swp_r;
U[i]=H_cond_mean[i]-H_mean[i];
fileout<<"U["<<i<<"]="<<U[i]<<endl;
fileout<<"H_cond_mean["<<i<<"]="<<H_cond_mean[i]<<endl;
fileout<<"H_mean["<<i<<"]="<<H_mean[i]<<endl;
}
fileout<<"The information matrix I:"<<endl;
two_d_var_matrix_inv(fileout,swp_r);// for I update
//////////////////////////revised by 20140219
for(pi=0;pi<pN;pi++)
{
gamma_1[pi]=gamma_1[pi]/swp_r;
}
double gamma_1_sum=sigmasum(gamma_1,0,pN);
mu1_hat_pre=mu1_hat;
sigma1_sq_hat_pre=sigma1_sq_hat;
pEll_hat_pre=pEll_hat;
if(L==1)
{
pEll_hat[0]=1;
mu1_hat[0]=(sigmasum(gamma_1,y,0,pN)+TINY)/(gamma_1_sum+TINY);//+TINY to avoid the case gamma_1_sum=0
sigma1_sq_hat[0]=(sigmasum(gamma_1,y,mu1_hat[0],0,pN )+TINY)/(gamma_1_sum+TINY);
}
else
{
vector<double> f1y(pN);
vector<double> omega_sum(L);
vector<vector<double> > omega(L,vector<double>(pN));
for(pi=0;pi<pN;pi++)
{
for(ell=0;ell<L;ell++)
{
f1y[pi]+=pEll_hat[ell]*normpdf(y[pi], mu1_hat[ell], sigma1_sq_hat[ell]);
}
for(ell=0;ell<L;ell++)
{
omega[ell][pi]=gamma_1[pi]*pEll_hat[ell] \
*normpdf(y[pi],mu1_hat[ell],sigma1_sq_hat[ell])/f1y[pi];
omega_sum[ell]+=omega[ell][pi];
}
}
for(ell=0;ell<L;ell++)
{
mu1_hat[ell]=(sigmasum(omega[ell],y,0,pN )+TINY)/(omega_sum[ell]+TINY);// +TINY to avoid omega_sum[ell]=0
sigma1_sq_hat[ell]=(sigmasum(omega[ell],y,mu1_hat[ell],0,pN )+2*a)/(omega_sum[ell]+2*b);//PMLE
if(ell==(L-1))
{
pEll_hat[ell]=0;
for(int ell_sub=0;ell_sub<ell;ell_sub++)
{
pEll_hat[ell]+=pEll_hat[ell_sub];
}
pEll_hat[ell]=1-pEll_hat[ell];
///////if pELL_hat[L-1]<0, because of precision problem, to be like -2.22045e-16
if(pEll_hat[ell]<0)
{
pEll_hat[(ell-1)]+=pEll_hat[ell];
pEll_hat[ell]=0;
}
}
else {
pEll_hat[ell]=(omega_sum[ell]+TINY_L)/(gamma_1_sum+TINY);//+TINY to avoid gamma_1_sum=0
// fileout<<"omega_sum["<<ell<<"]="<<omega_sum[ell]<<endl;
// fileout<<"gamma_1_sum="<<gamma_1_sum<<endl;
}
}
}
/*
fileout<<"before backtrack"<<endl;
fileout<<"beta_hat="<<beta_hat<<endl;
fileout<<"h_hat="<<h_hat<<endl;
fileout<<"U[0]="<<U[0]<<endl;
fileout<<"U[1]="<<U[1]<<endl;
*/
fileout<<"Go to backtracking linear seach"<<endl;
backtrack(fileout,rnd2);
fileout<<"Back to GEM from backtracking linear seach"<<endl;
//m iteration is over
///////////////////////
diff[0]=fabs(beta_hat_pre-beta_hat)/(fabs(beta_hat_pre)+del);
diff[1]=fabs(h_hat_pre-h_hat)/(fabs(h_hat_pre)+del);
fileout<<"The estimates obtained in the iteration of GEM t="<<t<<":"<<endl;
fileout<<"beta_hat="<<beta_hat<<endl;
fileout<<"h_hat="<<h_hat<<endl;
for(ell=0;ell<L;ell++)
{
if(pEll_hat_pre[ell]<1e-3 && pEll_hat[ell]<1e-3)
{
fileout<<"pEll_hat_pre and pEll_hat are too close to 0 to be accounted into d1"<<endl;
diff[2+ell]=0;
diff[2+L+ell]=0;
diff[2+2*L+ell]=0;
}
else
{
diff[2+ell]=fabs(mu1_hat_pre[ell]-mu1_hat[ell])/(fabs(mu1_hat_pre[ell])+del );
diff[2+L+ell]=fabs(sigma1_sq_hat_pre[ell]-sigma1_sq_hat[ell])/(fabs(sigma1_sq_hat_pre[ell])+del );
diff[2+2*L+ell]=fabs(pEll_hat_pre[ell]-pEll_hat[ell])/(fabs(pEll_hat_pre[ell])+del );
fileout<<"mu1_hat["<<ell<<"]="<<mu1_hat[ell]<<endl;
fileout<<"sigma1_sq_hat["<<ell<<"]="<<sigma1_sq_hat[ell]<<endl;
fileout<<"pEll_hat["<<ell<<"]="<<pEll_hat[ell]<<endl;
}
}
d1=*max_element(diff.begin(), diff.end());
fileout<<"The scaled distance between previous and current estimates, d1="<<d1<<endl<<endl;
if(d1<1e-3 && alam==1 && stpover==0)
{
succ_t++;
}
else
{
succ_t=0;
}
if(succ_t==3)
{
n_tout=0;//satisfying the stopping rule for consecutive three times
}
}// t iteration is over
/////////////////////////////////////////////////////////
// ofstream fileout;
ofstream fileout2;
fileout2.open (filename);
fileout<<"Summary:"<<endl;
fileout<<"The number of voxels, pN="<<pN<<endl;
fileout<<"The iterations of GEM used="<<t<<endl;
fileout<<"The indicator of satisfying the stopping rule for consecutive three times="<<1-n_tout<<endl;
fileout<<"The scaled distance between previous and current estimates, d1="<<d1<<endl;
fileout<<"The proportion of the regular Newton step used in the final GEM iteration,lambda="<<alam<<endl;
fileout<<"The score function U:"<<endl;
fileout<<"U[0]="<<U[0]<<endl;
fileout<<"H_cond_mean[0]="<<H_cond_mean[0]<<endl;
fileout<<"H_mean[0]="<<H_mean[0]<<endl;
fileout<<"U[1]="<<U[1]<<endl;
fileout<<"H_cond_mean[1]="<<H_cond_mean[1]<<endl;
fileout<<"H_mean[1]="<<H_mean[1]<<endl;
fileout<<"The final estimates:"<<endl;
fileout<<"beta_hat="<<beta_hat<<endl;
fileout<<"h_hat="<<h_hat<<endl;
for(ell=0;ell<L;ell++)
{
fileout<<"mu1_hat["<<ell<<"]="<<mu1_hat[ell]<<endl;
fileout<<"sigma1_sq_hat["<<ell<<"]="<<sigma1_sq_hat[ell]<<endl;
fileout<<"pEll_hat["<<ell<<"]="<<pEll_hat[ell]<<endl;
}
fileout2<<beta_hat<<" "<<h_hat;
for(ell=0;ell<L;ell++)
{
fileout2<<" "<<mu1_hat[ell];
fileout2<<" "<<sigma1_sq_hat[ell];
fileout2<<" "<<pEll_hat[ell];
}
fileout2.close();
}
void braindt::gem(char* filename)//The Generalized EM Algorithm
{
int i,j,sum,sum_cond,sweep,ell,pi;
int t;///////
double TINY_L=TINY/L;
int error_time=0;
double temp_like;
int seed=0;
CRandomMersenne rnd1(seed);
StochasticLib1 rnd2(seed);
double a,b;//for penalized MLE (PMLE) when L>=2
a=1;
b=2;
double del=1e-3;
int succ_t=0;// to see the stopping rule successive time in the t iteration
double exp_x[7];
vector<double> exp_y(pN );
double beta_hat_pre,h_hat_pre;//for t iteration, i.e., the interations of GEM algorithm
double d1=1;//error in t iteration of GEM
vector<double> diff(2+3*L);
double p_1;
double temp2;
vector<double> constant(L );
vector<int> x(N );//x is the unobservable state 0 or 1
vector<int> x_cond(N );// conditional x on y
double constant_0=sqrt(2*M_PI*sigma0_sq);
vector<double> temp1(pN);
for(pi=0;pi<pN ;pi++)
{
initial[pi]=rnd1.Random()<0.5?1:0;
temp1[pi]=constant_0*exp(pow(y[pi]-mu0,2)/(2*sigma0_sq));
}
vector<double> log_q2_sub(swp_r); // -phi^T*H(x)
alam=0;
int n_tout=1;
for(t=0;n_tout && t<max_iter;t++)
{
beta_hat_pre=beta_hat;
h_hat_pre=h_hat;
vector<double> gamma_1(pN);
for(i=0;i<L;i++)
{
constant[i]=pEll_hat[i]/sqrt(2*M_PI*sigma1_sq_hat[i]);
}
for(pi=0;pi<pN;pi++)
{
temp2=0;
i=loca[pi];
x[i]=x_cond[i]=initial[pi];
for(j=0;j<L;j++)
{
temp2+=constant[j]*exp(-pow(y[pi]-mu1_hat[j],2)/(2*sigma1_sq_hat[j]));
}
exp_y[pi]=temp1[pi]*temp2;
}
for(i=0; i<7; i++)
{
exp_x[i]=exp(beta_hat*i+h_hat);
}
H_cond_mean=H_mean0;//set to be zeros
for(sweep=1;sweep<swp_b;sweep++)
{
for(pi=0; pi<pN;pi++)
{
i=loca[pi];
sum=x[i+1]+x[i-1]+x[i+Lx]+x[i-Lx]+x[i+YbyX]+x[i-YbyX];
sum_cond=x_cond[i+1]+x_cond[i-1]+x_cond[i+Lx]+x_cond[i-Lx]+x_cond[i+YbyX]+x_cond[i-YbyX];
p_1=1-1/(1+exp_x[sum]);
x[i]=rnd2.Bernoulli(p_1);
p_1=1-1/(1+exp_x[sum_cond]*exp_y[pi]);
x_cond[i]=rnd2.Bernoulli(p_1);
}
}
H_mean=H_mean0;//set to be zeros
H=H0;
for(sweep=0;sweep<swp_r;sweep++)
{
for(pi=0; pi<pN;pi++)
{
i=loca[pi];
sum=x[i+1]+x[i-1]+x[i+Lx]+x[i-Lx]+x[i+YbyX]+x[i-YbyX];
sum_cond=x_cond[i+1]+x_cond[i-1]+x_cond[i+Lx]+x_cond[i-Lx]+x_cond[i+YbyX]+x_cond[i-YbyX];
p_1=1-1/(1+exp_x[sum]);
x[i]=rnd2.Bernoulli(p_1);
p_1=1-1/(1+exp_x[sum_cond]*exp_y[pi]);
x_cond[i]=rnd2.Bernoulli(p_1);
}
for(pi=0; pi<pN;pi++)
{
i=loca[pi];
sum=x[i+1]+x[i-1]+x[i+Lx]+x[i-Lx]+x[i+YbyX]+x[i-YbyX];
sum_cond=x_cond[i+1]+x_cond[i-1]+x_cond[i+Lx]+x_cond[i-Lx]+x_cond[i+YbyX]+x_cond[i-YbyX];
H_cond_mean[0]+=x_cond[i]*sum_cond;//warning:1/2 in the later, for compute twice//Don't forgot to divide swp_r at the end
H_cond_mean[1]+=x_cond[i];
H[sweep][0]+=x[i]*sum*0.5;
H[sweep][1]+=x[i];
gamma_1[pi]+=x_cond[i];//r_[i][1]
//////////////////////////////////////////////////////////////
}
log_q2_sub[sweep]=-H[sweep][0]*beta_hat-H[sweep][1]*h_hat;
//////////////////////////////////////////////////////////////////
H_mean[0]+=H[sweep][0];//don't forgot to divide swp_r at the end
H_mean[1]+=H[sweep][1];
}
///////////////////////////////////////////////////////////////////////
temp_like=*max_element(log_q2_sub.begin(), log_q2_sub.end());
log_plikelihood=0;
for(sweep=0;sweep<swp_r;sweep++)
{
log_plikelihood+=exp(log_q2_sub[sweep]-temp_like);
}
q2=temp_like+log(log_plikelihood);//-logZ
///////////////// U
H_cond_mean[0]/=2;
for(i=0;i<2;i++)
{
H_cond_mean[i]/=swp_r;
H_mean[i]/=swp_r;
U[i]=H_cond_mean[i]-H_mean[i];
}
two_d_var_matrix_inv( swp_r);// for I update
//////////////////////////revised by 20140219
for(pi=0;pi<pN;pi++)
{
gamma_1[pi]=gamma_1[pi]/swp_r;
}
double gamma_1_sum=sigmasum(gamma_1,0,pN);
mu1_hat_pre=mu1_hat;
sigma1_sq_hat_pre=sigma1_sq_hat;
pEll_hat_pre=pEll_hat;
if(L==1)
{
pEll_hat[0]=1;
mu1_hat[0]=(sigmasum(gamma_1,y,0,pN)+TINY)/(gamma_1_sum+TINY);//+TINY to avoid the case gamma_1_sum=0
sigma1_sq_hat[0]=(sigmasum(gamma_1,y,mu1_hat[0],0,pN )+TINY)/(gamma_1_sum+TINY);
}
else
{
vector<double> f1y(pN);
vector<double> omega_sum(L);
vector<vector<double> > omega(L,vector<double>(pN));
for(pi=0;pi<pN;pi++)
{
for(ell=0;ell<L;ell++)
{
f1y[pi]+=pEll_hat[ell]*normpdf(y[pi], mu1_hat[ell], sigma1_sq_hat[ell]);
}
for(ell=0;ell<L;ell++)
{
omega[ell][pi]=gamma_1[pi]*pEll_hat[ell] \
*normpdf(y[pi],mu1_hat[ell],sigma1_sq_hat[ell])/f1y[pi];
omega_sum[ell]+=omega[ell][pi];
}
}
for(ell=0;ell<L;ell++)
{
mu1_hat[ell]=(sigmasum(omega[ell],y,0,pN )+TINY)/(omega_sum[ell]+TINY);// +TINY to avoid omega_sum[ell]=0
sigma1_sq_hat[ell]=(sigmasum(omega[ell],y,mu1_hat[ell],0,pN )+2*a)/(omega_sum[ell]+2*b);//PMLE
if(ell==(L-1))
{
pEll_hat[ell]=0;
for(int ell_sub=0;ell_sub<ell;ell_sub++)
{
pEll_hat[ell]+=pEll_hat[ell_sub];
}
pEll_hat[ell]=1-pEll_hat[ell];
///////if pELL_hat[L-1]<0, because of precision problem, to be like -2.22045e-16
if(pEll_hat[ell]<0)
{
pEll_hat[(ell-1)]+=pEll_hat[ell];
pEll_hat[ell]=0;
}
}
else {
pEll_hat[ell]=(omega_sum[ell]+TINY_L)/(gamma_1_sum+TINY);//+TINY to avoid gamma_1_sum=0
}
}
}
backtrack(rnd2);
//m iteration is over
///////////////////////
diff[0]=fabs(beta_hat_pre-beta_hat)/(fabs(beta_hat_pre)+del);
diff[1]=fabs(h_hat_pre-h_hat)/(fabs(h_hat_pre)+del);
for(ell=0;ell<L;ell++)
{
if(pEll_hat_pre[ell]<1e-3 && pEll_hat[ell]<1e-3)
{
diff[2+ell]=0;
diff[2+L+ell]=0;
diff[2+2*L+ell]=0;
}
else
{
diff[2+ell]=fabs(mu1_hat_pre[ell]-mu1_hat[ell])/(fabs(mu1_hat_pre[ell])+del );
diff[2+L+ell]=fabs(sigma1_sq_hat_pre[ell]-sigma1_sq_hat[ell])/(fabs(sigma1_sq_hat_pre[ell])+del );
diff[2+2*L+ell]=fabs(pEll_hat_pre[ell]-pEll_hat[ell])/(fabs(pEll_hat_pre[ell])+del );
}
}
d1=*max_element(diff.begin(), diff.end());
if(d1<1e-3 && alam==1 && stpover==0)
{
succ_t++;
}
else
{
succ_t=0;
}
if(succ_t==3)
{
n_tout=0;//satisfying the stopping rule for consecutive three times
}
}// t iteration is over
/////////////////////////////////////////////////////////
// ofstream fileout2;
ofstream fileout2;
fileout2.open (filename);
fileout2<<beta_hat<<" "<<h_hat;
for(ell=0;ell<L;ell++)
{
fileout2<<" "<<mu1_hat[ell];
fileout2<<" "<<sigma1_sq_hat[ell];
fileout2<<" "<<pEll_hat[ell];
}
fileout2.close();
}
void braindt::func(ofstream &fileout,StochasticLib1 &rnd2)//update U and q2 in the Armijo condition
{
vector<double> log_q2_sub(swp_r);// -phi^T*H(x)
double temp_like;
H_mean=H_mean0;
H=H0;
int pi,i,sweep,sum;
double p_1;
vector<int> x(N );
for(pi=0;pi<pN;pi++)
{
i=loca[pi];
x[i]=initial[pi];
}
double exp_x[7];
for(i=0; i<7; i++)
{
exp_x[i]=exp(beta_hat*i+h_hat);
}
for(sweep=1;sweep<swp_b;sweep++)
{
for(pi=0; pi<pN;pi++)
{
i=loca[pi];
sum=x[i+1]+x[i-1]+x[i+Lx]+x[i-Lx]+x[i+YbyX]+x[i-YbyX];
p_1=1-1/(1+exp_x[sum]);
x[i]=rnd2.Bernoulli(p_1);
}
}
for(sweep=0;sweep<swp_r;sweep++)
{
for(pi=0; pi<pN;pi++)
{
i=loca[pi];
sum=x[i+1]+x[i-1]+x[i+Lx]+x[i-Lx]+x[i+YbyX]+x[i-YbyX];
p_1=1-1/(1+exp_x[sum]);
x[i]=rnd2.Bernoulli(p_1);
}
for(pi=0;pi<pN;pi++)
{
i=loca[pi];
sum=x[i+1]+x[i-1]+x[i+Lx]+x[i-Lx]+x[i+YbyX]+x[i-YbyX];
H[sweep][0]+=x[i]*sum*0.5;
H[sweep][1]+=x[i];
//////////////////////////////////////////////////////////////
}
log_q2_sub[sweep]=-H[sweep][0]*beta_hat-H[sweep][1]*h_hat;
//////////////////////////////////////////////////////////////////
H_mean[0]+=H[sweep][0];//don't forgot to divide swp_r at the end
H_mean[1]+=H[sweep][1];
}
///////////////////////////////
temp_like=*max_element(log_q2_sub.begin(), log_q2_sub.end());
log_plikelihood=0;
for(sweep=0;sweep<swp_r;sweep++)
{
log_plikelihood+=exp(log_q2_sub[sweep]-temp_like);
}
q2=temp_like+log(log_plikelihood);//-logZ
//fileout<<"q2="<<q2<<endl;
fileout<<"The score function U:"<<endl;
for(i=0;i<2;i++)
{
H_mean[i]/=swp_r;
U[i]=H_cond_mean[i]-H_mean[i];
fileout<<"U["<<i<<"]="<<U[i]<<endl;
fileout<<"H_cond_mean["<<i<<"]="<<H_cond_mean[i]<<endl;
fileout<<"H_mean["<<i<<"]="<<H_mean[i]<<endl;
}
}
void braindt::func(StochasticLib1 &rnd2)//update U and q2 in the Armijo condition
{
vector<double> log_q2_sub(swp_r);// -phi^T*H(x)
double temp_like;
H_mean=H_mean0;
H=H0;
int pi,i,sweep,sum;
double p_1;
vector<int> x(N );
for(pi=0;pi<pN;pi++)
{
i=loca[pi];
x[i]=initial[pi];
}
double exp_x[7];
for(i=0; i<7; i++)
{
exp_x[i]=exp(beta_hat*i+h_hat);
}
for(sweep=1;sweep<swp_b;sweep++)
{
for(pi=0; pi<pN;pi++)
{
i=loca[pi];
sum=x[i+1]+x[i-1]+x[i+Lx]+x[i-Lx]+x[i+YbyX]+x[i-YbyX];
p_1=1-1/(1+exp_x[sum]);
x[i]=rnd2.Bernoulli(p_1);
}
}
for(sweep=0;sweep<swp_r;sweep++)
{
for(pi=0; pi<pN;pi++)
{
i=loca[pi];
sum=x[i+1]+x[i-1]+x[i+Lx]+x[i-Lx]+x[i+YbyX]+x[i-YbyX];
p_1=1-1/(1+exp_x[sum]);
x[i]=rnd2.Bernoulli(p_1);
}
for(pi=0;pi<pN;pi++)
{
i=loca[pi];
sum=x[i+1]+x[i-1]+x[i+Lx]+x[i-Lx]+x[i+YbyX]+x[i-YbyX];
H[sweep][0]+=x[i]*sum*0.5;
H[sweep][1]+=x[i];
//////////////////////////////////////////////////////////////
}
log_q2_sub[sweep]=-H[sweep][0]*beta_hat-H[sweep][1]*h_hat;
//////////////////////////////////////////////////////////////////
H_mean[0]+=H[sweep][0];//don't forgot to divide swp_r at the end
H_mean[1]+=H[sweep][1];
}
///////////////////////////////
temp_like=*max_element(log_q2_sub.begin(), log_q2_sub.end());
log_plikelihood=0;
for(sweep=0;sweep<swp_r;sweep++)
{
log_plikelihood+=exp(log_q2_sub[sweep]-temp_like);
}
q2=temp_like+log(log_plikelihood);//-logZ
}
void braindt::backtrack(ofstream &fileout,StochasticLib1 &rnd2)//The backtracking line search method
{//see P479 code of the book Numerical Recipe, 3rd
stpover=0;
q2_old=q2;
double del=1e-3;
beta_hat_old=beta_hat;
h_hat_old=h_hat;
const double ALF=1e-4;//alpha in the Armijo condition
const double tolx=1e-4;//tolerance for error in m iteration of backtracking
const double stpmax=5;// the maximum length of the newton step
double alamin;//the minimum lambda in the Armijo condition
double temp,test;
delta[0]=invI[0]*U[0]+invI[1]*U[1];
delta[1]=invI[1]*U[0]+invI[2]*U[1];
double sum=delta[0]*delta[0]+delta[1]*delta[1];//delta=p in the book, is the regular newton step
sum=sqrt(sum);
if(sum>stpmax)//Scale if attempted step is too big.
{
stpover=1;
delta[0]*=stpmax/sum;
delta[1]*=stpmax/sum;
//fileout<<"in backtrack firt delta1="<<delta[0]<<endl;
//fileout<<"in backtrack firt delta2="<<delta[1]<<endl;
fileout<<"Over the max step in backtracking,the step is thus scaled"<<endl;
}
double slope=U[0]*delta[0]+U[1]*delta[1];
//fileout<<"U*invI*U="<<slope<<endl;
temp=fabs(delta[0])/(fabs(beta_hat_old)+del);
test=fabs(delta[1])/(fabs(h_hat_old)+del);
test=test>temp?test:temp;
alamin=tolx/test;
alam=1;
int m=0;//
fileout<<"The backtracking iteration m=0"<<endl;
fileout<<"lambda="<<alam<<endl;
for(;;)
{
if(alam<alamin)
{
fileout<<"Converge at the case that lambda less than the minimum limit"<<endl;
beta_hat=beta_hat_old;
h_hat=h_hat_old;
break;
}
else {
beta_hat=beta_hat_old+alam*delta[0];
h_hat=h_hat_old+alam*delta[1];
fileout<<"beta_hat="<<beta_hat<<endl;
fileout<<"h_hat="<<h_hat<<endl;
func(fileout,rnd2);//update q2
// fileout<<"q2_old="<<q2_old<<endl;
// fileout<<"q2="<<q2<<endl;
delta_q2=(beta_hat-beta_hat_old)*H_cond_mean[0]+(h_hat-h_hat_old)*H_cond_mean[1]+q2-q2_old;//q2-q2old
fileout<<"delta_Q2="<<delta_q2<<endl;
if(delta_q2<(ALF*alam*slope))
{
m++;
fileout<<"The backtracking iteration m="<<m<<endl;
alam*=0.5;
fileout<<"lambda="<<alam<<endl;
}
else{
fileout<<"End backtracking for the Q2 increase"<<endl;
break;
}
}
}
}
void braindt::backtrack(StochasticLib1 &rnd2)//The backtracking line search method
{//see P479 code of the book Numerical Recipe, 3rd
stpover=0;
q2_old=q2;
double del=1e-3;
beta_hat_old=beta_hat;
h_hat_old=h_hat;
const double ALF=1e-4;//alpha in the Armijo condition
const double tolx=1e-4;//tolerance for error in m iteration of backtracking
const double stpmax=5;// the maximum length of the newton step
double alamin;//the minimum lambda in the Armijo condition
double temp,test;
delta[0]=invI[0]*U[0]+invI[1]*U[1];
delta[1]=invI[1]*U[0]+invI[2]*U[1];
double sum=delta[0]*delta[0]+delta[1]*delta[1];//delta=p in the book, is the regular newton step
sum=sqrt(sum);
if(sum>stpmax)//Scale if attempted step is too big.
{
stpover=1;
delta[0]*=stpmax/sum;
delta[1]*=stpmax/sum;
}
double slope=U[0]*delta[0]+U[1]*delta[1];
temp=fabs(delta[0])/(fabs(beta_hat_old)+del);
test=fabs(delta[1])/(fabs(h_hat_old)+del);
test=test>temp?test:temp;
alamin=tolx/test;
alam=1;
for(;;)
{
if(alam<alamin)
{
beta_hat=beta_hat_old;
h_hat=h_hat_old;
break;
}
else {
beta_hat=beta_hat_old+alam*delta[0];
h_hat=h_hat_old+alam*delta[1];
func(rnd2);//update q2
delta_q2=(beta_hat-beta_hat_old)*H_cond_mean[0]+(h_hat-h_hat_old)*H_cond_mean[1]+q2-q2_old;//q2-q2old
if(delta_q2<(ALF*alam*slope))
{
alam*=0.5;
}
else{
break;
}
}
}
}
////////////////computing LIS
void braindt::computeLIS(char* filename)
{
int seed=0;
CRandomMersenne rnd1(seed);
StochasticLib1 rnd2(seed);
int i,j,k,nn,sum,sum_cond,sweep,ell,pi;
int m;///////
double exp_x[7];
vector<double> exp_y(pN );
double p_1;
double temp,temp2;
vector<double> constant(L );
vector<int> x_cond(N );// conditional x
double constant_0=sqrt(2*M_PI*sigma0_sq);
vector<double> temp1(pN);
for(pi=0;pi<pN ;pi++)
{
temp1[pi]=constant_0*exp(pow(y[pi]-mu0,2)/(2*sigma0_sq));
}
for(i=0; i<7; i++)
{
exp_x[i]=exp(beta_hat*i+h_hat);
}
for(i=0;i<L;i++)
{
constant[i]=pEll_hat[i]/sqrt(2*M_PI*sigma1_sq_hat[i]);
}
for(pi=0;pi<pN;pi++)
{
i=loca[pi];
x_cond[i]=initial[pi];//
temp2=0;
for(j=0;j<L;j++)
{
temp2+=constant[j]*exp(-pow(y[pi]-mu1_hat[j],2)/(2*sigma1_sq_hat[j]));
}
exp_y[pi]=temp1[pi]*temp2;
}
for(int sweep=1;sweep<swp_b;sweep++)
{
for(pi=0; pi<pN;pi++)
{
i=loca[pi];
sum_cond=x_cond[i+1]+x_cond[i-1]+x_cond[i+Lx]+x_cond[i-Lx]+x_cond[i+YbyX]+x_cond[i-YbyX];
p_1=1-1/(1+exp_x[sum_cond]*exp_y[pi]);
x_cond[i]=rnd2.Bernoulli(p_1);
}
}
for(int sweep=0;sweep<swp_lis;sweep++)
{
for(pi=0; pi<pN;pi++)
{
i=loca[pi];
sum_cond=x_cond[i+1]+x_cond[i-1]+x_cond[i+Lx]+x_cond[i-Lx]+x_cond[i+YbyX]+x_cond[i-YbyX];
p_1=1-1/(1+exp_x[sum_cond]*exp_y[pi]);
x_cond[i]=rnd2.Bernoulli(p_1);
LIS[pi]+=x_cond[i];
}
}
ofstream file;
file.open (filename);
for(pi=0;pi<pN ;pi++)
{
LIS[pi]=1-LIS[pi]/swp_lis;
file<<LIS[pi]<<" ";
}
file.close();
}
Computing file changes ...
