Changeset 1346

Timestamp:

05/02/11 11:27:15 (14 years ago)

Author:

sindj

Message:

Oprava integrace, pokud se merguje i splituje. JS

Location:

applications/robust

Files:

• main.cpp (modified) (11 diffs)
• robustlib.h (modified) (15 diffs)

applications/robust/main.cpp

@@ -18 +18 @@
 
 const int emlig_size = 2;
-const int utility_constant = 2;
+const int utility_constant = 3;
 
 
@@ -139 +139 @@
         vector<vector<string>> strings;
 
-        char* file_strings[3] = {"c:\\dataADClosePercDiff","c:\\ar_student_single","c:\\ar_cauchy_single"};
+        char* file_strings[3] = {"c:\\Users\\Hontik\\Desktop\\dataGCClosePercDiff","c:\\ar_student_single","c:\\ar_cauchy_single"};
 
         for(int i = 0;i<3;i++)
@@ -171 +171 @@
                 vector<vec> conditions;
                 //emlig* emliga = new emlig(2);
-                RARX* my_rarx = new RARX(3,30,true);
+                RARX* my_rarx = new RARX(2,30,false);
                 
                 
@@ -211 +211 @@
                                         
                                 
-                                if(k>10)
+                                if(k>8)
                                 {
                                         //my_rarx->posterior->step_me(0);
 
-                                        mat samples = my_rarx->posterior->sample_mat(10);
-
-                                        pair<vec,mat> imp_samples = my_rarx->posterior->importance_sample(20);
-
-                                        cout << imp_samples.first << endl;
+                                        //mat samples = my_rarx->posterior->sample_mat(10);
+
+                                        pair<vec,mat> imp_samples = my_rarx->posterior->importance_sample(1000);
+
+                                        //cout << imp_samples.first << endl;
                                         
                                         vec sample_prediction;
-                                        for(int t = 0;t<samples.cols();t++)
+                                        for(int t = 0;t<imp_samples.first.size();t++)
                                         {
                                                 vec lap_sample = conditions[k-3].left(2);
-                                                lap_sample.ins(lap_sample.size(),1.0);
+                                                //lap_sample.ins(lap_sample.size(),1.0);
                                                 
                                                 lap_sample.ins(0,LapRNG());
 
-                                                sample_prediction.ins(0,lap_sample*samples.get_col(t));
+                                                sample_prediction.ins(0,lap_sample*imp_samples.second.get_col(t));
                                         }
 
@@ -237 +237 @@
                                         // cout << sample_prediction << endl;
                                         vec poly_coefs;
+                                        double prediction;
                                         bool stop_iteration = false;
                                         int en = 1;
                                         do
                                         {
-                                                double poly_coef = ones(sample_pow.size())*sample_pow/sample_pow.size();
+                                                double poly_coef = imp_samples.first*sample_pow/(imp_samples.first*ones(imp_samples.first.size()));
+
+                                                if(en==1)
+                                                {
+                                                        prediction = poly_coef;
+                                                }
 
                                                 poly_coef = poly_coef*en*fact(utility_constant-2+en)/fact(utility_constant-2);
@@ -278 +284 @@
                                         */
 
-                                        // cout << "Coefficients: " << poly_coefs << endl;
+                                        cout << "Coefficients: " << poly_coefs << endl;
                                                                                 
                                         /*
@@ -287 +293 @@
                                         */
                                         
+
+
                                         cvec actions = roots(poly_coefs);
                                         
@@ -295 +303 @@
                                                 if(actions[t].imag() == 0)
                                                 {
-                                                        
-
                                                         double second_derivative = 0;
                                                         for(int q = 1;q<poly_coefs.size();q++)
@@ -319 +325 @@
                                         // cout << "MaxLik coords:" << my_rarx->posterior->minimal_vertex->get_coordinates() << endl;
 
+                                        /*
                                         double prediction = 0;
                                         for(int s = 1;s<samples.rows();s++)
                                         {
                                                 
-                                                double avg_parameter = samples.get_row(s)*ones(samples.cols())/samples.cols();
+                                                double avg_parameter = imp_samples.get_row(s)*ones(samples.cols())/samples.cols();
 
                                                 prediction += avg_parameter*conditions[k-3][s-1];
@@ -355 +362 @@
                                                 myfile.close();
                                                 */
-                                        }
-
-                                        cout << "Prediction: "<< prediction << endl;
+
+                                        
+                                        //}
+
+                                        // cout << "Prediction: "<< prediction << endl;
                                         
                                         enorm<ldmat>* pred_mat = my_arx->epredictor(conditions[k-3].left(2));
@@ -401 +410 @@
                                         }
                                         myfile.close();
-                                        
+                                        //*/
 
                                 }                                       

