Index: applications/robust/robustlib.h
===================================================================
--- applications/robust/robustlib.h (revision 1394)
+++ applications/robust/robustlib.h (revision 1395)
@@ -76,4 +76,6 @@
 
 	double probability;
+
+	double volume;
 
 	vector<multimap<double,double>> positive_gamma_parameters;
@@ -2237,64 +2239,45 @@
 	{
 		double sigma = 0;
-		double pg_sum;
-		double ng_sum;
-		do
-		{			
-			double rnumber = randu();
-			
-						
-			double sum_g = 0;
-			for(int i = 0;i<sampled_simplex->positive_gamma_parameters.size();i++)
-			{
-				for(multimap<double,double>::iterator g_ref = sampled_simplex->positive_gamma_parameters[i].begin();g_ref != sampled_simplex->positive_gamma_parameters[i].end();g_ref++)
-				{
-					sum_g += (*g_ref).first/sampled_simplex->positive_gamma_sum;
-
-								
-					if(sum_g>rnumber)
-					{
-						//itpp::Gamma_RNG* gamma = new itpp::Gamma_RNG(conditions.size()-number_of_parameters,1/(*g_ref).second);
-						//sigma = 1/(*gamma)();
-						
-						// TODO: Tady je chyba, musi zaviset i na parametru order, coz je tady i, ale nevim, jestli +i
-						// nebo -i
-						GamRNG.setup(conditions.size()-number_of_parameters+3,(*g_ref).second);
-																	
-						sigma = 1/GamRNG();
-
-						// cout << "Sigma mean:   " << (*g_ref).second/(conditions.size()-number_of_parameters-1) << endl;								
-						break;
-					}							
-				}
-
-				if(sigma!=0)
-				{
+		double pg_sum;					
+		double rnumber = randu();			
+					
+		double sum_g = 0;
+		for(int i = 0;i<sampled_simplex->positive_gamma_parameters.size();i++)
+		{
+			for(multimap<double,double>::iterator g_ref = sampled_simplex->positive_gamma_parameters[i].begin();g_ref != sampled_simplex->positive_gamma_parameters[i].end();g_ref++)
+			{
+				sum_g += (*g_ref).first/sampled_simplex->positive_gamma_sum;
+
+							
+				if(sum_g>rnumber)
+				{					
+					GamRNG.setup(condition_order-number_of_parameters-1+i,1/(*g_ref).second);
+																
+					sigma = 1/GamRNG();
+					// cout << "Sigma mean:   " << (*g_ref).second/(conditions.size()-number_of_parameters-1) << endl;								
+
 					break;
-				}
+				}							
 			}
 
-			rnumber = randu();
-
-			pg_sum = 0;
-			for(vector<multimap<double,double>>::iterator v_ref = sampled_simplex->positive_gamma_parameters.begin();v_ref!=sampled_simplex->positive_gamma_parameters.end();v_ref++)
-			{
-				for(multimap<double,double>::iterator pg_ref = (*v_ref).begin();pg_ref!=(*v_ref).end();pg_ref++)
-				{
-					pg_sum += exp((sampled_simplex->min_beta-(*pg_ref).second)/sigma)*pow((*pg_ref).second/sigma,(int)conditions.size())*(*pg_ref).second/fact(conditions.size())*(*pg_ref).first;
-				}					
+			if(sigma!=0)
+			{
+				break;
 			}
-
-			ng_sum = 0;
-			for(vector<multimap<double,double>>::iterator v_ref = sampled_simplex->negative_gamma_parameters.begin();v_ref!=sampled_simplex->negative_gamma_parameters.end();v_ref++)
-			{
-				for(multimap<double,double>::iterator ng_ref = (*v_ref).begin();ng_ref!=(*v_ref).end();ng_ref++)
-				{
-					ng_sum += exp((sampled_simplex->min_beta-(*ng_ref).second)/sigma)*pow((*ng_ref).second/sigma,(int)conditions.size())*(*ng_ref).second/fact(conditions.size())*(*ng_ref).first;
-				}					
+		}	
+
+		pg_sum = 0;
+		int i = 0;
+		for(vector<multimap<double,double>>::iterator v_ref = sampled_simplex->positive_gamma_parameters.begin();v_ref!=sampled_simplex->positive_gamma_parameters.end();v_ref++)
+		{
+			for(multimap<double,double>::iterator pg_ref = (*v_ref).begin();pg_ref!=(*v_ref).end();pg_ref++)
+			{
+				pg_sum += exp(-(*pg_ref).second/sigma)*pow((*pg_ref).second/sigma,condition_order-number_of_parameters-1+i)/sigma/fact(condition_order-number_of_parameters-1+i)*(*pg_ref).first;  // pg_sum += exp((sampled_simplex->min_beta-(*pg_ref).second)/sigma)*pow((*pg_ref).second/sigma,(int)conditions.size())*(*pg_ref).second/fact(conditions.size())*(*pg_ref).first;
 			}
-		}
-		while(pg_sum-ng_sum<0);
-
-		return pair<double,double>((pg_sum-ng_sum)/pg_sum,sigma);
+
+			i++;
+		}	
+
+		return pair<double,double>(sum_g/pg_sum,sigma);
 	}
 
@@ -2656,6 +2639,5 @@
 
 			double exponent = extended_coords*condition_and_simplex.first;
-			double sample_prob = 1/condition_and_simplex.second->probability/pow(probability_and_sigma.second,(int)conditions.size()-number_of_parameters+3)*exp((-1)/probability_and_sigma.second*exponent);
-			sample_prob *= probability_and_sigma.first;
+			double sample_prob = 1*condition_and_simplex.second->volume/condition_and_simplex.second->probability/pow(2*probability_and_sigma.second,condition_order)*exp(-exponent/probability_and_sigma.second)*probability_and_sigma.first;			
 
 			sample_coords.ins(sample_coords.size(),probability_and_sigma.second);
@@ -2663,4 +2645,7 @@
 			samples.ins_col(0,sample_coords);
 			probabilities.ins(0,sample_prob);
+
+			cout << "C:" << sample_coords << "   p:" << sample_prob << endl;
+			pause(1 );
 		}
 	
