Advertisement
Mihai_Preda

Untitled

Feb 2nd, 2021
229
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 3.09 KB | None | 0 0
  1. # Marja de eroare pentru comparatii pe numere reale.
  2. EPSILON <- 10e-9
  3.  
  4. CONST_DISP <- 5
  5.  
  6. # De cate ori verificam sa fie PDF pe fiecare interval al functiei date.
  7. TESTS <- 500
  8.  
  9. equalToOne <- function (x)
  10. {
  11. return (x < 1 + EPSILON & x > 1 - EPSILON)
  12. }
  13.  
  14. calcIntegralSum <- function (myFunc, domain)
  15. {
  16. sum <- 0
  17. for (i in domain)
  18. sum <- sum + integrate (Vectorize (myFunc), i[1], i[2])$value
  19. return (sum)
  20. }
  21.  
  22. # Ca o functie sa fie densitate de probabilitate trebuie sa fie pozitiva, iar integrala ei pe R sa fie 1.
  23. checkDensProb <- function (myFunc, domain = list (c (-Inf, Inf)), constNorm = 1)
  24. {
  25. tryCatch (
  26. {
  27. # Calculam suma integralelor, similar cu cerinta 1.
  28. integralSum <- calcIntegralSum (myFunc, domain)
  29.  
  30. # Comparam cu 1 folosindu-ne de EPSILON.
  31. if (equalToOne (integralSum / constNorm))
  32. {
  33. # Pentru fiecare interval al functiei, verificam daca este PDF.
  34. for (i in domain)
  35. {
  36. # Verificam intervalul [st, dr].
  37. st <- max(i[1], -2e9)
  38. dr <- min(i[2], 2e9)
  39.  
  40. length <- dr / CONST_DISP - st / CONST_DISP
  41. nrs <- runif (TESTS, st, dr-length)
  42.  
  43. # Testam TESTS subintervale de lungime egala.
  44. for (mn in nrs)
  45. {
  46. mx = mn + length
  47. # Daca integrala este negativa, functia nu este PDF.
  48. if (integrate (Vectorize (myFunc), mn, mx)$value < 0) return (FALSE)
  49. }
  50. }
  51. return (TRUE)
  52. }
  53. return (FALSE)
  54. },
  55. error = function (e)
  56. {
  57. # Daca integrarea a esuat, returnam false.
  58. return (FALSE)
  59. })
  60. }
  61.  
  62. isPdf <- function(f, suport = list(c(-Inf, Inf)), normCt = 1, TCNT = 1000){
  63.  
  64. tryCatch({
  65.  
  66. integratedVal <- 0
  67.  
  68. for(interval in suport){
  69.  
  70. integratedVal <- integratedVal + integrate(Vectorize(f), lower = interval[1], upper = interval[2])$value
  71. }
  72.  
  73. if((integratedVal / normCt) < 1 + 10e-9 & (integratedVal / normCt) > 1 - 10e-9){
  74.  
  75. cnt <- TCNT / length(suport)
  76.  
  77. for(interval in suport){
  78.  
  79. lowRand <- 0
  80. upperRand <- 0
  81.  
  82. if(interval[1] == -Inf){
  83. lowRand <- -.Machine$integer.max
  84. }
  85. else{
  86. lowRand <- interval[1]
  87. }
  88.  
  89. if(interval[2] == Inf){
  90. upperRand <- .Machine$integer.max
  91. }
  92. else{
  93. upperRand <- interval[2]
  94. }
  95.  
  96. testLen <- upperRand / 10 - lowRand / 10
  97.  
  98. rands <- runif(cnt, lowRand, upperRand - testLen)
  99.  
  100. for(x in c(1:cnt)){
  101.  
  102. lowerBound <- rands[x]
  103. upperBound = lowerBound + testLen
  104.  
  105. if(lowerBound > upperBound){
  106.  
  107. swapAux <- upperBound
  108. upperBound <- lowerBound
  109. lowerBoud <- swapAux
  110. }
  111.  
  112. if(integrate(Vectorize(f), lower = lowerBound, upper = upperBound)$value < 0)
  113. return(FALSE)
  114. }
  115.  
  116. }
  117.  
  118. return(TRUE)
  119. }
  120. else
  121. return(FALSE)
  122.  
  123. },
  124. error = function(err){
  125.  
  126. message(paste("eroare a functiei isPdf: ", err))
  127. return(FALSE)
  128. })
  129.  
  130. }
  131.  
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement