Changeset 1379 for applications/robust

Show
Ignore:
Timestamp:
07/15/11 20:26:08 (13 years ago)
Author:
sindj
Message:

Prepracovani integrace v robustlib.cpp, zjednoduseni, zesymetrizovani. Nevim jestli funguje, nutno overit a vyzkouset testovani hypotez a integraci na jednicku (u normalizacnich faktoru).JS

Location:
applications/robust
Files:
3 modified

Legend:

Unmodified
Added
Removed
  • applications/robust/main.cpp

    r1376 r1379  
    2828 
    2929const int max_model_order = 2; 
    30 const double apriorno=0.02; 
     30const double apriorno=0.005; 
    3131 
    3232/* 
     
    7474        // problematic. 
    7575        RARX* my_rarx; //vzmenovane parametre pre triedu model 
    76         ARX* my_arx; 
     76        ARXwin* my_arx; 
    7777 
    7878        bool has_constant; 
     
    110110                { 
    111111                        my_rarx = NULL; 
    112                         my_arx  = new ARX(); 
     112                        my_arx  = new ARXwin(); 
    113113                        mat V0; 
    114114 
     
    117117                                V0  = apriorno * eye(ar_components.size()+2); //aj tu konst 
    118118                                V0(0,0) = 1; 
    119                                 my_arx->set_constant(true);      
    120                                  
     119                                my_arx->set_constant(true);                              
    121120                        } 
    122121                        else 
     
    129128                        } 
    130129 
    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); 
    133132                        my_arx->validate(); 
    134133                } 
     
    445444        vector<vector<string>> strings; 
    446445 
    447         char* file_string = "c:\\ar_cauchy_single";  
     446        char* file_string = "c:\\ar_normal_single"; // "c:\\dataTYClosePercDiff"; //  
    448447 
    449448        char dfstring[80]; 
     
    485484        for(set<set<pair<int,int>>>::iterator model_type = model_types.begin();model_type!=model_types.end();model_type++) 
    486485        {// prechadza rozne typy kanalov, a poctu regresorov 
    487                 for(int window_size = 100;window_size < 101;window_size++) 
     486                for(int window_size = 15;window_size < 16;window_size++) 
    488487                { 
    489488                        models.push_back(new model((*model_type),true,true,window_size,0,&data_matrix));   // to su len konstruktory, len inicializujeme rozne typy 
     
    493492                } 
    494493 
    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)); 
    497496        } 
    498497         
     
    508507                {//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 
    509508                        (*model_ref)->data_update(time); //pozret sa preco je toto tu nutne 
     509 
     510                        cout << "Updated." << endl; 
    510511                        //if (time = max_model_order) nazvy.push_back(models.model_ref]);// ako by som mohol dostat nazov modelu? 
    511512 
  • applications/robust/robustlib.cpp

    r1376 r1379  
    111111                        emlig* current_emlig; 
    112112                        simplex->clear_gammas(); 
     113                        simplex->probability = 0; 
    113114 
    114115                        if(this->my_emlig!=NULL) 
     
    128129 
    129130                        map<double,int> as; 
    130                         vertex* base_vertex = simplex->vertices.begin(); 
    131                                                          
     131                        vertex* base_vertex = *simplex->vertices.begin();                                                        
    132132                                 
     133                        int row_count = 0; 
    133134                        for(set<vertex*>::iterator vert_ref = simplex->vertices.begin(); vert_ref!=simplex->vertices.end();vert_ref++) 
    134135                        { 
     
    136137                                { 
    137138                                        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; 
    142146                                if(a<current_emlig->min_ll) 
    143147                                { 
     
    155159                                { 
    156160                                        (*returned.first).second++; 
    157                                 }                                                
    158                         } 
    159                          
     161                                } 
     162                                 
     163                                row_count++; 
     164                        }                
    160165                         
    161166                        /* 
    162                         cout << "a_0: " << a_0 << "    "; 
    163167                        int as_count = 1; 
    164168                        for(map<double,int>::iterator as_ref = as.begin();as_ref!=as.end();as_ref++) 
     
    167171                                as_count++; 
    168172                        } 
    169                         */ 
    170  
    171                         double int_value = 0;                    
    172  
     173                        */               
     174 
     175                        vector<double> gamma_facs; 
    173176                        // 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;                         
    180181                                 
    181182                                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;                                                                                          
    186186                                for(map<double,int>::iterator a2_ref = as.begin();a2_ref!=as.end();a2_ref++) 
    187187                                {                                        
     
    190190                                                for(int k = 0;k<(*a2_ref).second;k++) 
    191191                                                { 
    192                                                         factors.push_back((*a_ref).first-(*a2_ref).first); 
     192                                                        factors.push_back((*a2_ref).first-(*a_ref).first); 
    193193                                                } 
    194194 
    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);                                            
    197196                                        }                                                                                
     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);                               
    198248                                }        
    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; 
    258256                        //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;                
    265258                } 
    266259                else 
  • applications/robust/robustlib.h

    r1376 r1379  
    10771077                this->number_of_parameters = number_of_parameters; 
    10781078 
    1079                 condition_order = number_of_parameters+2; 
     1079                condition_order = number_of_parameters+3; 
    10801080                                                 
    10811081                create_statistic(number_of_parameters, soft_prior_parameter); 
     
    10971097                } 
    10981098 
    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); 
    11001100 
    11011101                /* 
     
    11051105                */ 
    11061106                 
    1107  
     1107                cout << "Prior constructed." << endl; 
    11081108        } 
    11091109 
     
    19051905                // cout << "Normalization factor: " << normalization_factor << endl;     
    19061906 
    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                /* 
    19091910                if(condition_order == 20) 
    19101911                                                        step_me(88); 
     1912                                                        */ 
    19111913 
    19121914                //cout << "Factorial factor: " << condition_order-number_of_parameters-2 << endl; 
     
    28182820                this->has_constant = has_constant; 
    28192821                 
    2820                 posterior = new emlig(number_of_parameters,0.001); 
     2822                posterior = new emlig(number_of_parameters,0.1); 
    28212823 
    28222824                this->window_size = window_size;