Changeset 1272 for applications/robust

Timestamp:

01/14/11 18:28:55 (13 years ago)

Author:

sindj

Message:

Dodelano pocitani, odstranena chyba se zapornymi pravdepodobnostmi. JS

Location:

applications/robust

Files:

• main.cpp (modified) (4 diffs)
• robustlib.cpp (modified) (8 diffs)
• robustlib.h (modified) (3 diffs)

applications/robust/main.cpp

@@ -14 +14 @@
 
 const int emlig_size = 4;
+
 
 int main ( int argc, char* argv[] ) {
@@ -51 +52 @@
                 for(int k = 0;k<=emlig_size;k++)
                 {
-                        condition[k] = rand()/1000.0;           
+                        condition[k] = (rand()-RAND_MAX/2)/1000.0;              
                 }
                         
@@ -58 +59 @@
                 emlig1->add_condition(*condition_vec);
 
+                /*
+                for(polyhedron* toprow_ref = emlig1->statistic.rows[emlig_size]; toprow_ref != emlig1->statistic.end_poly; toprow_ref = toprow_ref->next_poly)
+                {
+                        cout << ((toprow*)toprow_ref)->probability << endl;
+                }
+                */
+
                 cout << emlig1->statistic_rowsize(emlig_size) << endl << endl;
         
+                /*
                 if(i-emlig1->number_of_parameters >= 0)
                 {
                         pause(30);
                 }
+                */
 
                 // emlig1->step_me(i);
                 
+                /*
                 vector<int> sizevector;
                 for(int s = 0;s<=emlig1->number_of_parameters;s++)
@@ -72 +83 @@
                         sizevector.push_back(emlig1->statistic_rowsize(s));
                 }
+                */
         }
 

applications/robust/robustlib.cpp

@@ -55 +55 @@
                         vec cur_condition = as_toprow->condition.get(1,as_toprow->condition.size()-1);
 
+                        // cout << as_toprow->condition << endl;
+
                         vertex* base_vertex = (*simplex.begin());
 
@@ -61 +63 @@
                         mat jacobian(dimension,dimension);
 
-                        double a_0 = base_vertex->get_coordinates()*cur_condition+as_toprow->condition[0];
+                        // cout << base_vertex->get_coordinates() << endl;
+
+                        double a_0 = base_vertex->get_coordinates()*cur_condition-as_toprow->condition[0];
                         map<double,int> as;
 
@@ -69 +73 @@
                                 vec relative_coords = (*vert_ref)->get_coordinates()-base_vertex->get_coordinates();
 
+                                // cout << (*vert_ref)->get_coordinates() << endl;
                                 // cout << relative_coords << endl;
 
@@ -80 +85 @@
                                         (*returned.first).second++;
                                 }
+                                /*
                                 else
                                 {
                                         cout << "a[" << row_count << "] = " << new_a << endl;
                                 }
+                                */
+                                
                                 //as.ins(as.size(),new_a);                                                      
                                 
@@ -89 +97 @@
                         }
 
-                        // cout << endl;
+                        
 
                         double int_value = 0;
@@ -95 +103 @@
                         // cout << jacobian << endl;
 
-                        double det_jacobian = det(jacobian);
+                        double det_jacobian    = abs(det(jacobian));
+                        double correction_term = det_jacobian;                  
                         for(map<double,int>::iterator as_ref = as.begin();as_ref!=as.end();as_ref++)
                         {
                                 double multiplier = det_jacobian;
+                                
                                 if(a_0!=(*as_ref).first)
-                                {                                                               
+                                {
                                         int as1_order = (*as_ref).second;
-
+                                        
+                                        correction_term /= -pow((*as_ref).first,as1_order);                                     
+                                        
                                         current_emlig->set_correction_factors(as1_order);
 
@@ -126 +138 @@
                                                         multiplier        /= (*as_ref).first;
                                                         number_of_factors += 1;
-                                                }                                                                       
+                                                }
+                                                
                                         }                                       
 
-                                        double bracket = fact(condition_order-number_of_factors)/fact(as1_order-1)/pow((*as_ref).first-a_0,condition_order-number_of_factors);
+                                        double bracket = fact(condition_order-number_of_factors)/fact(as1_order-1)/pow(a_0+(*as_ref).first,condition_order-number_of_factors+1);
                                         for(int k = 0;k < as1_order-1;k++)
                                         {
-                                                double bracket_factor = pow((double)-1,k+1)*fact(condition_order-1-number_of_factors-k)/fact(as1_order-2-k)/pow((*as_ref).first-a_0,condition_order-1-number_of_factors-k);
+                                                double bracket_factor = pow((double)-1,k+1)*fact(condition_order-1-number_of_factors-k)/fact(as1_order-2-k)/pow(a_0+(*as_ref).first,condition_order-number_of_factors-k);
                                                 
                                                 for(set<my_ivec>::iterator combi_ref = this->my_emlig->correction_factors[k].begin();combi_ref!=this->my_emlig->correction_factors[k].end();combi_ref++)
@@ -144 +157 @@
                                                         bracket+=bracket_factor*bracket_combination;                                                                    
                                                 }                                                                       
-                                        }                                                               
+                                        }
 
-                                        int_value += multiplier*bracket;        
+                                        
 
+                                        int_value += multiplier*bracket;
                                 }
                                 else
                                 {
                                         throw new exception("Should this happen? a_i=a_0 in the formula for integrating over a simplex! I think it can happen with 0 probability(when phi*2*v_0=phi*v_i)");
-                                }
+                                }                                                                                                               
+                        }
 
-                                // cout << c << int_value << endl;
+                        
+                        correction_term *= fact(condition_order-my_emlig->number_of_parameters)/pow(a_0,condition_order-my_emlig->number_of_parameters+1);
 
-                                return int_value;                                                                               
-                        }
+                        // cout << c << int_value << endl;
+
+                        int_value += correction_term;
+
+                        // cout << int_value << endl;
+
+                        // cout  << "***************************"  << endl << endl;
+
+                        return int_value;
                         
                 }

applications/robust/robustlib.h

@@ -273 +273 @@
 class c_statistic
 {
+
+public:
         polyhedron* end_poly;
         polyhedron* start_poly;
-
-public: 
 
         vector<polyhedron*> rows;
@@ -567 +567 @@
         /// A statistic in a form of a Hasse diagram representing a complex of convex polyhedrons obtained as a result
         /// of data update from Bayesian estimation or set by the user if this emlig is a prior density
-        c_statistic statistic;
+        
 
         vector<list<polyhedron*>> for_splitting;
@@ -787 +787 @@
         
 public: 
+        c_statistic statistic;
 
         vector<set<my_ivec>> correction_factors;