Index: /applications/robust/main.cpp
===================================================================
--- /applications/robust/main.cpp (revision 1331)
+++ /applications/robust/main.cpp (revision 1336)
@@ -134,5 +134,5 @@
 	vector<vector<string>> strings;
 
-	char* file_strings[3] = {"c:\\ar_cauchy_single.txt", "c:\\ar_normal_single.txt","c:\\ar_student_single.txt"};
+	char* file_strings[3] = {"c:\\ar_student_single.txt", "c:\\ar_student_single.txt","c:\\ar_cauchy_single.txt"};
 
 	for(int i = 0;i<3;i++)
@@ -162,5 +162,5 @@
 		vector<vec> conditions;
 		//emlig* emliga = new emlig(2);
-		RARX* my_rarx = new RARX(2,70);
+		RARX* my_rarx = new RARX(3,20,true);
 
 		for(int k = 1;k<170;k++)
@@ -189,28 +189,43 @@
 					
 				
-				if(k>5)
+				if(k>10)
 				{
 					//my_rarx->posterior->step_me(0);
 
-					my_rarx->posterior->sample_mat(1);
-
+					mat samples = my_rarx->posterior->sample_mat(50);
+					
+					
 					cout << "MaxLik coords:" << my_rarx->posterior->minimal_vertex->get_coordinates() << endl;
 
-					ofstream myfile;
-					char fstring[80];
-					strcpy(fstring,file_strings[j]);
-					strcat(fstring,"_res.txt");
-
-					myfile.open(fstring,ios::app);
-					myfile << my_rarx->posterior->minimal_vertex->get_coordinates()[0];
-					if(k!=strings[j].size()-1)
+					for(int s = 0;s<samples.rows();s++)
 					{
-						myfile << ",";
+						
+						double avg_parameter = samples.get_row(s)*ones(samples.cols())/samples.cols();						
+						
+						ofstream myfile;
+						char fstring[80];
+						strcpy(fstring,file_strings[j]);
+
+						char es[5];
+						strcat(fstring,itoa(s,es,10));
+
+						strcat(fstring,"_res.txt");
+						
+
+						myfile.open(fstring,ios::app);
+						
+						//myfile << my_rarx->posterior->minimal_vertex->get_coordinates()[0];
+						myfile << avg_parameter;
+						
+						if(k!=strings[j].size()-1)
+						{
+							myfile << ",";
+						}
+						else
+						{
+							myfile << endl;
+						}
+						myfile.close();
 					}
-					else
-					{
-						myfile << endl;
-					}
-					myfile.close();
 				}					
 			}	
Index: /applications/robust/robustlib.h
===================================================================
--- /applications/robust/robustlib.h (revision 1335)
+++ /applications/robust/robustlib.h (revision 1336)
@@ -17,5 +17,5 @@
 #include <algorithm>
 	
-//using namespace bdm;
+using namespace bdm;
 using namespace std;
 using namespace itpp;
@@ -1831,5 +1831,7 @@
 			while(sample_mat.cols()<n)
 			{
-				cout << "*************************************" << endl;
+				//// cout << "*************************************" << endl;
+
+				
 				
 				double rnumber = randu()*sum_a;
@@ -1852,5 +1854,5 @@
 				}				
 
-				cout << "Toprow/Count: " << toprow_count << "/" << ordered_toprows.size() << endl;
+				//// cout << "Toprow/Count: " << toprow_count << "/" << ordered_toprows.size() << endl;
 				// cout << &sampled_toprow << ";";
 
@@ -1870,8 +1872,8 @@
 				}
 
-				cout << "Simplex/Count: " << simplex_count << "/" << sampled_toprow->triangulation.size() << endl;
-				cout << "Simplex factor: " << (*s_ref)->probability << endl;
-				cout << "Toprow factor:  " << sampled_toprow->probability << endl;
-				cout << "Emlig factor:   " << normalization_factor << endl;
+				//// cout << "Simplex/Count: " << simplex_count << "/" << sampled_toprow->triangulation.size() << endl;
+				//// cout << "Simplex factor: " << (*s_ref)->probability << endl;
+				//// cout << "Toprow factor:  " << sampled_toprow->probability << endl;
+				//// cout << "Emlig factor:   " << normalization_factor << endl;
 				// cout << &(*tri_ref) << endl;
 
@@ -1893,8 +1895,12 @@
 							if(sum_g>rnumber)
 							{
-								itpp::Gamma_RNG* gamma = new itpp::Gamma_RNG(conditions.size()-number_of_parameters,1/(*g_ref).second);
-								sigma = 1/(*gamma)();
+								//itpp::Gamma_RNG* gamma = new itpp::Gamma_RNG(conditions.size()-number_of_parameters,1/(*g_ref).second);
+								//sigma = 1/(*gamma)();
 								
-								cout << "Sigma mean:   " << (*g_ref).second/(conditions.size()-number_of_parameters-1) << endl;								
+								GamRNG.setup(conditions.size()-number_of_parameters,(*g_ref).second);
+																
+								sigma = 1/GamRNG();
+
+								// cout << "Sigma mean:   " << (*g_ref).second/(conditions.size()-number_of_parameters-1) << endl;								
 								break;
 							}							
