Changeset 1343 for applications/robust/robustlib.h

Timestamp:

04/29/11 13:26:04 (13 years ago)

Author:

sindj

Message:

Importance sampling. JS

Files:

• applications/robust/robustlib.h (modified) (8 diffs)

applications/robust/robustlib.h

@@ -113 +113 @@
                 if(weight>=0)
                 {
-                        if(positive_gamma_parameters.size()<order+1)
+                        while(positive_gamma_parameters.size()<order+1)
                         {
                                 multimap<double,double> map;
@@ -125 +125 @@
                 else
                 {
-                        if(negative_gamma_parameters.size()<order+1)
+                        while(negative_gamma_parameters.size()<order+1)
                         {
                                 multimap<double,double> map;
@@ -1577 +1577 @@
                         sizevector.push_back(statistic.row_size(s));
                         cout << statistic.row_size(s) << ", ";
-                }
-                */
+                }*/
+                
 
                 //cout << endl;
@@ -1739 +1739 @@
                 }
 
-                /*
-                sizevector.clear();
+                vector<int> sizevector;
+                //sizevector.clear();
                 for(int s = 0;s<statistic.size();s++)
                 {
@@ -1748 +1748 @@
                 
                 cout << endl;
-                */
+                
 
                 /*
@@ -1813 +1813 @@
                 }
 
-
+        pair<vec,simplex*> choose_simplex()
+        {
+                double rnumber = randu();
+
+                // cout << "RND:" << rnumber << endl;
+
+                // This could be more efficient (log n), but map::upper_bound() doesn't let me dereference returned iterator
+                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)
+                {
+                        // cout << "CDF:"<< (*top_ref).first << endl;
+
+                        toprow* current_toprow = ((toprow*)top_ref);
+
+                        prob_sum += current_toprow->probability;
+
+                        if(prob_sum >= rnumber*normalization_factor)
+                        {
+                                sampled_toprow = (toprow*)top_ref;
+                                break;
+                        }                                               
+                }                               
+
+                //// cout << "Toprow/Count: " << toprow_count << "/" << ordered_toprows.size() << endl;
+                // cout << &sampled_toprow << ";";
+
+                rnumber = randu();                              
+
+                set<simplex*>::iterator s_ref;
+                prob_sum = 0;           
+                for(s_ref = sampled_toprow->triangulation.begin();s_ref!=sampled_toprow->triangulation.end();s_ref++)
+                {               
+                        prob_sum += (*s_ref)->probability;
+
+                        if(prob_sum/sampled_toprow->probability >= rnumber)
+                                break;
+                }
+
+                return pair<vec,simplex*>(sampled_toprow->condition_sum,*s_ref);        
+        }
+
+        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)();
+                                                                
+                                        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;
+                        }
+                }
+
+                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;
+                        }                                       
+                }
+
+                return pair<double,double>((pg_sum-ng_sum)/pg_sum,sigma);
+        }
 
         mat sample_mat(int n)
@@ -2114 +2210 @@
                         }
 
-                        // cout << endl;
+                        cout << endl;
                         return sample_mat;
                 }
@@ -2125 +2221 @@
 
                 
+        }
+
+        pair<vec,mat> importance_sample(int n)
+        {
+                vec probabilities;
+                mat samples;
+                
+                for(int i = 0;i<n;i++)
+                {
+                        pair<vec,simplex*> condition_and_simplex = choose_simplex();
+
+                        pair<double,double> probability_and_sigma = choose_sigma(condition_and_simplex.second);
+
+                        int dimension = condition_and_simplex.second->vertices.size()-1;
+
+                        mat jacobian(dimension,dimension);
+                        vec gradient = condition_and_simplex.first.right(dimension);
+
+                        vertex* base_vert = *condition_and_simplex.second->vertices.begin();
+
+                        //// cout << "Base vertex coords(should be close to est. param.): " << base_vert->get_coordinates() << endl;
+                                
+                        int row_count = 0;
+
+                        for(set<vertex*>::iterator vert_ref = ++condition_and_simplex.second->vertices.begin();vert_ref!=condition_and_simplex.second->vertices.end();vert_ref++)
+                        {
+                                vec current_coords = (*vert_ref)->get_coordinates();
+
+                                //// cout << "Coords of vertex[" << row_count << "]: " << current_coords << endl; 
+                                        
+                                vec relative_coords = current_coords-base_vert->get_coordinates();                              
+
+                                jacobian.set_row(row_count,relative_coords);
+
+                                row_count++;
+                        }                               
+                                
+                        //// cout << "Jacobian: " << jacobian << endl;                  
+
+                        /// \todo Is this correct? Are the random coordinates really jointly uniform? I don't know.
+                        vec sample_coords;
+                        double sampling_diff = 1;
+                        for(int j = 0;j<number_of_parameters;j++)
+                        {
+                                double rnumber = randu()*sampling_diff;
+
+                                sample_coords.ins(0,rnumber);
+
+                                sampling_diff -= rnumber;
+                        }
+
+                        sample_coords = jacobian.T()*sample_coords+(*base_vert).get_coordinates();
+
+                        vec extended_coords = sample_coords;
+                        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));
+                        sample_prob *= probability_and_sigma.first;
+
+                        sample_coords.ins(0,probability_and_sigma.second);
+
+                        samples.ins_col(0,sample_coords);
+                        probabilities.ins(0,sample_prob);
+                }
+        
+                return pair<vec,mat>(probabilities,samples);
         }