#include <itpp/itbase.h>
#include <stat/libEF.h>

using namespace itpp;

//These lines are needed for use of cout and endl
using std::cout;
using std::endl;

void disp(const vec &tmu, const mat &tR,const mat &Smp){
	int N = Smp.cols();
	vec Emu = Smp*ones(N) /N ;
	mat Er = (Smp*Smp.transpose())/N - outer_product(Emu,Emu);
	cout << "True mu:" << tmu <<endl;
	cout << "Emp  mu:" << Emu <<endl;
	
	cout << "True R:" << tR <<endl;
	cout << "Emp  R:" << Er <<endl;
}

int main() {

	RNG_randomize();

	RV rv("1","{x }","2","0");
	int N = 10000; //number of samples
	vec mu0 = "1.5 1.7";
	mat V0("1.2 0.3; 0.3 5");
	ldmat R = ldmat(V0);
	
	cout << "====== ENorm ====== " <<endl;
	enorm<ldmat> eN(rv);
	eN.set_parameters(mu0,R);
	mat Smp = eN.sample(N);

	disp(mu0,R.to_mat(),Smp);

	cout << "====== MlNorm ====== " <<endl;
	mat I = eye(2);
	vec lik(N);
	mlnorm<ldmat> ML(rv,rv);
	ML.set_parameters(I,R);
	Smp = ML.samplecond(mu0,lik,N);
	
	disp(mu0,R.to_mat(),Smp);

	cout << "====== EGamma ====== " <<endl;	
	vec a = "100000,10000";
	vec b = a/10.0;
	egamma eG(rv);
	eG.set_parameters(a,b);
	
	Smp = eG.sample(N);

	vec g_mu = elem_div(a,b);
	vec g_var = elem_div(a,pow(b,2.0));
	disp(g_mu,diag(g_var),Smp);

	cout << "====== MGamma ====== " <<endl;	
	mgamma mG(rv,rv);
	double k = 10.0;
	mG.set_parameters(k);

	Smp=mG.samplecond(mu0,lik,N);
	disp(mu0,pow(mu0,2.0)/k,Smp);

	//Exit program:
	return 0;

}