Changeset 1362 for applications

Timestamp:

05/09/11 19:07:15 (14 years ago)

Author:

sindj

Message:

Pridani soft apriorna, boj s integraly. JS

Location:

applications/robust

Files:

• main.cpp (modified) (2 diffs)
• robustlib.cpp (modified) (3 diffs)
• robustlib.h (modified) (16 diffs)

applications/robust/main.cpp

@@ -111 +111 @@
                         if(has_constant)
                         {                               
-                                V0  = 0.0001 * eye(ar_components.size()+2);
-                                //V0(0,0) = 0.1;
+                                V0  = 0.01 * eye(ar_components.size()+2);
+                                V0(0,0) = 1;
                                 my_arx->set_constant(true);     
                                 
@@ -119 +119 @@
                         {
                                 
-                                V0  = 0.0001 * eye(ar_components.size()+1);
-                                //V0(0,0) = 0.1;
+                                V0  = 0.01 * eye(ar_components.size()+1);
+                                V0(0,0) = 1;
                                 my_arx->set_constant(false);                            
                                 

applications/robust/robustlib.cpp

@@ -96 +96 @@
         {
                 // cout << ((toprow*)this)->condition << endl;          
-
-                int condition_order = ((toprow*)this)->condition_order-2; // -2 by bylo, pokud chceme uniformni apriorno
-
-                int sigma_order = condition_order-my_emlig->number_of_parameters;
+                int sigma_order = ((toprow*)this)->condition_order-simplex->vertices.size()-1;
 
                 // cout << "C:" << condition_order << "  N:" << my_emlig->number_of_parameters << "  C+N:" << condition_order-my_emlig->number_of_parameters << endl;
@@ -147 +144 @@
 
                                 base_vertex = (*base_ref);
-
-                        
-                                // cout << endl << "Base coords:" << base_vertex->get_coordinates() << endl;
+                        
+                                //cout << "Base coords: " << base_vertex->get_coordinates() << endl;
+                                //cout << "Condition: " << as_toprow->condition_sum << endl;
+                                //pause(10);
 
                                 a_0 = base_vertex->get_coordinates()*cur_condition-as_toprow->condition_sum[0];
@@ -340 +338 @@
                 else
                 {
+                        cout << "Improper probability density." << endl;
                         return 0.0;
                 }

applications/robust/robustlib.h

@@ -22 +22 @@
 using namespace itpp;
 
-const double max_range = 20;//numeric_limits<double>::max()/10e-10;
+const double max_range = 5;//numeric_limits<double>::max()/10e-10;
 
 /// An enumeration of possible actions performed on the polyhedrons. We can merge them or split them.
@@ -740 +740 @@
         int condition_order;
 
+        double last_log_nc;
+
         
 
@@ -826 +828 @@
                                         if(is_last)
                                         {                                               
-                                                if(parent_state == 1)
-                                                {
-                                                        ((toprow*)current_parent)->condition_sum-=toremove->value;                                                      
-                                                }
-
-                                                if(parent_state == -1)
-                                                {
-                                                        ((toprow*)current_parent)->condition_sum+=toremove->value;                                                      
+                                                if(level == number_of_parameters-1)
+                                                {
+                                                        if(parent_state == 1)
+                                                        {
+                                                                ((toprow*)current_parent)->condition_sum-=toremove->value;                                                      
+                                                        }
+
+                                                        if(parent_state == -1)
+                                                        {
+                                                                ((toprow*)current_parent)->condition_sum+=toremove->value;                                                      
+                                                        }
                                                 }
 
@@ -964 +969 @@
                                                                 current_parent->set_state(-1, SPLIT);
 
-                                                                ((toprow*)current_parent)->condition_sum-=toadd->value;
+                                                                if(level == number_of_parameters-1)
+                                                                {
+                                                                        ((toprow*)current_parent)->condition_sum-=toadd->value;
+                                                                }
                                                                         
                                                         }
@@ -971 +979 @@
                                                                 current_parent->set_state(1, SPLIT);
 
-                                                                ((toprow*)current_parent)->condition_sum+=toadd->value;                                                                 
+                                                                if(level == number_of_parameters-1)
+                                                                {
+                                                                        ((toprow*)current_parent)->condition_sum+=toadd->value;                                                                 
+                                                                }
                                                         }
                                                         else
@@ -1044 +1055 @@
         double log_nc;
 
+        
+
         vector<multiset<my_ivec>> correction_factors;
 
@@ -1050 +1063 @@
         /// A default constructor creates an emlig with predefined statistic representing only the range of the given
         /// parametric space, where the number of parameters of the needed model is given as a parameter to the constructor.
-        emlig(int number_of_parameters)
+        emlig(int number_of_parameters, bool is_proper, double soft_prior_parameter)
         {       
                 this->number_of_parameters = number_of_parameters;
 
-                create_statistic(number_of_parameters);
-
-                min_ll = numeric_limits<double>::max();
-
-                condition_order = 0;
+                condition_order = number_of_parameters+2;
+
+                if(is_proper)
+                        condition_order++;
+                
+                create_statistic(number_of_parameters, soft_prior_parameter);
+
+                //step_me(10);
+
+                min_ll = numeric_limits<double>::max();         
+
+                if(is_proper)
+                {
+                        double normalization_factor = 0;
+                        int counter = 0;
+                        for(polyhedron* top_ref = statistic.rows[number_of_parameters];top_ref!=statistic.get_end();top_ref=top_ref->next_poly)
+                        {
+                                counter++;
+                                toprow* cur_toprow = (toprow*)top_ref;
+                                
+                                set<simplex*>::iterator cur_simplex = cur_toprow->triangulation.begin();
+                                normalization_factor += cur_toprow->integrate_simplex(*cur_simplex,'X');
+                        }
+
+                        log_nc = log(normalization_factor) + logfact(condition_order-number_of_parameters-2);
+
+                        /*
+                        cout << "part1: " << log(normalization_factor) << endl;
+                        cout << "part2: " << logfact(condition_order-number_of_parameters-2) << endl;
+                        pause(1);
+                        */
+                }
+
         }
 
@@ -1085 +1126 @@
                         {
                                 
-                                /*
+                                
                                 if(i==statistic.size()-1)
                                 {
                                         cout << ((toprow*)horiz_ref)->condition_sum << "   " << ((toprow*)horiz_ref)->probability << endl;
+                                        cout << "Condition: " << ((toprow*)horiz_ref)->condition_sum << endl;
                                         cout << "Order:" << ((toprow*)horiz_ref)->condition_order << endl;
                                 }
-                                */
+                                
 
                                 // cout << "Stepped." << endl;
@@ -1212 +1254 @@
         void add_and_remove_condition(vec toadd, vec toremove)
         {
+                
                 //step_me(0);
                 normalization_factor = 0;
@@ -1717 +1760 @@
 
                                         // cout << ((toprow*)current_polyhedron)->condition << endl << toadd << endl;
-
-                                        toprow* positive_poly = new toprow(((toprow*)current_polyhedron)->condition_sum+toadd, ((toprow*)current_polyhedron)->condition_order+1);
-                                        toprow* negative_poly = new toprow(((toprow*)current_polyhedron)->condition_sum-toadd, ((toprow*)current_polyhedron)->condition_order+1);
+                                        vec cur_pos_condition = ((toprow*)current_polyhedron)->condition_sum;
+                                        vec cur_neg_condition = ((toprow*)current_polyhedron)->condition_sum;
+                                        
+                                        if(k == number_of_parameters)
+                                        {
+                                                cur_pos_condition = cur_pos_condition + toadd;
+                                                cur_neg_condition = cur_neg_condition - toadd;
+                                        }
+
+                                        toprow* positive_poly = new toprow(cur_pos_condition, ((toprow*)current_polyhedron)->condition_order+1);
+                                        toprow* negative_poly = new toprow(cur_neg_condition, ((toprow*)current_polyhedron)->condition_order+1);
 
                                         positive_poly->my_emlig = this;
@@ -1835 +1886 @@
                 log_nc = log(normalization_factor) + logfact(condition_order-number_of_parameters-2);
 
+                
+                cout << "part1: " << log(normalization_factor) << endl;
+                cout << "part2: " << logfact(condition_order-number_of_parameters-2) << endl;
+                pause(1);
+                
+
                 /*
                 for(polyhedron* topr_ref = statistic.rows[statistic.size()-1];topr_ref!=statistic.row_ends[statistic.size()-1]->next_poly;topr_ref=topr_ref->next_poly)
@@ -1842 +1899 @@
                 */
 
-                step_me(101);
+                // step_me(101);
 
         }
@@ -2403 +2460 @@
 
         /// A method for creating plain default statistic representing only the range of the parameter space.
-    void create_statistic(int number_of_parameters)
+    void create_statistic(int number_of_parameters, double soft_prior_parameter)
         {
                 /*
@@ -2448 +2505 @@
                         // of newly added parameter. Therefore they will have all coordinates except the last one zero. We get the 
                         // right amount of zero cooridnates by reading them from the origin
-                        vec origin_coord = origin->get_coordinates();                                           
+                        vec origin_coord = origin->get_coordinates();   
+
+                        
 
                         // And we incorporate the nonzero coordinates into the new cooordinate vectors
@@ -2534 +2593 @@
                                         // this condition will be a vector of zeros. There are two vectors, because we need two copies of 
                                         // the original Hasse diagram.
-                                        vec vec1(number_of_parameters+1);
-                                        vec1.zeros();
-
-                                        vec vec2(number_of_parameters+1);
-                                        vec2.zeros();
+                                        vec vec1;
+                                        vec vec2;
+                                        if(!horiz_ref->kids_rel_addresses.empty())
+                                        {                                       
+                                                vec1 = ((toprow*)horiz_ref)->condition_sum;
+                                                vec1.ins(vec1.size(),-soft_prior_parameter);
+
+                                                vec2 = ((toprow*)horiz_ref)->condition_sum;
+                                                vec2.ins(vec2.size(),soft_prior_parameter);
+                                        }
+                                        else
+                                        {                                               
+                                                vec1.ins(0,-soft_prior_parameter);
+                                                vec2.ins(0,soft_prior_parameter);
+
+                                                vec1.ins(0,-1);
+                                                vec2.ins(0,-1);
+                                        }
+                                        
+                                        cout << vec1 << endl;
+                                        cout << vec2 << endl;
+
 
                                         // We create a new toprow with the previously specified condition.
-                                        toprow* current_copy1 = new toprow(vec1, 0);
-                                        toprow* current_copy2 = new toprow(vec2, 0);
+                                        toprow* current_copy1 = new toprow(vec1, this->condition_order);
+                                        toprow* current_copy2 = new toprow(vec2, this->condition_order);
 
                                         current_copy1->my_emlig = this;
@@ -2715 +2791 @@
                 this->has_constant = has_constant;
                 
-                posterior = new emlig(number_of_parameters);
+                posterior = new emlig(number_of_parameters,true,0.001);
 
                 this->window_size = window_size;