Newer
Older
Lecture_repo / Lectures_my / EMPP / 2017 / Lecture4 / Exercises / Markov / mark.cc
@Marcin CHrzaszcz Marcin CHrzaszcz on 3 Oct 2017 1 KB started the lectores 2017
  1. #include <cmath>
  2. #include<iostream>
  3. #include <stdio.h>
  4. #include <stdlib.h>
  5. #include <TRandom1.h>
  6.  
  7.  
  8. using namespace std;
  9.  
  10. int main(int argc, char **argv)
  11. {
  12. double h[3][4]={ {0.4, 0.2, 0.3,0.1},{0.1,0.4,0.3,0.2},{0.5,0.3,0.1,0.1}};
  13. double a[3]={1.5,-1.0,0.7};
  14. int j = atoi(argv[1]);
  15. int original =j;
  16. int n_of_samples = atoi(argv[2]);
  17. TRandom1 *rand = new TRandom1(2131);
  18. double result =0;
  19. double result2=0;
  20. for(int i =0;i< n_of_samples;++i)
  21. {
  22. int last;
  23. j=original;
  24. while(1)
  25. {
  26. last =j;
  27. double tmp =rand->Rndm();
  28. if(tmp<h[j-1][0]) j=0;
  29. else if(tmp > h[j-1][0] && tmp<h[j-1][1]+h[j-1][0]) j=1;
  30. else if(tmp > h[j-1][1]+h[j-1][0] && tmp<h[j-1][2]+h[j-1][1]+h[j-1][0]) j=2;
  31. else if(tmp > h[j-1][2]+h[j-1][1]+h[j-1][0] && tmp<h[j-1][3]+h[j-1][2]+h[j-1][1]+h[j-1][0]) j=3;
  32. else
  33. {
  34. cout<<"rand= "<<tmp<<endl;
  35. cout<<h[j-1][3]+h[j-1][2]+h[j-1][1]+h[j-1][0]<<endl;
  36. }
  37. if(j==0) break;
  38. }
  39. result+=a[last-1]/(double)(h[last-1][0]);
  40. //errors
  41. }
  42.  
  43. cout<<result/ n_of_samples<<endl;
  44.  
  45.  
  46. return 1;
  47. }