#include <estim/arx.h>
using namespace bdm;

int main() {
	// Setup model : ARX for 1D Gaussian	
	//Test constructor
	mat V0 = 0.00001*eye(2); V0(0,0)= 0.1; //
	ARX Ar; Ar.set_statistics(1, V0, -1.0);
				
	mat mu(1,1);
	mat R(1,1);
	Ar._e()->mean_mat(mu,R);
	cout << "Prior moments: mu="<< mu << ", R=" << R <<endl; 
	
	int ndat = 200;
	vec smp=randn(ndat);
	//
	mat Smp=ones(2,ndat);
	Smp.set_row(0,smp);
	//
	Ar.bayesB(Smp);
	// Ar is now filled with estimates of N(0,1);
	cout << "Empirical moments: mu=" << sum(smp)/ndat << ", R=" << sum_sqr(smp)/ndat - pow(sum(smp)/ndat,2) << endl;
	Ar._e()->mean_mat(mu,R);
	cout << "Posterior moments: mu="<< mu << ", R=" << R <<endl; 
	
	//////// TEST prediction
	vec x=linspace(-3.0,3.0,100);
	double xstep = 6.0/100.0;
	mat X(1,100);
	mat X2(2,100);
	X.set_row(0,x);
	X2.set_row(0,x);
	
	mlstudent* Ap = Ar.predictor_student();
	vec Ap_x=Ap->evallogcond_m(X,vec_1(1.0));
	vec ll_x = Ar.logpred_m(X2);
	
	cout << "normalize : " << xstep*sum(exp(Ap_x)) << endl;
	cout << "normalize : " << xstep*sum(exp(ll_x)) << endl;
	
	it_file it("arx_elem_test.it");
	it << Name("Ap_x") << Ap_x;
	it << Name("ll_x") << ll_x;
}