applications/robust/robustlib.h

@@ -847 +847 @@
                                                         //current_parent->set_state(0,MERGE);   
 
-                                                        if(level == number_of_parameters - 1)
+                                                        if((level == number_of_parameters - 1) && (!shouldsplit))
                                                         {
                                                                 toprow* cur_par_toprow = ((toprow*)current_parent);
@@ -863 +863 @@
                                                                 }
 
-                                                                cur_par_toprow->my_emlig->normalization_factor += cur_par_toprow->probability;
+                                                                normalization_factor += cur_par_toprow->probability;
 
                                                                 //current_parent->triangulation.clear();
@@ -873 +873 @@
                                                 {
                                                         for_merging[level+1].push_back(current_parent);
-                                                        //current_parent->parentconditions.erase(toremove);
+                                                        //current_parent->parentconditions.erase(toremove);                                                     
                                                 }                                               
 
@@ -953 +953 @@
                                                         ((toprow*)current_parent)->condition_order++;
 
-                                                        if(level == number_of_parameters - 1)
+                                                        if(level == number_of_parameters - 1 && current_parent->mergechild == NULL)
                                                         {
                                                                 toprow* cur_par_toprow = ((toprow*)current_parent);
@@ -969 +969 @@
                                                                 }
 
-                                                                cur_par_toprow->my_emlig->normalization_factor += cur_par_toprow->probability;
+                                                                normalization_factor += cur_par_toprow->probability;
 
                                                                 //current_parent->triangulation.clear();
@@ -1154 +1154 @@
                 vec null_vector = "";
 
-                add_and_remove_condition(null_vector, toremove);
-        
+                add_and_remove_condition(null_vector, toremove);        
         }
 
@@ -1338 +1337 @@
                                         else
                                         {
+                                                bool will_be_split = false;
+                                                
                                                 toprow* current_positive = (toprow*)(*merge_ref)->positiveparent;
                                                 toprow* current_negative = (toprow*)(*merge_ref)->negativeparent;
@@ -1395 +1396 @@
                                                         current_positive->set_state(0,SPLIT);
                                                         for_splitting[k].push_back(current_positive);
+                                                        will_be_split = true;
                                                 }
                                                 
@@ -1428 +1430 @@
                                                         current_positive->negativeneutralvertices.insert(current_negative->negativeneutralvertices.begin(),current_negative->negativeneutralvertices.end());
                                                         current_positive->positiveneutralvertices.insert(current_negative->positiveneutralvertices.begin(),current_negative->positiveneutralvertices.end());
+
+                                                        will_be_split = true;
                                                 }
                                                 else
@@ -1505 +1509 @@
                                                 current_positive->totally_neutral = (current_positive->totally_neutral && current_negative->totally_neutral);
 
-                                                current_positive->my_emlig->normalization_factor += current_positive->triangulate(k==for_splitting.size()-1);
+                                                normalization_factor += current_positive->triangulate(k==for_splitting.size()-1 && !will_be_split);
                                                 
                                                 statistic.delete_polyhedron(k,current_negative);
@@ -1722 +1726 @@
                                         new_totally_neutral_child->triangulate(false);
 
