userddi.c

he AutoMix package is a C program for Unix-like systems, implementing the automatic reversible jum

	  free(Umin1[j]);
	}
      }
      free(C);
      free(R1min1);
      free(tempmat);
      free(tempmat2);
      free(U);
      free(Umin1);
      lp=-10000000.0;
      *llh1=lp;
      return lp;
    }
  }

  /* work out V or U (will be used in likelihood) */ 

  if(k==0){

    /* tempmat is matrix sqrt of V^(-1). tempmat2 will be inverse of this
       so that ((tempmat2)^T)(tempmat2)=V 
       and (tempmat)((tempmat)^T)=V^(-1) (same with U below) */ 

    for(j=0;j<3;j++){
      for(j1=0;j1<3;j1++){
	tempmat2[j][j1]=0.0;
      }
      tempmat2[j][j]=1.0;
    }
    for(j=0;j<3;j++){
      for(j1=0;j1<=j;j1++){
	for(j2=0;j2<j;j2++){
	  tempmat2[j][j1]-=tempmat[j][j2]*tempmat2[j2][j1];
	}
	tempmat2[j][j1]/=tempmat[j][j];
      }
    }
    for(j=0;j<3;j++){
      for(j1=0;j1<3;j1++){
        V[j][j1]=0.0;
        for(j2=0;j2<3;j2++){
	  V[j][j1]+=tempmat2[j2][j]*tempmat2[j2][j1];
	}
      }
    }
  }
  else{
    for(j=0;j<2;j++){
      for(j1=0;j1<2;j1++){
	tempmat2[j][j1]=0.0;
      }
      tempmat2[j][j]=1.0;
    }
    for(j=0;j<2;j++){
      for(j1=0;j1<=j;j1++){
	for(j2=0;j2<j;j2++){
	  tempmat2[j][j1]-=tempmat[j][j2]*tempmat2[j2][j1];
	}
	tempmat2[j][j1]/=tempmat[j][j];
      }
    }
    for(j=0;j<2;j++){
      for(j1=0;j1<2;j1++){
        U[j][j1]=0.0;
        for(j2=0;j2<2;j2++){
	  U[j][j1]+=tempmat2[j2][j]*tempmat2[j2][j1];
	}
      }
    }
  }



  /* Start with prior contribution */
  
  if(k==0){
    lp=1.0; 
    /* Normal for alpha */
    for(j=0;j<9;j++){
      lp*=D0min1[j];
    }
    lp=0.5*log(lp);
    lp-=4.5*log(tpi); 
    for(j=0;j<9;j++){
      lp-=0.5*(alpha[j]-c0[j])*(alpha[j]-c0[j])*D0min1[j];      
    }

    /* Wishart for V^(-1) */

    detVmin1=pow(det2(3,tempmat),2.0);

    lp+=((rho-3.0-1.0)/2.0)*log(detVmin1);
    for(j=0;j<3;j++){
      lp-=0.5*rho*R0[j][j]*Vmin1[j][j];
    }
    lp-=(rho/2.0)*log(pow(rho,-3.0)*det2(3,R0min1));
    lp-=(rho*3.0/2.0)*log(2.0);
    lp-=(3.0*2.0/4.0)*log(pi); 
    for(j=0;j<3;j++){
      lp-=loggamma((double)(rho-j)/2.0);
    }
    /* Inverse gamma for sigmasq */
    lp+=(-(a+1.0)*log(sigmasq)-a*log(b)-1.0/(b*sigmasq)-loggamma(a));
  }
  else{
    lp=1.0; 
    /* Normal for gamma */
    for(j=0;j<6;j++){
      lp*=D1min1[j];
    }
    lp=0.5*log(lp);
    lp-=3.0*log(tpi); 
    for(j=0;j<6;j++){
      lp-=0.5*(gamma[j]-c1[j])*(gamma[j]-c1[j])*D1min1[j];      
    }

    /* Wishart for U^(-1) */
    
    detUmin1=pow(det2(2,tempmat),2.0);
  
    lp+=((rho-2.0-1.0)/2.0)*log(detUmin1);
    for(j=0;j<2;j++){
      lp-=0.5*rho*R1[j][j]*Umin1[j][j];
    }
    lp-=(rho/2.0)*log(pow(rho,-2.0)*det2(2,R1min1));
    lp-=(rho*2.0/2.0)*log(2.0);
    lp-=(2.0*1.0/4.0)*log(pi); 
    for(j=0;j<2;j++){
      lp-=loggamma((double)(rho-j)/2.0);
    }
    /* Inverse gamma for tausq */
    lp+=(-(a+1.0)*log(tausq)-a*log(b)-1.0/(b*tausq)-loggamma(a));
  }


  /* Look at likelihood contribution */
  llh=0.0;
  if(k==0){
    for(i=0;i<467;i++){
      for(j=0;j<3;j++){
        for(j1=0;j1<5;j1++){
	  temp[j][j1]=0.0;
	}
	for(j1=0;j1<S[i];j1++){
	  for(j2=0;j2<3;j2++){
	    temp[j][j1]+=V[j][j2]*W[i][j1][j2];
	  }
	}
      }
      for(j=0;j<5;j++){
	for(j1=0;j1<5;j1++){
	  C[j][j1]=0.0;
	}
      }
      for(j=0;j<S[i];j++){
	for(j1=0;j1<S[i];j1++){
	  for(j2=0;j2<3;j2++){
	    C[j][j1]+=W[i][j][j2]*temp[j2][j1];
	  }
	}
      }
      for(j=0;j<S[i];j++){
	C[j][j]+=sigmasq;
      }
      posdef=1;
      chol2(S[i],C,&posdef);
      if(posdef==0){
	for(j=0;j<5;j++){
	  free(C[j]);
	  if(j<3){
	    free(R0min1[j]);
	    free(tempmat[j]);
	    free(tempmat2[j]);
	    free(V[j]);
	    free(Vmin1[j]);
	  }
	}
	free(C);
	free(R0min1);
	free(tempmat);
	free(tempmat2);
	free(V);
	free(Vmin1);	
	lp=-10000000.0;
	*llh1=lp;
	return lp;
      }
  
      for(j=0;j<S[i];j++){
	mu[j]=0;
        for(j1=0;j1<9;j1++){
	  mu[j]+=X[i][j][j1]*alpha[j1];
	}
      }

      j1=0;
      for(j=0;j<5;j++){
	if(Y[i][j]<90){
	  datai[j1]=Y[i][j];
	  j1++;
	}
      }

      llh+=lnormprob2(S[i],C,mu,datai);     
    }
  }
  else{
    for(i=0;i<467;i++){
      for(j=0;j<2;j++){
        for(j1=0;j1<5;j1++){
	  temp[j][j1]=0.0;
	}
	for(j1=0;j1<S[i];j1++){
	  for(j2=0;j2<2;j2++){
	    temp[j][j1]+=U[j][j2]*Q[i][j1][j2];
	  }
	}
      }
      for(j=0;j<5;j++){
	for(j1=0;j1<5;j1++){
	  C[j][j1]=0.0;
	}
      }
      for(j=0;j<S[i];j++){
	for(j1=0;j1<S[i];j1++){
	  for(j2=0;j2<2;j2++){
	    C[j][j1]+=Q[i][j][j2]*temp[j2][j1];
	  }
	}
      }
      for(j=0;j<S[i];j++){
	C[j][j]+=tausq;
      }
      posdef=1;
      chol2(S[i],C,&posdef);
      if(posdef==0){
	for(j=0;j<5;j++){
	  free(C[j]);
	  if(j<2){
	    free(R1min1[j]);
	    free(tempmat[j]);
	    free(tempmat2[j]);
	    free(U[j]);
	    free(Umin1[j]);
	  }
	}
	free(C);
	free(R1min1);
	free(tempmat);
	free(tempmat2);
	free(U);
	free(Umin1);
	lp=-10000000.0;
	*llh1=lp;
	return lp;
      }

      for(j=0;j<S[i];j++){
	mu[j]=0;
        for(j1=0;j1<6;j1++){
	  mu[j]+=P[i][j][j1]*gamma[j1];
	}
      }
      j1=0;
      for(j=0;j<5;j++){
	if(Y[i][j]<90){
	  datai[j1]=Y[i][j];
	  j1++;
	}
      }

      llh+=lnormprob2(S[i],C,mu,datai);     
    }
  }
  *llh1=llh;
  lp+=llh;

  if(k==0){
    for(j=0;j<5;j++){
      free(C[j]);
      if(j<3){
	free(R0min1[j]);
	free(tempmat[j]);
	free(tempmat2[j]);
	free(V[j]);
        free(Vmin1[j]);
      }
    }
    free(C);
    free(R0min1);
    free(tempmat);
    free(tempmat2);
    free(V);
    free(Vmin1);
  }
  else{
    for(j=0;j<5;j++){
      free(C[j]);
      if(j<2){
	free(R1min1[j]);
	free(tempmat[j]);
	free(tempmat2[j]);
	free(U[j]);
	free(Umin1[j]);
      }
    }
    free(C);
    free(R1min1);
    free(tempmat);
    free(tempmat2);
    free(U);
    free(Umin1);
  } 
	  
  return lp;
      
} 

