Changeset 1379 for applications/robust
- Timestamp:
- 07/15/11 20:26:08 (13 years ago)
- Location:
- applications/robust
- Files:
-
- 3 modified
-
main.cpp (modified) (9 diffs)
-
robustlib.cpp (modified) (6 diffs)
-
robustlib.h (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
applications/robust/main.cpp
r1376 r1379 28 28 29 29 const int max_model_order = 2; 30 const double apriorno=0.0 2;30 const double apriorno=0.005; 31 31 32 32 /* … … 74 74 // problematic. 75 75 RARX* my_rarx; //vzmenovane parametre pre triedu model 76 ARX * my_arx;76 ARXwin* my_arx; 77 77 78 78 bool has_constant; … … 110 110 { 111 111 my_rarx = NULL; 112 my_arx = new ARX ();112 my_arx = new ARXwin(); 113 113 mat V0; 114 114 … … 117 117 V0 = apriorno * eye(ar_components.size()+2); //aj tu konst 118 118 V0(0,0) = 1; 119 my_arx->set_constant(true); 120 119 my_arx->set_constant(true); 121 120 } 122 121 else … … 129 128 } 130 129 131 my_arx->set_statistics(1, V0 );132 //my_arx->set_parameters(window_size);130 my_arx->set_statistics(1, V0, V0.rows()+1); 131 my_arx->set_parameters(window_size); 133 132 my_arx->validate(); 134 133 } … … 445 444 vector<vector<string>> strings; 446 445 447 char* file_string = "c:\\ar_ cauchy_single";446 char* file_string = "c:\\ar_normal_single"; // "c:\\dataTYClosePercDiff"; // 448 447 449 448 char dfstring[80]; … … 485 484 for(set<set<pair<int,int>>>::iterator model_type = model_types.begin();model_type!=model_types.end();model_type++) 486 485 {// prechadza rozne typy kanalov, a poctu regresorov 487 for(int window_size = 1 00;window_size < 101;window_size++)486 for(int window_size = 15;window_size < 16;window_size++) 488 487 { 489 488 models.push_back(new model((*model_type),true,true,window_size,0,&data_matrix)); // to su len konstruktory, len inicializujeme rozne typy … … 493 492 } 494 493 495 set<pair<int,int>> empty_list;496 models.push_back(new model(empty_list,false,true,100,0,&data_matrix));494 //set<pair<int,int>> empty_list; 495 //models.push_back(new model(empty_list,false,true,100,0,&data_matrix)); 497 496 } 498 497 … … 508 507 {//posuvam s apo models, co je pole modelov urobene o cyklus vyssie. Teda som v case time a robim to tam pre vsetky typy modelov, kombinace regresorov 509 508 (*model_ref)->data_update(time); //pozret sa preco je toto tu nutne 509 510 cout << "Updated." << endl; 510 511 //if (time = max_model_order) nazvy.push_back(models.model_ref]);// ako by som mohol dostat nazov modelu? 511 512 -
applications/robust/robustlib.cpp
r1376 r1379 111 111 emlig* current_emlig; 112 112 simplex->clear_gammas(); 113 simplex->probability = 0; 113 114 114 115 if(this->my_emlig!=NULL) … … 128 129 129 130 map<double,int> as; 130 vertex* base_vertex = simplex->vertices.begin(); 131 131 vertex* base_vertex = *simplex->vertices.begin(); 132 132 133 int row_count = 0; 133 134 for(set<vertex*>::iterator vert_ref = simplex->vertices.begin(); vert_ref!=simplex->vertices.end();vert_ref++) 134 135 { … … 136 137 { 137 138 vec relative_coords = (*vert_ref)->get_coordinates()-base_vertex->get_coordinates(); 138 jacobian.set_row(row_count,relative_coords); 139 } 140 141 double a = -((*vert_ref)->get_coordinates().ins(0,-1)*cur_condition); 139 jacobian.set_row(row_count-1,relative_coords); 140 } 141 142 vec extended_coords = (*vert_ref)->get_coordinates(); 143 extended_coords.ins(0,-1); 144 145 double a = extended_coords*as_toprow->condition_sum; 142 146 if(a<current_emlig->min_ll) 143 147 { … … 155 159 { 156 160 (*returned.first).second++; 157 } 158 } 159 161 } 162 163 row_count++; 164 } 160 165 161 166 /* 162 cout << "a_0: " << a_0 << " ";163 167 int as_count = 1; 164 168 for(map<double,int>::iterator as_ref = as.begin();as_ref!=as.end();as_ref++) … … 167 171 as_count++; 168 172 } 169 */ 170 171 double int_value = 0; 172 173 */ 174 175 vector<double> gamma_facs; 173 176 // cout << jacobian << endl; 174 175 double det_jacobian = abs(det(jacobian)); 176 double correction_term; 177 for(map<double,int>::iterator a_ref = as.begin();as_ref!=as.end();as_ref++) 178 { 179 double multiplier = det_jacobian/fact(jacobian.rows()); 177 double det_jacobian = abs(det(jacobian)); 178 for(map<double,int>::iterator a_ref = as.begin();a_ref!=as.end();a_ref++) 179 { 180 double k_multiplier = 1; 180 181 181 182 int a_order = (*a_ref).second; 182 current_emlig->set_correction_factors(a_order); 183 184 vector<double> factors; 185 int number_of_factors = 0; 183 current_emlig->set_correction_factors(a_order-1); 184 185 vector<double> factors; 186 186 for(map<double,int>::iterator a2_ref = as.begin();a2_ref!=as.end();a2_ref++) 187 187 { … … 190 190 for(int k = 0;k<(*a2_ref).second;k++) 191 191 { 192 factors.push_back((*a _ref).first-(*a2_ref).first);192 factors.push_back((*a2_ref).first-(*a_ref).first); 193 193 } 194 194 195 multiplier /= pow((*a_ref).first-(*a2_ref).first,(*a2_ref).second); 196 number_of_factors += (*a2_ref).second; 195 k_multiplier /= pow((*a2_ref).first-(*a_ref).first,(*a2_ref).second); 197 196 } 197 } 198 199 double bracket_combination = 0; 200 201 for(int k = 0;k < a_order;k++) 202 { 203 if(k==gamma_facs.size()) 204 { 205 gamma_facs.push_back(0.0); 206 } 207 208 double factor_multiplier = pow((double)-1,k)/pow((*a_ref).first,sigma_order-k);//pow((double)-1,k+1) 209 210 if(k!=0) 211 { 212 ivec control_vec = ivec(); 213 control_vec.ins(0,my_emlig->number_of_parameters-a_order+1); 214 for(multiset<my_ivec>::iterator combi_ref = this->my_emlig->correction_factors[k-1].begin();combi_ref!=this->my_emlig->correction_factors[k-1].upper_bound(my_ivec(control_vec));combi_ref++) 215 { 216 double bracket_factor = 1/fact(a_order-1-k); 217 218 for(int j = 0;j<(*combi_ref).size();j++) 219 { 220 //cout << "Factor vector:" << (*combi_ref) << endl; 221 222 bracket_factor /= factors[(*combi_ref)[j]-1]; 223 } 224 225 double value = bracket_factor*factor_multiplier*k_multiplier; 226 227 simplex->insert_gamma(k,value,(*a_ref).first); 228 gamma_facs[k] += value; 229 } 230 } 231 else 232 { 233 double value = k_multiplier*factor_multiplier/fact(a_order-1); 234 235 simplex->insert_gamma(0,value,(*a_ref).first); 236 gamma_facs[k] += value; 237 } 238 } 239 } 240 241 double int_value = 0; 242 int gamma_size = (int)gamma_facs.size(); 243 if(sigma_order-gamma_size>=0) 244 { 245 for(int k = 0;k<gamma_size;k++) 246 { 247 int_value += det_jacobian*fact(sigma_order-gamma_size+k)/fact(jacobian.rows())*gamma_facs[k]/pow(2.0,sigma_order); 198 248 } 199 200 double cur_as_ref = (*a_ref).first; 201 202 double gamma_multiplier = -cur_as_ref-a_0; 203 204 double bracket = fact(as1_order-1)/pow(gamma_multiplier,sigma_order); 205 206 simplex->insert_gamma(0,bracket*multiplier,gamma_multiplier); 207 208 // bracket *= fact(sigma_order); 209 210 for(int k = 0;k < as1_order-1;k++) 211 { 212 double local_bracket = 0; 213 double bracket_factor = 1/fact(as1_order-1-k)/pow(gamma_multiplier,sigma_order-k-1);//pow((double)-1,k+1) 214 215 ivec control_vec = ivec(); 216 control_vec.ins(0,my_emlig->number_of_parameters-as1_order+1); 217 218 for(multiset<my_ivec>::iterator combi_ref = this->my_emlig->correction_factors[k].begin();combi_ref!=this->my_emlig->correction_factors[k].upper_bound(my_ivec(control_vec));combi_ref++) 219 { 220 double bracket_combination = 1; 221 for(int j = 0;j<(*combi_ref).size();j++) 222 { 223 //cout << "Factor vector:" << (*combi_ref) << endl; 224 225 bracket_combination /= factors[(*combi_ref)[j]-1]; 226 } 227 228 local_bracket += bracket_factor*bracket_combination; 229 } 230 231 simplex->insert_gamma(k+1,local_bracket*multiplier,gamma_multiplier); 232 233 int division_factor = 1; 234 for(int s = 0;s<k;s++) 235 { 236 division_factor *= k-s; 237 } 238 239 bracket += local_bracket/division_factor; //local_bracket*fact(sigma_order-k); 240 } 241 242 int_value += multiplier*bracket; 243 244 245 246 } 247 248 double correction_term_base = correction_term/pow(-a_0,sigma_order); 249 250 simplex->insert_gamma(0,correction_term_base,-a_0); 251 252 correction_term = correction_term_base;//fact(sigma_order)*correction_term_base; 253 254 //cout << c << int_value << endl; 255 256 int_value += correction_term; 257 249 } 250 else 251 { 252 int_value = -1; 253 } 254 255 simplex->probability = int_value; 258 256 //cout << "Probability:" << int_value << endl; 259 //pause(0.100); 260 261 simplex->probability = int_value; 262 263 return int_value; 264 257 return int_value; 265 258 } 266 259 else -
applications/robust/robustlib.h
r1376 r1379 1077 1077 this->number_of_parameters = number_of_parameters; 1078 1078 1079 condition_order = number_of_parameters+ 2;1079 condition_order = number_of_parameters+3; 1080 1080 1081 1081 create_statistic(number_of_parameters, soft_prior_parameter); … … 1097 1097 } 1098 1098 1099 log_nc = log(normalization_factor) + logfact(condition_order-number_of_parameters-2)-(condition_order-number_of_parameters-2)*log(2.0);1099 log_nc = log(normalization_factor); 1100 1100 1101 1101 /* … … 1105 1105 */ 1106 1106 1107 1107 cout << "Prior constructed." << endl; 1108 1108 } 1109 1109 … … 1905 1905 // cout << "Normalization factor: " << normalization_factor << endl; 1906 1906 1907 log_nc = log(normalization_factor) + logfact(condition_order-number_of_parameters-2)-(condition_order-number_of_parameters-2)*log(2.0); 1908 1907 log_nc = log(normalization_factor); // + logfact(condition_order-number_of_parameters-2)-(condition_order-number_of_parameters-2)*log(2.0); 1908 1909 /* 1909 1910 if(condition_order == 20) 1910 1911 step_me(88); 1912 */ 1911 1913 1912 1914 //cout << "Factorial factor: " << condition_order-number_of_parameters-2 << endl; … … 2818 2820 this->has_constant = has_constant; 2819 2821 2820 posterior = new emlig(number_of_parameters,0. 001);2822 posterior = new emlig(number_of_parameters,0.1); 2821 2823 2822 2824 this->window_size = window_size;