-                                        positive_poly->my_emlig->normalization_factor += positive_poly->triangulate(k==for_splitting.size()-1);
-                                        negative_poly->my_emlig->normalization_factor += negative_poly->triangulate(k==for_splitting.size()-1);
+                                        normalization_factor += positive_poly->triangulate(k==for_splitting.size()-1);
+                                        normalization_factor += negative_poly->triangulate(k==for_splitting.size()-1);
                                         
                                         statistic.append_polyhedron(k-1, new_totally_neutral_child);                                    
@@ -1822 +1826 @@
                 double  prob_sum     = 0;       
                 toprow* sampled_toprow;                         
-                for(polyhedron* top_ref = statistic.rows[statistic.size()-1];top_ref!=statistic.row_ends[statistic.size()-1];top_ref=top_ref->next_poly)
+                for(polyhedron* top_ref = statistic.rows[number_of_parameters];top_ref!=statistic.end_poly;top_ref=top_ref->next_poly)
                 {
                         // cout << "CDF:"<< (*top_ref).first << endl;
@@ -1834 +1838 @@
                                 sampled_toprow = (toprow*)top_ref;
                                 break;
-                        }                                               
+                        }
+                        else
+                        {
+                                if(top_ref->next_poly==statistic.end_poly)
+                                {
+                                        cout << "Error.";
+                                }
+                        }
                 }                               
 
@@ -1857 +1868 @@
         pair<double,double> choose_sigma(simplex* sampled_simplex)
         {
-                double rnumber = randu();
-                double sigma;
-                                        
-                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)();
+                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;
+
                                                                 
-                                        GamRNG.setup(conditions.size()-number_of_parameters,(*g_ref).second);
-                                                                                                                                
-                                        sigma = 1/GamRNG();
-
-                                        // cout << "Sigma mean:   " << (*g_ref).second/(conditions.size()-number_of_parameters-1) << endl;                                                              
+                                        if(sum_g>rnumber)
+                                        {
+                                                //itpp::Gamma_RNG* gamma = new itpp::Gamma_RNG(conditions.size()-number_of_parameters,1/(*g_ref).second);
+                                                //sigma = 1/(*gamma)();
+                                                                        
+                                                GamRNG.setup(conditions.size()-number_of_parameters,(*g_ref).second);
+                                                                                                                                        
+                                                sigma = 1/GamRNG();
+
+                                                // cout << "Sigma mean:   " << (*g_ref).second/(conditions.size()-number_of_parameters-1) << endl;                                                              
+                                                break;
+                                        }                                                       
+                                }
+
+                                if(sigma!=0)
+                                {
                                         break;
-                                }                                                       
-                        }
-
-                        if(sigma!=0)
-                        {
-                                break;
-                        }
-                }
-
-                rnumber = randu();
-
-                double 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()-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 = 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()-number_of_parameters-1)*(*ng_ref).second/fact(conditions.size()-number_of_parameters-1)*(*ng_ref).first;
-                        }                                       
-                }
+                                }
+                        }
+
+                        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()-number_of_parameters-1)*(*pg_ref).second/fact(conditions.size()-number_of_parameters-1)*(*pg_ref).first;
+                                }                                       
+                        }
+
+                        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()-number_of_parameters-1)*(*ng_ref).second/fact(conditions.size()-number_of_parameters-1)*(*ng_ref).first;
+                                }                                       
+                        }
+                }
+                while(pg_sum-ng_sum<0);
 
                 return pair<double,double>((pg_sum-ng_sum)/pg_sum,sigma);
@@ -2277 +2295 @@
                         extended_coords.ins(0,-1.0);
 
-                        double sample_prob = 1/condition_and_simplex.second->probability/pow(probability_and_sigma.second,(int)conditions.size()-number_of_parameters)*exp(1/probability_and_sigma.second*(extended_coords*condition_and_simplex.first));
+                        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)*exp((-1)/probability_and_sigma.second*exponent);
                         sample_prob *= probability_and_sigma.first;