Changeset 1395 for applications/robust/robustlib.h

Timestamp:

09/21/11 20:35:02 (13 years ago)

Author:

sindj

Message:

Dodelavani toho zatracenyho samplingu. JS

Files:

• applications/robust/robustlib.h (modified) (4 diffs)

applications/robust/robustlib.h

@@ -76 +76 @@
 
         double probability;
+
+        double volume;
 
         vector<multimap<double,double>> positive_gamma_parameters;
@@ -2237 +2239 @@
         {
                 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 +2639 @@
 
                         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 +2645 @@
                         samples.ins_col(0,sample_coords);
                         probabilities.ins(0,sample_prob);
+
+                        cout << "C:" << sample_coords << "   p:" << sample_prob << endl;
+                        pause(1 );
                 }