Changeset 1343 for applications/robust

Timestamp:

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

Author:

sindj

Message:

Importance sampling. JS

Location:

applications/robust

Files:

• main.cpp (modified) (7 diffs)
• robustlib.h (modified) (8 diffs)

applications/robust/main.cpp

@@ -18 +18 @@
 
 const int emlig_size = 2;
+const int utility_constant = 2;
 
 
@@ -138 +139 @@
         vector<vector<string>> strings;
 
-        char* file_strings[3] = {"c:\\dataCDClosePercDiff", "c:\\ar_student_single","c:\\ar_cauchy_single"};
+        char* file_strings[3] = {"c:\\dataADClosePercDiff","c:\\ar_student_single","c:\\ar_cauchy_single"};
 
         for(int i = 0;i<3;i++)
@@ -170 +171 @@
                 vector<vec> conditions;
                 //emlig* emliga = new emlig(2);
-                RARX* my_rarx = new RARX(2,30,false);
+                RARX* my_rarx = new RARX(3,30,true);
                 
                 
@@ -199 +200 @@
                                 conditions[k-3].ins(0,strings[j][k]);
 
-                                //cout << "Condition:" << conditions[k-3] << endl;
+                                cout << "Condition:" << conditions[k-3] << endl;
 
                                 my_rarx->bayes(conditions[k-3]);
@@ -210 +211 @@
                                         
                                 
-                                if(k>8)
+                                if(k>10)
                                 {
                                         //my_rarx->posterior->step_me(0);
 
-                                        mat samples = my_rarx->posterior->sample_mat(50);
+                                        mat samples = my_rarx->posterior->sample_mat(10);
+
+                                        pair<vec,mat> imp_samples = my_rarx->posterior->importance_sample(20);
+
+                                        cout << imp_samples.first << endl;
                                         
                                         vec sample_prediction;
-                                        for(int t = 0;t<50;t++)
+                                        for(int t = 0;t<samples.cols();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());
@@ -229 +234 @@
                                         
                                         vec sample_pow = sample_prediction;
+                                        
+                                        // cout << sample_prediction << endl;
                                         vec poly_coefs;
                                         bool stop_iteration = false;
-                                        int en = 0;
+                                        int en = 1;
                                         do
                                         {
                                                 double poly_coef = ones(sample_pow.size())*sample_pow/sample_pow.size();
+
+                                                poly_coef = poly_coef*en*fact(utility_constant-2+en)/fact(utility_constant-2);
 
                                                 if(abs(poly_coef)>numeric_limits<double>::epsilon())
                                                 {
                                                         sample_pow = elem_mult(sample_pow,sample_prediction);
-                                                        poly_coefs.ins(poly_coefs.size(),((-1)^en)*poly_coef);
+                                                        poly_coefs.ins(0,pow(-1.0,en+1)*poly_coef);
                                                 }
                                                 else
@@ -247 +256 @@
                                                 
                                                 en++;
+
+                                                if(en>20)
+                                                {
+                                                        stop_iteration = true;
+                                                }
                                         }
                                         while(!stop_iteration);
+
+                                        /*
+                                        ofstream myfile_coef;                                           
+
+                                        myfile_coef.open("c:\\coefs.txt",ios::app);
+                                        
+                                        for(int t = 0;t<poly_coefs.size();t++)
+                                        {
+                                                myfile_coef << poly_coefs[t] << ",";                                    
+                                        }
+
+                                        myfile_coef << endl;
+                                        myfile_coef.close();
+                                        */
 
                                         // cout << "Coefficients: " << poly_coefs << endl;
                                                                                 
+                                        /*
+                                        vec bas_coef = vec("1.0 2.0 -8.0");
+                                        cout << "Coefs: " << bas_coef << endl;
+                                        cvec actions2 = roots(bas_coef);
+                                        cout << "Roots: " << actions2 << endl;
+                                        */
+                                        
                                         cvec actions = roots(poly_coefs);
+                                        
+
                                         bool is_max = false;
                                         for(int t = 0;t<actions.size();t++)
                                         {
-                                                if(actions[t].imag() == 0 && actions[t].real()>-1 && actions[t].real()<1)
+                                                if(actions[t].imag() == 0)
                                                 {
-                                                        cout << "Action:" << actions[t].real() << endl;
-                                                        is_max = true;
+                                                        
+
+                                                        double second_derivative = 0;
+                                                        for(int q = 1;q<poly_coefs.size();q++)
+                                                        {
+                                                                second_derivative+=poly_coefs[q]*pow(actions[t].real(),q-1)*q;
+                                                        }
+
+                                                        if(second_derivative<0)
+                                                        {
+                                                                cout << "Action:" << actions[t].real() << endl;
+
+                                                                is_max = true;
+                                                        }
                                                 }
                                         }

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);
         }