Changeset 1346 for applications
- Timestamp:
- 05/02/11 11:27:15 (13 years ago)
- Location:
- applications/robust
- Files:
-
- 2 modified
-
main.cpp (modified) (11 diffs)
-
robustlib.h (modified) (15 diffs)
Legend:
- Unmodified
- Added
- Removed
-
applications/robust/main.cpp
r1343 r1346 18 18 19 19 const int emlig_size = 2; 20 const int utility_constant = 2;20 const int utility_constant = 3; 21 21 22 22 … … 139 139 vector<vector<string>> strings; 140 140 141 char* file_strings[3] = {"c:\\ dataADClosePercDiff","c:\\ar_student_single","c:\\ar_cauchy_single"};141 char* file_strings[3] = {"c:\\Users\\Hontik\\Desktop\\dataGCClosePercDiff","c:\\ar_student_single","c:\\ar_cauchy_single"}; 142 142 143 143 for(int i = 0;i<3;i++) … … 171 171 vector<vec> conditions; 172 172 //emlig* emliga = new emlig(2); 173 RARX* my_rarx = new RARX( 3,30,true);173 RARX* my_rarx = new RARX(2,30,false); 174 174 175 175 … … 211 211 212 212 213 if(k> 10)213 if(k>8) 214 214 { 215 215 //my_rarx->posterior->step_me(0); 216 216 217 mat samples = my_rarx->posterior->sample_mat(10);218 219 pair<vec,mat> imp_samples = my_rarx->posterior->importance_sample( 20);220 221 cout << imp_samples.first << endl;217 //mat samples = my_rarx->posterior->sample_mat(10); 218 219 pair<vec,mat> imp_samples = my_rarx->posterior->importance_sample(1000); 220 221 //cout << imp_samples.first << endl; 222 222 223 223 vec sample_prediction; 224 for(int t = 0;t< samples.cols();t++)224 for(int t = 0;t<imp_samples.first.size();t++) 225 225 { 226 226 vec lap_sample = conditions[k-3].left(2); 227 lap_sample.ins(lap_sample.size(),1.0);227 //lap_sample.ins(lap_sample.size(),1.0); 228 228 229 229 lap_sample.ins(0,LapRNG()); 230 230 231 sample_prediction.ins(0,lap_sample* samples.get_col(t));231 sample_prediction.ins(0,lap_sample*imp_samples.second.get_col(t)); 232 232 } 233 233 … … 237 237 // cout << sample_prediction << endl; 238 238 vec poly_coefs; 239 double prediction; 239 240 bool stop_iteration = false; 240 241 int en = 1; 241 242 do 242 243 { 243 double poly_coef = ones(sample_pow.size())*sample_pow/sample_pow.size(); 244 double poly_coef = imp_samples.first*sample_pow/(imp_samples.first*ones(imp_samples.first.size())); 245 246 if(en==1) 247 { 248 prediction = poly_coef; 249 } 244 250 245 251 poly_coef = poly_coef*en*fact(utility_constant-2+en)/fact(utility_constant-2); … … 278 284 */ 279 285 280 //cout << "Coefficients: " << poly_coefs << endl;286 cout << "Coefficients: " << poly_coefs << endl; 281 287 282 288 /* … … 287 293 */ 288 294 295 296 289 297 cvec actions = roots(poly_coefs); 290 298 … … 295 303 if(actions[t].imag() == 0) 296 304 { 297 298 299 305 double second_derivative = 0; 300 306 for(int q = 1;q<poly_coefs.size();q++) … … 319 325 // cout << "MaxLik coords:" << my_rarx->posterior->minimal_vertex->get_coordinates() << endl; 320 326 327 /* 321 328 double prediction = 0; 322 329 for(int s = 1;s<samples.rows();s++) 323 330 { 324 331 325 double avg_parameter = samples.get_row(s)*ones(samples.cols())/samples.cols();332 double avg_parameter = imp_samples.get_row(s)*ones(samples.cols())/samples.cols(); 326 333 327 334 prediction += avg_parameter*conditions[k-3][s-1]; … … 355 362 myfile.close(); 356 363 */ 357 } 358 359 cout << "Prediction: "<< prediction << endl; 364 365 366 //} 367 368 // cout << "Prediction: "<< prediction << endl; 360 369 361 370 enorm<ldmat>* pred_mat = my_arx->epredictor(conditions[k-3].left(2)); … … 401 410 } 402 411 myfile.close(); 403 412 //*/ 404 413 405 414 } -
applications/robust/robustlib.h
r1343 r1346 847 847 //current_parent->set_state(0,MERGE); 848 848 849 if( level == number_of_parameters - 1)849 if((level == number_of_parameters - 1) && (!shouldsplit)) 850 850 { 851 851 toprow* cur_par_toprow = ((toprow*)current_parent); … … 863 863 } 864 864 865 cur_par_toprow->my_emlig->normalization_factor += cur_par_toprow->probability;865 normalization_factor += cur_par_toprow->probability; 866 866 867 867 //current_parent->triangulation.clear(); … … 873 873 { 874 874 for_merging[level+1].push_back(current_parent); 875 //current_parent->parentconditions.erase(toremove); 875 //current_parent->parentconditions.erase(toremove); 876 876 } 877 877 … … 953 953 ((toprow*)current_parent)->condition_order++; 954 954 955 if(level == number_of_parameters - 1 )955 if(level == number_of_parameters - 1 && current_parent->mergechild == NULL) 956 956 { 957 957 toprow* cur_par_toprow = ((toprow*)current_parent); … … 969 969 } 970 970 971 cur_par_toprow->my_emlig->normalization_factor += cur_par_toprow->probability;971 normalization_factor += cur_par_toprow->probability; 972 972 973 973 //current_parent->triangulation.clear(); … … 1154 1154 vec null_vector = ""; 1155 1155 1156 add_and_remove_condition(null_vector, toremove); 1157 1156 add_and_remove_condition(null_vector, toremove); 1158 1157 } 1159 1158 … … 1338 1337 else 1339 1338 { 1339 bool will_be_split = false; 1340 1340 1341 toprow* current_positive = (toprow*)(*merge_ref)->positiveparent; 1341 1342 toprow* current_negative = (toprow*)(*merge_ref)->negativeparent; … … 1395 1396 current_positive->set_state(0,SPLIT); 1396 1397 for_splitting[k].push_back(current_positive); 1398 will_be_split = true; 1397 1399 } 1398 1400 … … 1428 1430 current_positive->negativeneutralvertices.insert(current_negative->negativeneutralvertices.begin(),current_negative->negativeneutralvertices.end()); 1429 1431 current_positive->positiveneutralvertices.insert(current_negative->positiveneutralvertices.begin(),current_negative->positiveneutralvertices.end()); 1432 1433 will_be_split = true; 1430 1434 } 1431 1435 else … … 1505 1509 current_positive->totally_neutral = (current_positive->totally_neutral && current_negative->totally_neutral); 1506 1510 1507 current_positive->my_emlig->normalization_factor += current_positive->triangulate(k==for_splitting.size()-1);1511 normalization_factor += current_positive->triangulate(k==for_splitting.size()-1 && !will_be_split); 1508 1512 1509 1513 statistic.delete_polyhedron(k,current_negative); … … 1722 1726 new_totally_neutral_child->triangulate(false); 1723 1727 1724 positive_poly->my_emlig->normalization_factor += positive_poly->triangulate(k==for_splitting.size()-1);1725 n egative_poly->my_emlig->normalization_factor += negative_poly->triangulate(k==for_splitting.size()-1);1728 normalization_factor += positive_poly->triangulate(k==for_splitting.size()-1); 1729 normalization_factor += negative_poly->triangulate(k==for_splitting.size()-1); 1726 1730 1727 1731 statistic.append_polyhedron(k-1, new_totally_neutral_child); … … 1822 1826 double prob_sum = 0; 1823 1827 toprow* sampled_toprow; 1824 for(polyhedron* top_ref = statistic.rows[ statistic.size()-1];top_ref!=statistic.row_ends[statistic.size()-1];top_ref=top_ref->next_poly)1828 for(polyhedron* top_ref = statistic.rows[number_of_parameters];top_ref!=statistic.end_poly;top_ref=top_ref->next_poly) 1825 1829 { 1826 1830 // cout << "CDF:"<< (*top_ref).first << endl; … … 1834 1838 sampled_toprow = (toprow*)top_ref; 1835 1839 break; 1836 } 1840 } 1841 else 1842 { 1843 if(top_ref->next_poly==statistic.end_poly) 1844 { 1845 cout << "Error."; 1846 } 1847 } 1837 1848 } 1838 1849 … … 1857 1868 pair<double,double> choose_sigma(simplex* sampled_simplex) 1858 1869 { 1859 double rnumber = randu();1860 double sigma;1861 1862 do uble sum_g = 0;1863 for(int i = 0;i<sampled_simplex->positive_gamma_parameters.size();i++)1864 {1865 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++)1866 {1867 sum_g += (*g_ref).first/sampled_simplex->positive_gamma_sum;1868 1869 1870 if(sum_g>rnumber)1871 { 1872 //itpp::Gamma_RNG* gamma = new itpp::Gamma_RNG(conditions.size()-number_of_parameters,1/(*g_ref).second);1873 //sigma = 1/(*gamma)(); 1870 double sigma = 0; 1871 double pg_sum; 1872 double ng_sum; 1873 do 1874 { 1875 double rnumber = randu(); 1876 1877 1878 double sum_g = 0; 1879 for(int i = 0;i<sampled_simplex->positive_gamma_parameters.size();i++) 1880 { 1881 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++) 1882 { 1883 sum_g += (*g_ref).first/sampled_simplex->positive_gamma_sum; 1884 1874 1885 1875 GamRNG.setup(conditions.size()-number_of_parameters,(*g_ref).second); 1876 1877 sigma = 1/GamRNG(); 1878 1879 // cout << "Sigma mean: " << (*g_ref).second/(conditions.size()-number_of_parameters-1) << endl; 1886 if(sum_g>rnumber) 1887 { 1888 //itpp::Gamma_RNG* gamma = new itpp::Gamma_RNG(conditions.size()-number_of_parameters,1/(*g_ref).second); 1889 //sigma = 1/(*gamma)(); 1890 1891 GamRNG.setup(conditions.size()-number_of_parameters,(*g_ref).second); 1892 1893 sigma = 1/GamRNG(); 1894 1895 // cout << "Sigma mean: " << (*g_ref).second/(conditions.size()-number_of_parameters-1) << endl; 1896 break; 1897 } 1898 } 1899 1900 if(sigma!=0) 1901 { 1880 1902 break; 1881 } 1882 } 1883 1884 if(sigma!=0) 1885 { 1886 break; 1887 } 1888 } 1889 1890 rnumber = randu(); 1891 1892 double pg_sum = 0; 1893 for(vector<multimap<double,double>>::iterator v_ref = sampled_simplex->positive_gamma_parameters.begin();v_ref!=sampled_simplex->positive_gamma_parameters.end();v_ref++) 1894 { 1895 for(multimap<double,double>::iterator pg_ref = (*v_ref).begin();pg_ref!=(*v_ref).end();pg_ref++) 1896 { 1897 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; 1898 } 1899 } 1900 1901 double ng_sum = 0; 1902 for(vector<multimap<double,double>>::iterator v_ref = sampled_simplex->negative_gamma_parameters.begin();v_ref!=sampled_simplex->negative_gamma_parameters.end();v_ref++) 1903 { 1904 for(multimap<double,double>::iterator ng_ref = (*v_ref).begin();ng_ref!=(*v_ref).end();ng_ref++) 1905 { 1906 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; 1907 } 1908 } 1903 } 1904 } 1905 1906 rnumber = randu(); 1907 1908 pg_sum = 0; 1909 for(vector<multimap<double,double>>::iterator v_ref = sampled_simplex->positive_gamma_parameters.begin();v_ref!=sampled_simplex->positive_gamma_parameters.end();v_ref++) 1910 { 1911 for(multimap<double,double>::iterator pg_ref = (*v_ref).begin();pg_ref!=(*v_ref).end();pg_ref++) 1912 { 1913 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; 1914 } 1915 } 1916 1917 ng_sum = 0; 1918 for(vector<multimap<double,double>>::iterator v_ref = sampled_simplex->negative_gamma_parameters.begin();v_ref!=sampled_simplex->negative_gamma_parameters.end();v_ref++) 1919 { 1920 for(multimap<double,double>::iterator ng_ref = (*v_ref).begin();ng_ref!=(*v_ref).end();ng_ref++) 1921 { 1922 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; 1923 } 1924 } 1925 } 1926 while(pg_sum-ng_sum<0); 1909 1927 1910 1928 return pair<double,double>((pg_sum-ng_sum)/pg_sum,sigma); … … 2277 2295 extended_coords.ins(0,-1.0); 2278 2296 2279 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)); 2297 double exponent = extended_coords*condition_and_simplex.first; 2298 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); 2280 2299 sample_prob *= probability_and_sigma.first; 2281 2300
