|
Other R functions
I have developed many other functions that, although being useful, are too specific to be part of an R Package published on CRAN. I have coded most of the functions I have in mind during the Fall semester of 2016, while I was TA'ing Probability and Stochastic Processes (PSP) at Barcelona Tech. Here I publish the code for one of them that the PSP professor included in the course materials for the students; I think it's particularly interesting and fun: it shows the expected number of coin tosses needed in order to obtain a certain Heads-Tails pattern. Enjoy!
################################################# # Expectations of Geometric R.V's with Heads-Tails patterns # # (C) Albert Dorador, November 2016 # ################################################# We'll perform a Monte Carlo simulation and, by the Law of Large Numbers, we will obtain a very good approximation to the true expected value of the number of coin tosses needed to achieve the specified pattern. hat.Exp.Geom = function(objective.seq, reps = 1e5){ toss.vector = integer(reps) #integer is 4 bytes, numeric is 8. len.seq = length(objective.seq) for (i in 1:reps){ seq_1 = rbinom(n=len.seq,size=1,prob=0.5) tosses= len.seq #at least we'll need n tosses if (all(seq_1 == objective.seq)){ #maybe we are lucky success = T #we forbid entering the loop }else{ success = F seq_i = seq_1 #so that later in the while loop we call seq_i, instead of seq_1 while(!success){ xi = rbinom(n=1,size=1,prob=0.5) #1 Bernoulli experiment tosses = tosses + 1 seq_i = seq_i[-1] #objective.seq has length n, clearly a sequence of n+1 cannot be correct seq_i = append(seq_i, xi) if (all(seq_i == objective.seq)){ success = T #only in case all characters coincide we give permission to exit the while loop } } } toss.vector[i] = tosses } hat.Expect = mean(toss.vector) #Law of Large Numbers return(hat.Expect) } #Let's test set.seed(123) hat.Exp.Geom(c(1, 0)) #1 is Heads, 0 is Tails #4.00101, correct hat.Exp.Geom(c(1, 1)) #6.01517, correct hat.Exp.Geom(c(1, 0, 1)) # 9.96092, correct hat.Exp.Geom(c(1, 0, 0)) # 7.9962, correct # Now let's try something really hard, which would take us hours # to compute analytically (and with some chance of error): set.seed(321) hat.Exp.Geom(c(1, 0, 0, 1, 0, 1, 1, 1, 0)) # 515.9843 ... probably correct (I haven't calculated this one analytically!) |