| 2281 | | return pair<double,double>(sum_g/pg_sum,sigma); |
| 2282 | | } |
| 2283 | | |
| 2284 | | mat sample_mat(int n) |
| 2285 | | { |
| 2286 | | |
| 2287 | | /// \TODO tady je to spatne, tady nesmi byt conditions.size(), viz RARX.bayes() |
| 2288 | | if(conditions.size()-2-number_of_parameters>=0) |
| 2289 | | { |
| 2290 | | mat sample_mat; |
| 2291 | | map<double,toprow*> ordered_toprows; |
| 2292 | | double sum_a = 0; |
| 2293 | | |
| 2294 | | //cout << "Likelihoods of toprows:" << endl; |
| 2295 | | |
| 2296 | | for(polyhedron* top_ref = statistic.rows[number_of_parameters];top_ref!=statistic.end_poly;top_ref=top_ref->next_poly) |
| 2297 | | { |
| 2298 | | toprow* current_top = (toprow*)top_ref; |
| 2299 | | |
| 2300 | | sum_a+=current_top->probability; |
| 2301 | | /* |
| 2302 | | cout << current_top->probability << " "; |
| 2303 | | |
| 2304 | | for(set<vertex*>::iterator vert_ref = (*top_ref).vertices.begin();vert_ref!=(*top_ref).vertices.end();vert_ref++) |
| 2305 | | { |
| 2306 | | cout << round(100*(*vert_ref)->get_coordinates())/100 << " ; "; |
| 2307 | | } |
| 2308 | | */ |
| 2309 | | |
| 2310 | | // cout << endl; |
| 2311 | | ordered_toprows.insert(pair<double,toprow*>(sum_a,current_top)); |
| 2312 | | } |
| 2313 | | |
| 2314 | | // cout << "Sum N: " << normalization_factor << endl; |
| 2315 | | |
| 2316 | | while(sample_mat.cols()<n) |
| 2317 | | { |
| 2318 | | //// cout << "*************************************" << endl; |
| 2319 | | |
| 2320 | | |
| 2321 | | |
| 2322 | | double rnumber = randu()*sum_a; |
| 2323 | | |
| 2324 | | // cout << "RND:" << rnumber << endl; |
| 2325 | | |
| 2326 | | // This could be more efficient (log n), but map::upper_bound() doesn't let me dereference returned iterator |
| 2327 | | int toprow_count = 0; |
| 2328 | | toprow* sampled_toprow; |
| 2329 | | for(map<double,toprow*>::iterator top_ref = ordered_toprows.begin();top_ref!=ordered_toprows.end();top_ref++) |
| 2330 | | { |
| 2331 | | // cout << "CDF:"<< (*top_ref).first << endl; |
| 2332 | | toprow_count++; |
| 2333 | | |
| 2334 | | if((*top_ref).first >= rnumber) |
| 2335 | | { |
| 2336 | | sampled_toprow = (*top_ref).second; |
| 2337 | | break; |
| 2338 | | } |
| 2339 | | } |
| 2340 | | |
| 2341 | | //// cout << "Toprow/Count: " << toprow_count << "/" << ordered_toprows.size() << endl; |
| 2342 | | // cout << &sampled_toprow << ";"; |
| 2343 | | |
| 2344 | | rnumber = randu(); |
| 2345 | | |
| 2346 | | set<simplex*>::iterator s_ref; |
| 2347 | | double sum_b = 0; |
| 2348 | | int simplex_count = 0; |
| 2349 | | for(s_ref = sampled_toprow->triangulation.begin();s_ref!=sampled_toprow->triangulation.end();s_ref++) |
| 2350 | | { |
| 2351 | | simplex_count++; |
| 2352 | | |
| 2353 | | sum_b += (*s_ref)->probability; |
| 2354 | | |
| 2355 | | if(sum_b/sampled_toprow->probability >= rnumber) |
| 2356 | | break; |
| 2357 | | } |
| 2358 | | |
| 2359 | | //// cout << "Simplex/Count: " << simplex_count << "/" << sampled_toprow->triangulation.size() << endl; |
| 2360 | | //// cout << "Simplex factor: " << (*s_ref)->probability << endl; |
| 2361 | | //// cout << "Toprow factor: " << sampled_toprow->probability << endl; |
| 2362 | | //// cout << "Emlig factor: " << normalization_factor << endl; |
| 2363 | | // cout << &(*tri_ref) << endl; |
| 2364 | | |
| 2365 | | int number_of_runs = 0; |
| 2366 | | bool have_sigma = false; |
| 2367 | | double sigma = 0; |
| 2368 | | do |
| 2369 | | { |
| 2370 | | rnumber = randu(); |
| 2371 | | |
| 2372 | | double sum_g = 0; |
| 2373 | | for(int i = 0;i<(*s_ref)->positive_gamma_parameters.size();i++) |
| 2374 | | { |
| 2375 | | for(multimap<double,double>::iterator g_ref = (*s_ref)->positive_gamma_parameters[i].begin();g_ref != (*s_ref)->positive_gamma_parameters[i].end();g_ref++) |
| 2376 | | { |
| 2377 | | sum_g += (*g_ref).first/(*s_ref)->positive_gamma_sum; |
| 2378 | | |
| 2379 | | |
| 2380 | | if(sum_g>rnumber) |
| 2381 | | { |
| 2382 | | //itpp::Gamma_RNG* gamma = new itpp::Gamma_RNG(conditions.size()-number_of_parameters,1/(*g_ref).second); |
| 2383 | | //sigma = 1/(*gamma)(); |
| 2384 | | |
| 2385 | | GamRNG.setup(conditions.size()-number_of_parameters,(*g_ref).second); |
| 2386 | | |
| 2387 | | sigma = 1/GamRNG(); |
| 2388 | | |
| 2389 | | // cout << "Sigma mean: " << (*g_ref).second/(conditions.size()-number_of_parameters-1) << endl; |
| 2390 | | break; |
| 2391 | | } |
| 2392 | | } |
| 2393 | | |
| 2394 | | if(sigma!=0) |
| 2395 | | { |
| 2396 | | break; |
| 2397 | | } |
| 2398 | | } |
| 2399 | | |
| 2400 | | rnumber = randu(); |
| 2401 | | |
| 2402 | | double pg_sum = 0; |
| 2403 | | for(vector<multimap<double,double>>::iterator v_ref = (*s_ref)->positive_gamma_parameters.begin();v_ref!=(*s_ref)->positive_gamma_parameters.end();v_ref++) |
| 2404 | | { |
| 2405 | | for(multimap<double,double>::iterator pg_ref = (*v_ref).begin();pg_ref!=(*v_ref).end();pg_ref++) |
| 2406 | | { |
| 2407 | | pg_sum += exp(((*s_ref)->min_beta-(*pg_ref).second)/sigma)*pow((*pg_ref).second/sigma,(int)conditions.size()-number_of_parameters-1)*(*pg_ref).second/fact(conditions.size()-number_of_parameters-1)*(*pg_ref).first; |
| 2408 | | } |
| 2409 | | } |
| 2410 | | |
| 2411 | | double ng_sum = 0; |
| 2412 | | for(vector<multimap<double,double>>::iterator v_ref = (*s_ref)->negative_gamma_parameters.begin();v_ref!=(*s_ref)->negative_gamma_parameters.end();v_ref++) |
| 2413 | | { |
| 2414 | | for(multimap<double,double>::iterator ng_ref = (*v_ref).begin();ng_ref!=(*v_ref).end();ng_ref++) |
| 2415 | | { |
| 2416 | | ng_sum += exp(((*s_ref)->min_beta-(*ng_ref).second)/sigma)*pow((*ng_ref).second/sigma,(int)conditions.size()-number_of_parameters-1)*(*ng_ref).second/fact(conditions.size()-number_of_parameters-1)*(*ng_ref).first; |
| 2417 | | } |
| 2418 | | } |
| 2419 | | |
| 2420 | | if((pg_sum-ng_sum)/pg_sum>rnumber) |
| 2421 | | { |
| 2422 | | have_sigma = true; |
| 2423 | | } |
| 2424 | | |
| 2425 | | number_of_runs++; |
| 2426 | | } |
| 2427 | | while(!have_sigma); |
| 2428 | | |
| 2429 | | //// cout << "Sigma: " << sigma << endl; |
| 2430 | | //// cout << "Nr. of sigma runs: " << number_of_runs << endl; |
| 2431 | | |
| 2432 | | int dimension = (*s_ref)->vertices.size()-1; |
| 2433 | | |
| 2434 | | mat jacobian(dimension,dimension); |
| 2435 | | vec gradient = sampled_toprow->condition_sum.right(dimension); |
| 2436 | | |
| 2437 | | vertex* base_vert = *(*s_ref)->vertices.begin(); |
| 2438 | | |
| 2439 | | //// cout << "Base vertex coords(should be close to est. param.): " << base_vert->get_coordinates() << endl; |
| 2440 | | |
| 2441 | | int row_count = 0; |
| 2442 | | |
| 2443 | | for(set<vertex*>::iterator vert_ref = ++(*s_ref)->vertices.begin();vert_ref!=(*s_ref)->vertices.end();vert_ref++) |
| 2444 | | { |
| 2445 | | vec current_coords = (*vert_ref)->get_coordinates(); |
| 2446 | | |
| 2447 | | //// cout << "Coords of vertex[" << row_count << "]: " << current_coords << endl; |
| 2448 | | |
| 2449 | | vec relative_coords = current_coords-base_vert->get_coordinates(); |
| 2450 | | |
| 2451 | | jacobian.set_row(row_count,relative_coords); |
| 2452 | | |
| 2453 | | row_count++; |
| 2454 | | } |
| 2455 | | |
| 2456 | | //// cout << "Jacobian: " << jacobian << endl; |
| 2457 | | |
| 2458 | | //// cout << "Gradient before trafo:" << gradient << endl; |
| 2459 | | |
| 2460 | | gradient = jacobian*gradient; |
| 2461 | | |
| 2462 | | //// cout << "Gradient after trafo:" << gradient << endl; |
| 2463 | | |
| 2464 | | // vec normal_gradient = gradient/sqrt(gradient*gradient); |
| 2465 | | // cout << gradient << endl; |
| 2466 | | // cout << normal_gradient << endl; |
| 2467 | | // cout << sqrt(gradient*gradient) << endl; |
| 2468 | | |
| 2469 | | mat rotation_matrix = eye(dimension); |
| 2470 | | |
| 2471 | | |
| 2472 | | |
| 2473 | | for(int i = 1;i<dimension;i++) |
| 2474 | | { |
| 2475 | | vec x_axis = zeros(dimension); |
| 2476 | | x_axis.set(0,1); |
| 2477 | | |
| 2478 | | x_axis = rotation_matrix*x_axis; |
| 2479 | | |
| 2480 | | double t = abs(gradient[i]/gradient*x_axis); |
| 2481 | | |
| 2482 | | double sin_theta = sign(gradient[i])*t/sqrt(1+pow(t,2)); |
| 2483 | | double cos_theta = sign(gradient*x_axis)/sqrt(1+pow(t,2)); |
| 2484 | | |
| 2485 | | mat partial_rotation = eye(dimension); |
| 2486 | | |
| 2487 | | partial_rotation.set(0,0,cos_theta); |
| 2488 | | partial_rotation.set(i,i,cos_theta); |
| 2489 | | |
| 2490 | | partial_rotation.set(0,i,sin_theta); |
| 2491 | | partial_rotation.set(i,0,-sin_theta); |
| 2492 | | |
| 2493 | | rotation_matrix = rotation_matrix*partial_rotation; |
| 2494 | | |
| 2495 | | } |
| 2496 | | |
| 2497 | | // cout << rotation_matrix << endl; |
| 2498 | | |
| 2499 | | mat extended_rotation = rotation_matrix; |
| 2500 | | extended_rotation.ins_col(0,zeros(extended_rotation.rows())); |
| 2501 | | |
| 2502 | | //// cout << "Extended rotation: " << extended_rotation << endl; |
| 2503 | | |
| 2504 | | vec minima = itpp::min(extended_rotation,2); |
| 2505 | | vec maxima = itpp::max(extended_rotation,2); |
| 2506 | | |
| 2507 | | //// cout << "Minima: " << minima << endl; |
| 2508 | | //// cout << "Maxima: " << maxima << endl; |
| 2509 | | |
| 2510 | | vec sample_coordinates; |
| 2511 | | bool is_inside = true; |
| 2512 | | |
| 2513 | | vec new_sample; |
| 2514 | | sample_coordinates = new_sample; |
| 2515 | | |
| 2516 | | for(int j = 0;j<number_of_parameters;j++) |
| 2517 | | { |
| 2518 | | rnumber = randu(); |
| 2519 | | |
| 2520 | | double coordinate; |
| 2521 | | |
| 2522 | | if(j==0) |
| 2523 | | { |
| 2524 | | vec new_gradient = rotation_matrix*gradient; |
| 2525 | | |
| 2526 | | //// cout << "New gradient(should have only first component nonzero):" << new_gradient << endl; |
| 2527 | | |
| 2528 | | // cout << "Max: " << maxima[0] << " Min: " << minima[0] << " Grad:" << new_gradient[0] << endl; |
| 2529 | | |
| 2530 | | double log_bracket = 1-rnumber*(1-exp(new_gradient[0]/sigma*(minima[0]-maxima[0]))); |
| 2531 | | |
| 2532 | | coordinate = minima[0]-sigma/new_gradient[0]*log(log_bracket); |
| 2533 | | } |
| 2534 | | else |
| 2535 | | { |
| 2536 | | coordinate = minima[j]+rnumber*(maxima[j]-minima[j]); |
| 2537 | | } |
| 2538 | | |
| 2539 | | sample_coordinates.ins(j,coordinate); |
| 2540 | | } |
| 2541 | | |
| 2542 | | //// cout << "Sampled coordinates(gradient direction): " << sample_coordinates << endl; |
| 2543 | | |
| 2544 | | sample_coordinates = rotation_matrix.T()*sample_coordinates; |
| 2545 | | |
| 2546 | | //// cout << "Sampled coordinates(backrotated direction):" << sample_coordinates << endl; |
| 2547 | | |
| 2548 | | |
| 2549 | | for(int j = 0;j<sample_coordinates.size();j++) |
| 2550 | | { |
| 2551 | | if(sample_coordinates[j]<0) |
| 2552 | | { |
| 2553 | | is_inside = false; |
| 2554 | | } |
| 2555 | | } |
| 2556 | | |
| 2557 | | double above_criterion = ones(sample_coordinates.size())*sample_coordinates; |
| 2558 | | |
| 2559 | | if(above_criterion>1) |
| 2560 | | { |
| 2561 | | is_inside = false; |
| 2562 | | } |
| 2563 | | |
| 2564 | | if(is_inside) |
| 2565 | | { |
| 2566 | | sample_coordinates = jacobian.T()*sample_coordinates+(*base_vert).get_coordinates(); |
| 2567 | | |
| 2568 | | sample_coordinates.ins(0,sigma); |
| 2569 | | |
| 2570 | | //// cout << "Sampled coordinates(parameter space):" << sample_coordinates << endl; |
| 2571 | | |
| 2572 | | sample_mat.ins_col(0,sample_coordinates); |
| 2573 | | |
| 2574 | | // cout << sample_mat.cols() << ","; |
| 2575 | | } |
| 2576 | | |
| 2577 | | // cout << sampled_toprow->condition_sum.right(sampled_toprow->condition_sum.size()-1)*min_grad->get_coordinates()-sampled_toprow->condition_sum[0] << endl; |
| 2578 | | // cout << sampled_toprow->condition_sum.right(sampled_toprow->condition_sum.size()-1)*max_grad->get_coordinates()-sampled_toprow->condition_sum[0] << endl; |
| 2579 | | |
| 2580 | | |
| 2581 | | } |
| 2582 | | |
| 2583 | | cout << endl; |
| 2584 | | return sample_mat; |
| 2585 | | } |
| 2586 | | else |
| 2587 | | { |
| 2588 | | throw new exception("You are trying to sample from density that is not determined (parameters can't be integrated out)!"); |
| | 2282 | return pair<double,double>(sampled_simplex->positive_gamma_sum/pg_sum,sigma); |
| | 2283 | } |
| | 2284 | |
| | 2285 | pair<vec,mat> sample(int n, bool rejection) |
| | 2286 | { |
| | 2287 | vec probabilities; |
| | 2288 | mat samples; |
