Changeset 1331 for applications/robust

Timestamp:

04/14/11 19:24:43 (14 years ago)

Author:

sindj

Message:

Dodelan sampling sigma, zbyva doladit sampling alpha.JS

Location:

applications/robust

Files:

• main.cpp (modified) (1 diff)
• robustlib.cpp (modified) (4 diffs)
• robustlib.h (modified) (5 diffs)

applications/robust/main.cpp

@@ -193 +193 @@
                                         //my_rarx->posterior->step_me(0);
 
-                                        //my_rarx->posterior->sample_mat(1);
+                                        my_rarx->posterior->sample_mat(1);
 
                                         cout << "MaxLik coords:" << my_rarx->posterior->minimal_vertex->get_coordinates() << endl;

applications/robust/robustlib.cpp

@@ -97 +97 @@
 
                         emlig* current_emlig;
-                        simplex->gamma_parameters.clear();
-                        simplex->gamma_sum = 0;
+                        simplex->clear_gammas();
 
                         if(this->my_emlig!=NULL)
@@ -257 +256 @@
 
                                 double bracket = fact(as1_order-1)/pow(a_0-(*as_ref).first,sigma_order+1);
-                                
-                                simplex->gamma_sum += bracket*multiplier;
-                                
-                                multimap<double,double> map;
-                                simplex->gamma_parameters.push_back(map);
-                                simplex->gamma_parameters[0].insert(pair<double,double>(bracket*multiplier,a_0-(*as_ref).first));
+                                                                
+                                simplex->insert_gamma(0,bracket*multiplier,a_0-(*as_ref).first);                                                                
 
                                 bracket *= fact(sigma_order);
@@ -287 +282 @@
                                         }
 
-                                        simplex->gamma_sum += local_bracket*multiplier;
-                                        
-                                        if(simplex->gamma_parameters.size()<=k+1)
-                                        {
-                                                multimap<double,double> loc_map;
-                                                simplex->gamma_parameters.push_back(loc_map);
-                                        }
-
-                                        simplex->gamma_parameters[k+1].insert(pair<double,double>(local_bracket*multiplier,a_0-(*as_ref).first));
-
+                                        simplex->insert_gamma(k+1,local_bracket*multiplier,a_0-(*as_ref).first);
+                                        
                                         bracket += local_bracket*fact(sigma_order-k);
                                 }                               
@@ -306 +293 @@
                         double correction_term_base = correction_term/pow(a_0,sigma_order+1);
 
-                        simplex->gamma_sum += correction_term_base;
-                        simplex->gamma_parameters[0].insert(pair<double,double>(correction_term_base,a_0));
+                        simplex->insert_gamma(0,correction_term_base,a_0);                      
 
                         correction_term = fact(sigma_order)*correction_term_base;

applications/robust/robustlib.h

@@ -67 +67 @@
 class simplex
 {
+        
+
 public:
 
@@ -73 +75 @@
         double probability;
 
-        vector<multimap<double,double>> gamma_parameters;
-
-        double gamma_sum;
+        vector<multimap<double,double>> positive_gamma_parameters;
+
+        vector<multimap<double,double>> negative_gamma_parameters;
+
+        double positive_gamma_sum;
+
+        double negative_gamma_sum;
         
 
@@ -88 +94 @@
                 this->vertices.insert(vertex);
                 probability = 0;
+        }
+
+        void clear_gammas()
+        {
+                positive_gamma_parameters.clear();
+                negative_gamma_parameters.clear();
+                
+                
+                positive_gamma_sum = 0;
+                negative_gamma_sum = 0;
+        }
+
+        void insert_gamma(int order, double weight, double beta)
+        {
+                if(weight>=0)
+                {
+                        if(positive_gamma_parameters.size()<order+1)
+                        {
+                                multimap<double,double> map;
+                                positive_gamma_parameters.push_back(map);
+                        }
+
+                        positive_gamma_sum += weight;
+
+                        positive_gamma_parameters[order].insert(pair<double,double>(weight,beta));              
+                }
+                else
+                {
+                        if(negative_gamma_parameters.size()<order+1)
+                        {
+                                multimap<double,double> map;
+                                negative_gamma_parameters.push_back(map);
+                        }
+
+                        negative_gamma_sum -= weight;
+
+                        negative_gamma_parameters[order].insert(pair<double,double>(-weight,beta));
+                }
         }
 };
@@ -1742 +1786 @@
 
         mat sample_mat(int n)
-        {
-                
-                
+        {               
 
                 /// \TODO tady je to spatne, tady nesmi byt conditions.size(), viz RARX.bayes()
@@ -1796 +1838 @@
                                 // cout << &(*tri_ref) << endl;
 
-                                rnumber = randu();
-
+                                bool have_sigma = false;
                                 double sigma = 0;
-                                double sum_g = 0;
-                                for(int i = 0;i<(*s_ref)->gamma_parameters.size();i++)
-                                {
-                                        for(multimap<double,double>::iterator g_ref = (*s_ref)->gamma_parameters[i].begin();g_ref != (*s_ref)->gamma_parameters[i].end();g_ref++)
-                                        {
-                                                sum_g += (*g_ref).first/(*s_ref)->gamma_sum;
-
-                                                if(sum_g>rnumber)
-                                                {
-                                                        itpp::Gamma_RNG* gamma = new itpp::Gamma_RNG(conditions.size()-number_of_parameters,(*g_ref).second);
-                                                        sigma = (*gamma)();
+                                do
+                                {
+                                        rnumber = randu();
+                                        
+                                        double sum_g = 0;
+                                        for(int i = 0;i<(*s_ref)->positive_gamma_parameters.size();i++)
+                                        {
+                                                for(multimap<double,double>::iterator g_ref = (*s_ref)->positive_gamma_parameters[i].begin();g_ref != (*s_ref)->positive_gamma_parameters[i].end();g_ref++)
+                                                {
+                                                        sum_g += (*g_ref).first/(*s_ref)->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)();
+                                                                break;
+                                                        }                                               
+                                                }
+
+                                                if(sigma!=0)
+                                                {
                                                         break;
-                                                }                                               
-                                        }
-
-                                        if(sigma!=0)
-                                        {
-                                                break;
-                                        }
-                                }
+                                                }
+                                        }
+
+                                        rnumber = randu();
+
+                                        double pg_sum = 0;
+                                        for(vector<multimap<double,double>>::iterator v_ref = (*s_ref)->positive_gamma_parameters.begin();v_ref!=(*s_ref)->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,(int)conditions.size()-number_of_parameters-1)*(*pg_ref).second/fact(conditions.size()-number_of_parameters-1)*(*pg_ref).first;
+                                                }                                       
+                                        }
+
+                                        double ng_sum = 0;
+                                        for(vector<multimap<double,double>>::iterator v_ref = (*s_ref)->negative_gamma_parameters.begin();v_ref!=(*s_ref)->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(-(*ng_ref).second/sigma)*pow((*ng_ref).second/sigma,(int)conditions.size()-number_of_parameters-1)*(*ng_ref).second/fact(conditions.size()-number_of_parameters-1)*(*ng_ref).first;
+                                                }                                       
+                                        }
+                                        
+                                        if((pg_sum-ng_sum)/pg_sum>rnumber)
+                                        {
+                                                have_sigma = true;
+                                        }
+                                }
+                                while(!have_sigma);
 
                                 int dimension = (*s_ref)->vertices.size()-1;