As you can realize, if you know how to solve the forward problem, you can solve the backward problem by brute force

You try several starting points, apply the forward method and see if you find the backward end point

For example, to solve 78/6 you can try

(1:15) * 6

 [1]  6 12 18 24 30 36 42 48 54 60 66 72 78 84 90

Experiment

system that produces a result. We do not know the result before

Example: survey, dice, counting cells, a plant did grow or not, measuring cytotoxicity.

Outcome

result of an experiment. Only one outcome for each experiment.

Population

Set of all possible outcomes. Very big, we assume infinite size.

Sample

Set of outcomes from a few experiments. Much smaller than the population.

Population is a very very big set of things

• All people in Turkey
• All humans living today
• All humans that ever lived
• All living organisms in the past and in the future
• All the infinite possible results of an experiment

You can assume that population size is infinite

For example

⚀ ⚁ ⚂ ⚃ ⚄ ⚅

A, C, G, T

Each experiment gives us an outcome

Several experiments are a sample

They are randomly taken from the population

A Sample is a small part of population, chosen randomly

• small number of elements
• changes every time
  • if we take two samples they will be different

All letters have the same chance of being the first one chosen

sample(LETTERS, size=5)

[1] "X" "R" "T" "A" "K"

But the first letter cannot be the second, o 3rd, 4th, 5th.

Every time we get an outcome, the chances change

In real life, we usually sample from huge sets

• We survey some people from a city, country or planet
• We test some plants among the billions of plants
• We do some experiments among the trillions of possible experiments

In practice, our sample is small compared to the total

Therefore the effect is small and chances do not change

If the population is big and our sample is small, then chances will not change

We can simulate that using a big population (bad idea)

sample(rep(LETTERS, 1e6), size=5)

Don’t do this: It takes too much time and memory.

A possible solution is this:

ans <- rep(NA, 5)
for(i in 1:5) {
    ans[i] <- sample(LETTERS, size=1)
}

Taking one letter each time does not change the chances

Other way of sampling without changing the chances is this

sample(LETTERS, size=5, replace=TRUE)

In this case we replace (put back) each letter we get

replace=TRUE means chances do not change

In other cases different outcomes may have different proportions

For example

(A, C, C, G, G, T)

(A, A, C, C, C, G, G, G, T)

If we know the probabilities, we can tell R to use them to simulate

To simulate (A, C, C, G, G, T) we use

sample(c("A","C","G","T"), size=N, replace=TRUE,
        prob=c(1/6, 1/3, 1/3, 1/6))

For example, when we throw two dice and count points, we have

• zero ways of getting 1,
• one way of getting 2,
• six ways of getting 7,
• one way of getting 12,
• etc.

The technique for finding these numbers is called combinatorics

We will not do combinatorics here

To simulate two dice, we have

sample(1:12, size=N, replace=TRUE,
        prob=c(0,1,2,3,4,5,6,5,4,3,2,1)/36)

If the system we want to simulate is complex, we usually do not know the probabilities

In that case we can decompose the random system

• small parts are easier to simulate
• combine them to make the complex system

We roll two dice. What is the sum 🎲 +🎲 ?

roll2 <- function() {
    dice <- sample(1:6, size=2, replace=TRUE)
    return(sum(dice))
}
roll2()

[1] 7

The output is a random outcome

Try it. What is your result?

One experiment is meaningless

We need to replicate the experiment

ans <- rep(NA, 15)
for(i in 1:15) {
    ans[i] <- roll2()
}

We can use the function replicate(n, expression)

replicate(15, roll2() )

 [1]  4  4  6  2  5  8  4  6  6  6  9  5  5  8 10

table(deney)

deney
 2  3  4  5  6  7  8  9 10 11 12 
 6 10 17 19 37 33 28 20 20  7  3 

barplot(table(deney))

deney <- replicate(200, roll2() )
barplot(table(deney) / 200)

deney <- replicate(1000, roll2() )
barplot(table(deney) / 1000)

deney <- replicate(40000, roll2() )
barplot(table(deney) / 40000)

