Changeset 1336 for applications
- Timestamp:
- 04/20/11 20:11:40 (13 years ago)
- Location:
- applications/robust
- Files:
-
- 2 modified
-
main.cpp (modified) (3 diffs)
-
robustlib.h (modified) (17 diffs)
Legend:
- Unmodified
- Added
- Removed
-
applications/robust/main.cpp
r1331 r1336 134 134 vector<vector<string>> strings; 135 135 136 char* file_strings[3] = {"c:\\ar_ cauchy_single.txt", "c:\\ar_normal_single.txt","c:\\ar_student_single.txt"};136 char* file_strings[3] = {"c:\\ar_student_single.txt", "c:\\ar_student_single.txt","c:\\ar_cauchy_single.txt"}; 137 137 138 138 for(int i = 0;i<3;i++) … … 162 162 vector<vec> conditions; 163 163 //emlig* emliga = new emlig(2); 164 RARX* my_rarx = new RARX( 2,70);164 RARX* my_rarx = new RARX(3,20,true); 165 165 166 166 for(int k = 1;k<170;k++) … … 189 189 190 190 191 if(k> 5)191 if(k>10) 192 192 { 193 193 //my_rarx->posterior->step_me(0); 194 194 195 my_rarx->posterior->sample_mat(1); 196 195 mat samples = my_rarx->posterior->sample_mat(50); 196 197 197 198 cout << "MaxLik coords:" << my_rarx->posterior->minimal_vertex->get_coordinates() << endl; 198 199 199 ofstream myfile; 200 char fstring[80]; 201 strcpy(fstring,file_strings[j]); 202 strcat(fstring,"_res.txt"); 203 204 myfile.open(fstring,ios::app); 205 myfile << my_rarx->posterior->minimal_vertex->get_coordinates()[0]; 206 if(k!=strings[j].size()-1) 200 for(int s = 0;s<samples.rows();s++) 207 201 { 208 myfile << ","; 202 203 double avg_parameter = samples.get_row(s)*ones(samples.cols())/samples.cols(); 204 205 ofstream myfile; 206 char fstring[80]; 207 strcpy(fstring,file_strings[j]); 208 209 char es[5]; 210 strcat(fstring,itoa(s,es,10)); 211 212 strcat(fstring,"_res.txt"); 213 214 215 myfile.open(fstring,ios::app); 216 217 //myfile << my_rarx->posterior->minimal_vertex->get_coordinates()[0]; 218 myfile << avg_parameter; 219 220 if(k!=strings[j].size()-1) 221 { 222 myfile << ","; 223 } 224 else 225 { 226 myfile << endl; 227 } 228 myfile.close(); 209 229 } 210 else211 {212 myfile << endl;213 }214 myfile.close();215 230 } 216 231 } -
applications/robust/robustlib.h
r1335 r1336 17 17 #include <algorithm> 18 18 19 //using namespace bdm;19 using namespace bdm; 20 20 using namespace std; 21 21 using namespace itpp; … … 1831 1831 while(sample_mat.cols()<n) 1832 1832 { 1833 cout << "*************************************" << endl; 1833 //// cout << "*************************************" << endl; 1834 1835 1834 1836 1835 1837 double rnumber = randu()*sum_a; … … 1852 1854 } 1853 1855 1854 cout << "Toprow/Count: " << toprow_count << "/" << ordered_toprows.size() << endl;1856 //// cout << "Toprow/Count: " << toprow_count << "/" << ordered_toprows.size() << endl; 1855 1857 // cout << &sampled_toprow << ";"; 1856 1858 … … 1870 1872 } 1871 1873 1872 cout << "Simplex/Count: " << simplex_count << "/" << sampled_toprow->triangulation.size() << endl;1873 cout << "Simplex factor: " << (*s_ref)->probability << endl;1874 cout << "Toprow factor: " << sampled_toprow->probability << endl;1875 cout << "Emlig factor: " << normalization_factor << endl;1874 //// cout << "Simplex/Count: " << simplex_count << "/" << sampled_toprow->triangulation.size() << endl; 1875 //// cout << "Simplex factor: " << (*s_ref)->probability << endl; 1876 //// cout << "Toprow factor: " << sampled_toprow->probability << endl; 1877 //// cout << "Emlig factor: " << normalization_factor << endl; 1876 1878 // cout << &(*tri_ref) << endl; 1877 1879 … … 1893 1895 if(sum_g>rnumber) 1894 1896 { 1895 itpp::Gamma_RNG* gamma = new itpp::Gamma_RNG(conditions.size()-number_of_parameters,1/(*g_ref).second);1896 sigma = 1/(*gamma)();1897 //itpp::Gamma_RNG* gamma = new itpp::Gamma_RNG(conditions.size()-number_of_parameters,1/(*g_ref).second); 1898 //sigma = 1/(*gamma)(); 1897 1899 1898 cout << "Sigma mean: " << (*g_ref).second/(conditions.size()-number_of_parameters-1) << endl; 1900 GamRNG.setup(conditions.size()-number_of_parameters,(*g_ref).second); 1901 1902 sigma = 1/GamRNG(); 1903 1904 // cout << "Sigma mean: " << (*g_ref).second/(conditions.size()-number_of_parameters-1) << endl; 1899 1905 break; 1900 1906 } … … 1936 1942 while(!have_sigma); 1937 1943 1938 cout << "Sigma: " << sigma << endl;1939 cout << "Nr. of runs: " << number_of_runs << endl;1944 //// cout << "Sigma: " << sigma << endl; 1945 //// cout << "Nr. of runs: " << number_of_runs << endl; 1940 1946 1941 1947 int dimension = (*s_ref)->vertices.size()-1; … … 1946 1952 vertex* base_vert = *(*s_ref)->vertices.begin(); 1947 1953 1948 cout << "Base vertex coords(should be close to est. param.): " << base_vert->get_coordinates() << endl;1954 //// cout << "Base vertex coords(should be close to est. param.): " << base_vert->get_coordinates() << endl; 1949 1955 1950 1956 int row_count = 0; … … 1954 1960 vec current_coords = (*vert_ref)->get_coordinates(); 1955 1961 1956 cout << "Coords of vertex[" << row_count << "]: " << current_coords << endl;1962 //// cout << "Coords of vertex[" << row_count << "]: " << current_coords << endl; 1957 1963 1958 1964 vec relative_coords = current_coords-base_vert->get_coordinates(); … … 1963 1969 } 1964 1970 1965 cout << "Jacobian: " << jacobian << endl;1966 1967 cout << "Gradient before trafo:" << gradient << endl;1971 //// cout << "Jacobian: " << jacobian << endl; 1972 1973 //// cout << "Gradient before trafo:" << gradient << endl; 1968 1974 1969 1975 gradient = jacobian*gradient; 1970 1976 1971 cout << "Gradient after trafo:" << gradient << endl;1977 //// cout << "Gradient after trafo:" << gradient << endl; 1972 1978 1973 1979 // vec normal_gradient = gradient/sqrt(gradient*gradient); … … 2009 2015 extended_rotation.ins_col(0,zeros(extended_rotation.rows())); 2010 2016 2011 cout << "Extended rotation: " << extended_rotation << endl;2017 //// cout << "Extended rotation: " << extended_rotation << endl; 2012 2018 2013 2019 vec minima = itpp::min(extended_rotation,2); 2014 2020 vec maxima = itpp::max(extended_rotation,2); 2015 2021 2016 cout << "Minima: " << minima << endl;2017 cout << "Maxima: " << maxima << endl;2022 //// cout << "Minima: " << minima << endl; 2023 //// cout << "Maxima: " << maxima << endl; 2018 2024 2019 2025 vec sample_coordinates; … … 2033 2039 vec new_gradient = rotation_matrix*gradient; 2034 2040 2035 cout << "New gradient(should have only first component nonzero):" << new_gradient << endl;2041 //// cout << "New gradient(should have only first component nonzero):" << new_gradient << endl; 2036 2042 2037 2043 // cout << "Max: " << maxima[0] << " Min: " << minima[0] << " Grad:" << new_gradient[0] << endl; … … 2049 2055 } 2050 2056 2051 cout << "Sampled coordinates(gradient direction): " << sample_coordinates << endl;2057 //// cout << "Sampled coordinates(gradient direction): " << sample_coordinates << endl; 2052 2058 2053 2059 sample_coordinates = rotation_matrix.T()*sample_coordinates; 2054 2060 2055 cout << "Sampled coordinates(backrotated direction):" << sample_coordinates << endl;2061 //// cout << "Sampled coordinates(backrotated direction):" << sample_coordinates << endl; 2056 2062 2057 2063 … … 2077 2083 sample_coordinates.ins(0,sigma); 2078 2084 2079 cout << "Sampled coordinates(parameter space):" << sample_coordinates << endl;2085 //// cout << "Sampled coordinates(parameter space):" << sample_coordinates << endl; 2080 2086 2081 2087 sample_mat.ins_col(0,sample_coordinates); 2088 2089 cout << sample_mat.cols() << ","; 2082 2090 } 2083 2091 … … 2088 2096 } 2089 2097 2098 cout << endl; 2090 2099 return sample_mat; 2091 2100 } … … 2405 2414 { 2406 2415 private: 2407 2408 2416 bool has_constant; 2409 2417 2410 2418 int window_size; … … 2415 2423 emlig* posterior; 2416 2424 2417 RARX(int number_of_parameters, const int window_size)//:BM() 2418 { 2425 RARX(int number_of_parameters, const int window_size, bool has_constant)//:BM() 2426 { 2427 this->has_constant = has_constant; 2428 2419 2429 posterior = new emlig(number_of_parameters); 2420 2430 … … 2422 2432 }; 2423 2433 2424 void bayes(const itpp::vec &yt, const itpp::vec &cond = "") 2425 { 2426 conditions.push_back(yt); 2434 void bayes(itpp::vec yt) 2435 { 2436 if(has_constant) 2437 { 2438 int c_size = yt.size(); 2439 2440 yt.ins(c_size,1.0); 2441 } 2442 2443 if(yt.size() == posterior->number_of_parameters+1) 2444 { 2445 conditions.push_back(yt); 2446 } 2447 else 2448 { 2449 throw new exception("Wrong condition size for bayesian data update!"); 2450 } 2427 2451 2428 2452 //posterior->step_me(0);
