Changeset 1325 for applications/robust

Timestamp:

04/12/11 18:11:50 (13 years ago)

Author:

sindj

Message:

Pokracuji v dodelavani samplovani, problem se zapornymi vahami pred gamma rozdelenimi sigmy. JS

Location:

applications/robust

Files:

• robustlib.cpp (modified) (4 diffs)
• robustlib.h (modified) (5 diffs)

applications/robust/robustlib.cpp

@@ -97 +97 @@
 
                         emlig* current_emlig;
+                        simplex->gamma_parameters.clear();
+                        simplex->gamma_sum = 0;
 
                         if(this->my_emlig!=NULL)
@@ -250 +252 @@
                                         }
                                         
-                                }       
-
-                                
-
-                                double bracket = fact(sigma_order)/fact(as1_order-1)/pow(a_0-(*as_ref).first,sigma_order+1);
+                                }
+                                
+                                                                
+
+                                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));
+
+                                bracket *= fact(sigma_order);
+
                                 for(int k = 0;k < as1_order-1;k++)
                                 {
-                                        double bracket_factor = -pow((double)-1,k+1)*fact(sigma_order-k)/fact(as1_order-1-k)/pow(a_0-(*as_ref).first,sigma_order-k);
+                                        double local_bracket = 0;
+                                        double bracket_factor = -pow((double)-1,k+1)/fact(as1_order-1-k)/pow(a_0-(*as_ref).first,sigma_order-k);
                                         
                                         ivec control_vec = ivec();
-                                        control_vec.ins(0,my_emlig->number_of_parameters-as1_order+1);          
-                                        
+                                        control_vec.ins(0,my_emlig->number_of_parameters-as1_order+1);                                  
                                         
                                         for(multiset<my_ivec>::iterator combi_ref = this->my_emlig->correction_factors[k].begin();combi_ref!=this->my_emlig->correction_factors[k].upper_bound(my_ivec(control_vec));combi_ref++)
@@ -271 +282 @@
                                                         
                                                         bracket_combination /= factors[(*combi_ref)[j]-1];
-                                                }
-
-                                                bracket+=bracket_factor*bracket_combination;                                                                    
-                                        }                                                                       
-                                }
-
-                                
+                                                }                                               
+                                                                                                
+                                                local_bracket += bracket_factor*bracket_combination;                                                                    
+                                        }
+
+                                        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));
+
+                                        bracket += local_bracket*fact(sigma_order-k);
+                                }                               
 
                                 int_value += multiplier*bracket;
@@ -283 +304 @@
                         }
 
-                        
-                        correction_term *= fact(sigma_order)/abs(pow(a_0,sigma_order+1));
-
-                        // cout << c << int_value << endl;
+                        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));
+
+                        correction_term = fact(sigma_order)*correction_term_base;
+
+                        //cout << c << int_value << endl;
 
                         int_value += correction_term;
                 
 
-                        // cout << "Probability:" << int_value << endl;
-                        // pause(0.100);
+                        //cout << "Probability:" << int_value << endl;
+                        //pause(0.100);
 
                         simplex->probability = int_value;

applications/robust/robustlib.h

@@ -8 +8 @@
 #define ROBUST_H
 
-//#include <stat/exp_family.h>
+#include <stat/exp_family.h>
 #include <itpp/itbase.h>
 #include <map>
@@ -73 +73 @@
         double probability;
 
-        vector<list<pair<double,double>>> gamma_parameters;
+        vector<multimap<double,double>> gamma_parameters;
+
+        double gamma_sum;
+        
 
         simplex(set<vertex*> vertices)
@@ -1749 +1752 @@
                         map<double,toprow*> ordered_toprows;                    
                         double sum_a = 0;
-                        
 
                         for(polyhedron* top_ref = statistic.rows[number_of_parameters];top_ref!=statistic.end_poly;top_ref=top_ref->next_poly)
@@ -1784 +1786 @@
                                 set<simplex*>::iterator s_ref;
                                 double sum_b = 0;
-
                                 for(s_ref = sampled_toprow->triangulation.begin();s_ref!=sampled_toprow->triangulation.end();s_ref++)
                                 {
@@ -1794 +1795 @@
 
                                 // cout << &(*tri_ref) << endl;
+
+                                rnumber = randu();
+
+                                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)();
+                                                        break;
+                                                }                                               
+                                        }
+
+                                        if(sigma!=0)
+                                        {
+                                                break;
+                                        }
+                                }
 
                                 int dimension = (*s_ref)->vertices.size()-1;