Changeset 1336 for applications/robust/robustlib.h
- Timestamp:
- 04/20/11 20:11:40 (13 years ago)
- Files:
-
- 1 modified
Legend:
- Unmodified
- Added
- Removed
-
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);