double boxm2(){

  /*Function for returning single N(0,1) variable */

  double u1,u2,out;

  u1=sdrand();
  u2=sdrand();
  u1=tpi*u1;
  u2=sqrt(-2.0*log(u2));
  out=u2*cos(u1);

  return out;

}

void chol2(int nkk, double **A,int *posdef){

  /* Function for performing cholesky decompositon of B and returns result 
   in the same matrix - adapted from PJG Fortran function*/  
    
  int j1,j2,j3;
  double sum;

  for(j1=0;j1<nkk;j1++){
    sum=A[j1][j1];
    for(j2=0;j2<j1;j2++){
      sum-=pow(A[j1][j2],2);
    }
    if(sum<0){
      *(posdef)=0;
	return;
    }
    A[j1][j1]=sqrt(sum);
  
    for(j2=j1+1;j2<nkk;j2++){
      sum=A[j2][j1];
      for(j3=0;j3<j1;j3++){
	sum-=A[j2][j3]*A[j1][j3];
      }
      A[j2][j1]=sum/A[j1][j1];
    }
  }
}

double det2(int nkk,double **B){

  /* Function for evaluating determinant of an nkk by nkk matrix B */

  int j1;
  double out;
  out=1.0;
  for(j1=0;j1<nkk;j1++){
    out*=B[j1][j1];
  }
  return out;
}

double lnormprob2(int nkk, double **B, double *mu, double *datai){

  /* Function for evaluating the log of p.d.f. of Multivariate normal
   with mean mu and sqrt of cov matrix B, at pt datai */

  int j1,j2;
  double work[nkk];
  double out;
  
  for(j1=0;j1<nkk;j1++){
    work[j1]=datai[j1]-mu[j1];
  }
  for(j1=0;j1<nkk;j1++){
    for(j2=0;j2<j1;j2++){
      (work[j1])-=B[j1][j2]*work[j2];
    }
    (work[j1])/=B[j1][j1]; 
  }
  out=0.0;
  for(j1=0;j1<nkk;j1++){
    out+=(work[j1]*work[j1]);
  }
  out=-0.5*out;
  out-=(nkk/2.0)*log(tpi);
  out-=log(det2(nkk,B));    
  return out;

}