#include "k.h"
/*  C implementation of Gibbs sampler  
	x (input) is a vector of chi statistics, n - length of x
    output - deltot a vector of prob obs is an outlier (length n) */
void prog(double *x,int n,double *deltot,long *idum,double *A,int *delta) {
double	xsumsq=0.,lambda,p,num,r,fun1,fun2,cc,dd,rsig,
		zeta,kappa,tau2,epsilon;
int i,j,
	nu=0,
	df=0,				/* rchisq degrees of freedom*/
	itmp=0,				
	ngib,
	burn;
	
/* hardcoded constants													*/
zeta=1000.;kappa=4.;tau2=1.;epsilon=.01;ngib=500;burn=50;

rsig=median(x,&n)/.675;
kappa=kappa*rsig;
tau2=tau2*rsig*rsig;
for (i=0;i<n;i++) xsumsq+=x[i]*x[i];
lambda=xsumsq/n;
dd=1./(1./tau2+1./lambda);
for (i=0;i<n;i++) { 
	cc=dd*(x[i]/lambda+kappa/tau2);
	A[i]=dd*rnorm(idum)+cc;
	fun1=(x[i]-A[i])/sqrt(lambda);
	fun2=x[i]/sqrt(lambda);
	p=dnorm(fun1)*epsilon/(dnorm(fun1)*epsilon+dnorm(fun2)*(1.-epsilon));
	delta[i]=(p>runif(idum));
}
/* gibbs update */
for (i=0;i<burn+ngib;i++) {
	/* calc num <- sum((x-delta*A)^2+nu*zeta  */
	/* and df <- sum(1-delta)+nu for rchisq		*/				 
	num=(double)nu*zeta;
	df=nu;
	for (j=0;j<n;j++) {
		num+=pow(x[j]-delta[j]*A[j],2); 
	    df+=(1-delta[j]);
	} 
	lambda=num/rchisq((long)df,idum);
	dd=1/(1/tau2+1/lambda); /* dd is used in loop below */
	r=0;
	for (j=0;j<n;j++) {
		fun1=(x[j]-A[j])/sqrt(lambda);
		fun2=x[j]/sqrt(lambda);
		p=dnorm(fun1)*epsilon/(dnorm(fun1)*epsilon+dnorm(fun2)*(1.-epsilon));
		delta[j]=(p>runif(idum));
		if (i>=burn) deltot[j]+=p;
		/* dd=1/(1/tau2+1/lambda) doesn't change so it's before the loop */
		cc=dd*(x[j]/lambda+kappa/tau2);
		A[j]=delta[j]*(sqrt(dd)*rnorm(idum)+cc)+
			(1-delta[j])*(sqrt(tau2)*rnorm(idum)+kappa);
		r+=delta[j];
	} 
}
for (j=0;j<n;j++) deltot[j]=deltot[j]/ngib;
return;
}		   
double dnorm(double x) {
  return(exp(-.5*x*x)/sqrt(2.*PI));
}