@@ -1936,6 +1942,6 @@
 				while(!have_sigma);
 
-				cout << "Sigma: " << sigma << endl;
-				cout << "Nr. of runs: " << number_of_runs << endl;
+				//// cout << "Sigma: " << sigma << endl;
+				//// cout << "Nr. of runs: " << number_of_runs << endl;
 
 				int dimension = (*s_ref)->vertices.size()-1;
@@ -1946,5 +1952,5 @@
 				vertex* base_vert = *(*s_ref)->vertices.begin();
 
-				cout << "Base vertex coords(should be close to est. param.): " << base_vert->get_coordinates() << endl;
+				//// cout << "Base vertex coords(should be close to est. param.): " << base_vert->get_coordinates() << endl;
 				
 				int row_count = 0;
@@ -1954,5 +1960,5 @@
 					vec current_coords = (*vert_ref)->get_coordinates();
 
-					cout << "Coords of vertex[" << row_count << "]: " << current_coords << endl; 
+					//// cout << "Coords of vertex[" << row_count << "]: " << current_coords << endl; 
 					
 					vec relative_coords = current_coords-base_vert->get_coordinates();				
@@ -1963,11 +1969,11 @@
 				}				
 				
-				cout << "Jacobian: " << jacobian << endl;
-
-				cout << "Gradient before trafo:" << gradient << endl;
+				//// cout << "Jacobian: " << jacobian << endl;
+
+				//// cout << "Gradient before trafo:" << gradient << endl;
 								
 				gradient = jacobian*gradient;	
 
-				cout << "Gradient after trafo:" << gradient << endl;
+				//// cout << "Gradient after trafo:" << gradient << endl;
 
 				// vec normal_gradient = gradient/sqrt(gradient*gradient);
@@ -2009,11 +2015,11 @@
 				extended_rotation.ins_col(0,zeros(extended_rotation.rows()));
 
-				cout << "Extended rotation: " << extended_rotation << endl;
+				//// cout << "Extended rotation: " << extended_rotation << endl;
 				
 				vec minima = itpp::min(extended_rotation,2);
 				vec maxima = itpp::max(extended_rotation,2);
 
-				cout << "Minima: " << minima << endl;
-				cout << "Maxima: " << maxima << endl;
+				//// cout << "Minima: " << minima << endl;
+				//// cout << "Maxima: " << maxima << endl;
 
 				vec sample_coordinates;		
@@ -2033,5 +2039,5 @@
 						vec new_gradient = rotation_matrix*gradient;
 						
-						cout << "New gradient(should have only first component nonzero):" << new_gradient << endl;
+						//// cout << "New gradient(should have only first component nonzero):" << new_gradient << endl;
 
 						// cout << "Max: " << maxima[0] << "  Min: " << minima[0] << "  Grad:" << new_gradient[0] << endl;
@@ -2049,9 +2055,9 @@
 				}
 
-				cout << "Sampled coordinates(gradient direction): " << sample_coordinates << endl;
+				//// cout << "Sampled coordinates(gradient direction): " << sample_coordinates << endl;
 
 				sample_coordinates = rotation_matrix.T()*sample_coordinates;
 
-				cout << "Sampled coordinates(backrotated direction):" << sample_coordinates << endl;
+				//// cout << "Sampled coordinates(backrotated direction):" << sample_coordinates << endl;
 
 				
@@ -2077,7 +2083,9 @@
 					sample_coordinates.ins(0,sigma);
 					
-					cout << "Sampled coordinates(parameter space):" << sample_coordinates << endl;
+					//// cout << "Sampled coordinates(parameter space):" << sample_coordinates << endl;
 
 					sample_mat.ins_col(0,sample_coordinates);
+
+					cout << sample_mat.cols() << ",";
 				}
 
@@ -2088,4 +2096,5 @@
 			}
 
+			cout << endl;
 			return sample_mat;
 		}
@@ -2405,6 +2414,5 @@
 {
 private:
-
-	
+	bool has_constant;
 
 	int window_size;	
@@ -2415,6 +2423,8 @@
 	emlig* posterior;
 
-	RARX(int number_of_parameters, const int window_size)//:BM()
-	{
+	RARX(int number_of_parameters, const int window_size, bool has_constant)//:BM()
+	{
+		this->has_constant = has_constant;
+		
 		posterior = new emlig(number_of_parameters);
 
@@ -2422,7 +2432,21 @@
 	};
 
-	void bayes(const itpp::vec &yt, const itpp::vec &cond = "")
-	{
-		conditions.push_back(yt);		
+	void bayes(itpp::vec yt)
+	{
+		if(has_constant)
+		{
+			int c_size = yt.size();
+			
+			yt.ins(c_size,1.0);
+		}
+		
+		if(yt.size() == posterior->number_of_parameters+1)
+		{
+			conditions.push_back(yt);		
+		}
+		else
+		{
+			throw new exception("Wrong condition size for bayesian data update!");
+		}
 
 		//posterior->step_me(0);
