Changeset 1395 for applications/robust

Show
Ignore:
Timestamp:
09/21/11 20:35:02 (13 years ago)
Author:
sindj
Message:

Dodelavani toho zatracenyho samplingu. JS

Location:
applications/robust
Files:
3 modified

Legend:

Unmodified
Added
Removed

  • applications/robust/main.cpp

    r1393 r1395  
    123123                        if(has_constant) 
    124124                        { 
    125                                 my_rarx = new RARX(ar_components.size()+1,window_size,true,pow(apriorno,2.0),pow(apriorno,2.0),ar_components.size()+4); 
     125                                my_rarx = new RARX(ar_components.size()+1,window_size,true,sqrt(2*apriorno),sqrt(2*apriorno),ar_components.size()+4); 
    126126                                my_arx  = NULL; 
    127127                        } 
    128128                else 
    129129                        { 
    130                                 my_rarx = new RARX(ar_components.size(),window_size,false,pow(apriorno,2.0),pow(apriorno,2.0),ar_components.size()+3); 
     130                                my_rarx = new RARX(ar_components.size(),window_size,false,sqrt(2*apriorno),sqrt(2*apriorno),ar_components.size()+3); 
    131131                                my_arx  = NULL; 
    132132                        } 
    … …  
    298298        vector<vector<string>> strings; 
    299299 
    300         char* file_string =  "C:\\results\\ar_student_single"; // "C:\\dataADClosePercDiff"; //   
     300        char* file_string =  "C:\\ar_student_single"; // "C:\\dataADClosePercDiff"; //   
    301301 
    302302        char dfstring[80]; 

  • applications/robust/robustlib.cpp

    r1384 r1395  
    210210                                        } 
    211211                                         
    212                                         double factor_multiplier = pow((double)-1,k)/pow((*a_ref).first,sigma_order-k);//pow((double)-1,k+1) 
     212                                        double factor_multiplier = pow((double)-1,k)/pow((*a_ref).first,sigma_order-a_order+k+1);//pow((double)-1,k+1) 
    213213                                                                                 
     214                                        double value = 0; 
    214215                                        if(k!=0) 
    215                                         {        
    216                                                 double value = 0; 
     216                                        {                                                        
    217217                                                ivec control_vec = ivec(); 
    218218                                                control_vec.ins(0,my_emlig->number_of_parameters-a_order+1);                                             
    219                                                 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++) 
     219                                                for(multiset<my_ivec>::iterator combi_ref = this->my_emlig->correction_factors[a_order-k-1].begin();combi_ref!=this->my_emlig->correction_factors[a_order-k-1].upper_bound(my_ivec(control_vec));combi_ref++) 
    220220                                                { 
    221                                                         double bracket_factor = 1/fact(a_order-1-k); 
     221                                                        double bracket_factor = 1/fact(k); 
    222222                                                          
    223223                                                        for(int j = 0;j<(*combi_ref).size();j++) 
    … …  
    229229                                                                                                         
    230230                                                        value += bracket_factor*factor_multiplier*k_multiplier;                                                  
    231                                                 } 
    232  
    233                                                 simplex->insert_gamma(k,value,(*a_ref).first); 
    234                                                 gamma_facs[k] += value; 
     231                                                }                                                
    235232                                        } 
    236233                                        else 
    237234                                        { 
    238                                                 double value = k_multiplier*factor_multiplier/fact(a_order-1);                                           
    239                                                  
    240                                                 simplex->insert_gamma(0,value,(*a_ref).first); 
    241                                                 gamma_facs[0] += value; 
    242                                         }                                                
     235                                                value = k_multiplier*factor_multiplier;                                                  
     236                                        } 
     237 
     238                                        simplex->insert_gamma(k,value,(*a_ref).first); 
     239                                        gamma_facs[k] += value; 
    243240                                }                                
    244241                        } 
    245242 
    246243                        double int_value = 0; 
     244                        double volume = det_jacobian/fact(jacobian.rows()); 
    247245                        int gamma_size = (int)gamma_facs.size(); 
    248246                        if(sigma_order-gamma_size>=0) 
    … …  
    250248                                for(int k = 0;k<gamma_size;k++) 
    251249                                {        
    252                                         int_value += det_jacobian*fact(sigma_order-gamma_size+k)/fact(jacobian.rows())*gamma_facs[k]/pow(2.0,sigma_order);                               
     250                                        int_value += volume*fact(sigma_order+k)*gamma_facs[k]/pow(2.0,((toprow*)this)->condition_order);                                 
    253251                                }        
    254252                        } 
    … …  
    259257                         
    260258                        simplex->probability = int_value; 
     259                        simplex->volume = volume; 
    261260                        //cout << "Probability:" << int_value << endl; 
     261                        //pause(0.1); 
    262262                        return int_value;                
    263263                } 

  • applications/robust/robustlib.h

    r1394 r1395  
    7676 
    7777        double probability; 
     78 
     79        double volume; 
    7880 
    7981        vector<multimap<double,double>> positive_gamma_parameters; 
    … …  
    22372239        { 
    22382240                double sigma = 0; 
    2239                 double pg_sum; 
    2240                 double ng_sum; 
    2241                 do 
    2242                 {                        
    2243                         double rnumber = randu(); 
    2244                          
    2245                                                  
    2246                         double sum_g = 0; 
    2247                         for(int i = 0;i<sampled_simplex->positive_gamma_parameters.size();i++) 
    2248                         { 
    2249                                 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++) 
    2250                                 { 
    2251                                         sum_g += (*g_ref).first/sampled_simplex->positive_gamma_sum; 
    2252  
    2253                                                                  
    2254                                         if(sum_g>rnumber) 
    2255                                         { 
    2256                                                 //itpp::Gamma_RNG* gamma = new itpp::Gamma_RNG(conditions.size()-number_of_parameters,1/(*g_ref).second); 
    2257                                                 //sigma = 1/(*gamma)(); 
    2258                                                  
    2259                                                 // TODO: Tady je chyba, musi zaviset i na parametru order, coz je tady i, ale nevim, jestli +i 
    2260                                                 // nebo -i 
    2261                                                 GamRNG.setup(conditions.size()-number_of_parameters+3,(*g_ref).second); 
    2262                                                                                                                                          
    2263                                                 sigma = 1/GamRNG(); 
    2264  
    2265                                                 // cout << "Sigma mean:   " << (*g_ref).second/(conditions.size()-number_of_parameters-1) << endl;                                                               
    2266                                                 break; 
    2267                                         }                                                        
    2268                                 } 
    2269  
    2270                                 if(sigma!=0) 
    2271                                 { 
     2241                double pg_sum;                                   
     2242                double rnumber = randu();                        
     2243                                         
     2244                double sum_g = 0; 
     2245                for(int i = 0;i<sampled_simplex->positive_gamma_parameters.size();i++) 
     2246                { 
     2247                        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++) 
     2248                        { 
     2249                                sum_g += (*g_ref).first/sampled_simplex->positive_gamma_sum; 
     2250 
     2251                                                         
     2252                                if(sum_g>rnumber) 
     2253                                {                                        
     2254                                        GamRNG.setup(condition_order-number_of_parameters-1+i,1/(*g_ref).second); 
     2255                                                                                                                                 
     2256                                        sigma = 1/GamRNG(); 
     2257                                        // cout << "Sigma mean:   " << (*g_ref).second/(conditions.size()-number_of_parameters-1) << endl;                                                               
     2258 
    22722259                                        break; 
    2273                                 } 
     2260                                }                                                        
    22742261                        } 
    22752262 
    2276                         rnumber = randu(); 
    2277  
    2278                         pg_sum = 0; 
    2279                         for(vector<multimap<double,double>>::iterator v_ref = sampled_simplex->positive_gamma_parameters.begin();v_ref!=sampled_simplex->positive_gamma_parameters.end();v_ref++) 
    2280                         { 
    2281                                 for(multimap<double,double>::iterator pg_ref = (*v_ref).begin();pg_ref!=(*v_ref).end();pg_ref++) 
    2282                                 { 
    2283                                         pg_sum += exp((sampled_simplex->min_beta-(*pg_ref).second)/sigma)*pow((*pg_ref).second/sigma,(int)conditions.size())*(*pg_ref).second/fact(conditions.size())*(*pg_ref).first; 
    2284                                 }                                        
     2263                        if(sigma!=0) 
     2264                        { 
     2265                                break; 
    22852266                        } 
    2286  
    2287                         ng_sum = 0; 
    2288                         for(vector<multimap<double,double>>::iterator v_ref = sampled_simplex->negative_gamma_parameters.begin();v_ref!=sampled_simplex->negative_gamma_parameters.end();v_ref++) 
    2289                         { 
    2290                                 for(multimap<double,double>::iterator ng_ref = (*v_ref).begin();ng_ref!=(*v_ref).end();ng_ref++) 
    2291                                 { 
    2292                                         ng_sum += exp((sampled_simplex->min_beta-(*ng_ref).second)/sigma)*pow((*ng_ref).second/sigma,(int)conditions.size())*(*ng_ref).second/fact(conditions.size())*(*ng_ref).first; 
    2293                                 }                                        
     2267                }        
     2268 
     2269                pg_sum = 0; 
     2270                int i = 0; 
     2271                for(vector<multimap<double,double>>::iterator v_ref = sampled_simplex->positive_gamma_parameters.begin();v_ref!=sampled_simplex->positive_gamma_parameters.end();v_ref++) 
     2272                { 
     2273                        for(multimap<double,double>::iterator pg_ref = (*v_ref).begin();pg_ref!=(*v_ref).end();pg_ref++) 
     2274                        { 
     2275                                pg_sum += exp(-(*pg_ref).second/sigma)*pow((*pg_ref).second/sigma,condition_order-number_of_parameters-1+i)/sigma/fact(condition_order-number_of_parameters-1+i)*(*pg_ref).first;  // pg_sum += exp((sampled_simplex->min_beta-(*pg_ref).second)/sigma)*pow((*pg_ref).second/sigma,(int)conditions.size())*(*pg_ref).second/fact(conditions.size())*(*pg_ref).first; 
    22942276                        } 
    2295                 } 
    2296                 while(pg_sum-ng_sum<0); 
    2297  
    2298                 return pair<double,double>((pg_sum-ng_sum)/pg_sum,sigma); 
     2277 
     2278                        i++; 
     2279                }        
     2280 
     2281                return pair<double,double>(sum_g/pg_sum,sigma); 
    22992282        } 
    23002283 
    … …  
    26562639 
    26572640                        double exponent = extended_coords*condition_and_simplex.first; 
    2658                         double sample_prob = 1/condition_and_simplex.second->probability/pow(probability_and_sigma.second,(int)conditions.size()-number_of_parameters+3)*exp((-1)/probability_and_sigma.second*exponent); 
    2659                         sample_prob *= probability_and_sigma.first; 
     2641                        double sample_prob = 1*condition_and_simplex.second->volume/condition_and_simplex.second->probability/pow(2*probability_and_sigma.second,condition_order)*exp(-exponent/probability_and_sigma.second)*probability_and_sigma.first;                       
    26602642 
    26612643                        sample_coords.ins(sample_coords.size(),probability_and_sigma.second); 
    … …  
    26632645                        samples.ins_col(0,sample_coords); 
    26642646                        probabilities.ins(0,sample_prob); 
     2647 
     2648                        cout << "C:" << sample_coords << "   p:" << sample_prob << endl; 
     2649                        pause(1 ); 
    26652650                } 
    26662651