The result of our experiments give us empirical frequencies

They are close to the theoretical probabilities

When size is bigger, the empirical frequencies are closer and closer to the real probabilities

We know for sure that when size grows we will get the probabilities

But size has to be really big

The main consequence of the Law of Large Numbers is

Samples tell us something about populations

Therefore we can learn about populations if we do experiments

In our course experiment means sample(x, size, replace=TRUE)

freq <- table(deney)
freq["7"]

   7 
6613 

This is a little tricky

It is easy to make a mistake here

Remember that we can use text and numbers as indices.

Here freq is

freq

deney
   2    3    4    5    6    7    8    9   10   11   12 
1093 2210 3411 4490 5659 6613 5577 4427 3257 2147 1116 

What is freq[12]? Be careful

If we want to know how may times we rolled a 7, we can do

sum(deney == 7)

[1] 6613

Remember that:

• Comparisons produce logic vectors
• sum() of a logic vector is the number of TRUE

Here we are in more trouble, since there is no freq["1"]

freq["1"]

<NA> 
  NA 

It is better to ask in a safe way

sum(deney == 1)

[1] 0

What is the GC content of a random gene?

You can see that GC content is a random variable

gene <- sample(c("A","C","T","G"), size,
                replace=TRUE)
gc_gene <- sum(gene=="C" | gene=="G")/length(gene)

size is the gene size. Symbol | means or

The outcome is the gene. The random variable is gc_gene

gc_one_gene <- function(M) {
    gene <- sample(c("A","C","T","G"), size=M, replace=TRUE)
    gc_gene <- sum(gene=="C" | gene=="G")/length(gene)
    return(gc_gene)
}
gc_genes <- replicate(10000, gc_one_gene(900))

hist(gc_genes, col = "grey")

hist(gc_genes, col = "grey", freq=FALSE)

We can have a better simulation if we use the real proportions for ("A","C","T","G")

Let’s evaluate them for Carsonella rudii

library(seqinr)
genome <- read.fasta("AP009180.fna")
table(genome[[1]])

         a          c          g          t 
0.41797046 0.08455988 0.08108379 0.41638587 

gc_one_gene <- function(M, prob) {
    gene <- sample(c("A","C","G","T"), size=M, prob=prob,
                replace=TRUE)
    gc_gene <- sum(gene=="C" | gene=="G")/length(gene)
    return(gc_gene)
}
gc_genes <- replicate(10000, gc_one_gene(900, prob))

hist(gc_genes, col = "grey", freq=FALSE)

• Average
• Number of trials until success
• Number plants that grow in the lab
• Number of people that answer the homework
• Number of people that passes the course

What other random variable do you see?

Depending on how many questions you ask and how many homework you do, you will have different chances of passing this course

To be neutral, let’s try with different probabilities

How many times do you need to do this course until you pass?

Assuming we have a function called fail(), that returns TRUE when the experiment fails; then this code gives us number of trials until sucess

k <- 1
while(fail()) {
    k <- k+1
}

A simple way to model fail() is to use a coin

fail <- function() {
    coin <- sample(c("H","T"), size=1)
    return(coin != "H")
}

The symbol != means not equal

It will be useful to assume that heads and tails have different probabilities

You can imagine a deck of cards with cards of two colors in different proportion

Since the sum of probabilities for all outcomes is 1, and there are two possible outcomes, then

\[\Pr(\text{Tails}) = 1- \Pr(\text{Heads})\]

This function takes one input: the probability of success

fail <- function(p) {
    coin <- sample(c("H","T"), size=1, prob=c(p, 1-p))
    return(coin != "H")
}

roll_until_head <- function(p) {
    i<-1
    while(sample(c("H","T"), size=1,
            prob=c(p, 1-p)) != "H") {
        i <- i+1
    }
    return(i)
}

sim <- replicate(1000, roll_until_head(0.5))
barplot(table(sim)/1000)

sim <- replicate(1000, roll_until_head(0.9))
barplot(table(sim)/1000)

sim <- replicate(1000, roll_until_head(0.1))
barplot(table(sim)/1000)