################################################################################
# Stat-Tutorial-01-deterministicsXrandomness                                   
#                                                                              
# Generation of deterministic apparent randomness.                             
#                                                                              
# R. Labouriau                                                                 
# Last revised: Spring 2020     
#
# Copyright © 2020 by Rodrigo Labouriau
################################################################################

# Background idea:
# We will explore examples of sequences of numbers genetated by repeatedly
# applying a function many times. That is, we will first define a function,
# for example f(x) = x + 1, then we choose a starting value, for example 3,
# next we apply the function to the starting value, e.g. f(3) = 3 + 1 = 4,
# we next apply the function to the result of the last calculation, f(4) = 5,
# and then apply the function again to the result last, f(5) = 6, and continue
# this process many times. In the example of the function f(x) = x + 1, using 
# the initial value 2 we get the sequence 2, 3, 4, 5, ... ; indeed a regular
# sequence.
# In many cases the sequence obtained is rather regular and predictable, but not
# always!
# 
# Please note that the sequence is deterministic, in the sence that knowing
# the initial value and the function we can reconstruct the sequence. Moreover,
# we can always predict the next element of the sequence. Furthermore, in many
# cases we can guess the function by inspecting a piece of the sequence.
# For example, if I present a sequence 2, 4, 16, 256, ... you may realize that 
# the function f(x) = x^2 with the initial value 2 would generate the sequence.
# In conclusion, the results of sequences generated in this way are NOT random!

# First we define a simple well behaved sequence generated by f(x)=x^2

# We define the generating function
f <- function(x){return(x^2)}

# We test the value of the function at several values
f(1/2)

f(1/4)

# We will next apply the function sucessively to its results, i.e. iterate
y <- f(1/2)
y
y <- f(2)
y

y <- f(y)
y

y <- f(y)
y

# We will iterate 10 times keeping the results in the vector Y
# and plot the results.
Y <- numeric(10)

Y[1] <- 0.5
for(i in 1:9){
  Y[i+1] <- f(Y[i])
}

Y

plot(Y)

# We improve the plot
plot(Y, type="b", lwd=2, cex=0.5, pch=19, xlab="Iteration")


# We iterate the function 10 times, but starting at a different values
Y[1] <- 0.9
for(i in 1:9){
  Y[i+1] <- f(Y[i])
}

Y

plot(Y, type="b", lwd=2, cex=0.5, pch=19, xlab="Iteration")

# Another starting value
Y[1] <- 0.99
for(i in 1:9){
  Y[i+1] <- f(Y[i])
}

plot(Y, type="b", lwd=2, cex=0.5, pch=19, ylim=c(0,1), xlab="Iteration")

# Now, we will superpose the representation of several sequence defined
# with different starting points.
Y[1] <- 0.9
for(i in 1:9){
  Y[i+1] <- f(Y[i])
}
points(Y, type="b", lwd=2, cex=0.5, pch=19, col=2)

Y[1] <- 0.5
for(i in 1:9){
  Y[i+1] <- f(Y[i])
}
points(Y, type="b", lwd=2, cex=0.5, pch=19, col=3)

#################################################


# We study now a sequence generated by a different function: 
# f(x) = 2x^2 - 1

# Next we define the function and evaluate it at several points
f <- function(x){return(2*x^2-1)}

f(1/2)

y <- f(1/2)
y

y <- f(y)
y

Y <- numeric(50)
Y[1] <- 1/2
for(i in 1:49){
  Y[i+1] <- f(Y[i])
}

Y

plot(Y , type='b', cex=0.5, pch=19)

# Starting at 1/2 the sequence behaves very regularly, indeed!
# What happens if we start at a value VERY close to 1/2, but NOT equal to 1/2 ?
# Let us see ... 

# Here are the first 50 values of the function
Y[1] <- 0.50001
for(i in 1:49){
  Y[i+1] <- f(Y[i])
}

Y

plot(Y , type='b', cex=0.5, pch=19)

plot(Y , type='l', lwd=2)

Y[50]

# Here we continue iterating getting the next 50 elements of the sequence
Y[1] <- Y[50]
for(i in 1:49){
  Y[i+1] <- f(Y[i])
}

plot(Y , type='l', lwd=2)

# Continue iterating by executing the last commands over and over and watch ...
# What happens? Do you get a regular pattern?

# What happens if you try other values different starting values than 1/2 ?
# (but between -1 and 1)
Y[1] <- 0.89
for(i in 1:49){
  Y[i+1] <- f(Y[i])
}
plot(Y , type='l', lwd=2)

Y[1] <- Y[50]
for(i in 1:49){
  Y[i+1] <- f(Y[i])
}

plot(Y , type='l', lwd=2)

# Try other initial values ...

# Suppose that someone just present the results for you
# whithout showing the function that generated the "regular" sequence.
# Would you be able to detect a regular pattern?

# THE END !

################################################################################

################################################################################
# OPTIONAL PRACTICE, here we explore a bit more ...
################################################################################

n <- 500
Y <- numeric(n)
Y[1] <- 0.62
for(i in 1:(n-1)){
  Y[i+1] <- f(Y[i])
}

plot(Y , type='l')

hist(Y)

#############
n <- 5000
Y <- numeric(n)
Y[1] <- 0.62
for(i in 1:(n-1)){
  Y[i+1] <- f(Y[i])
}

hist(Y)

#############
n <- 10000
Y <- numeric(n)
Y[1] <- 0.62
for(i in 1:(n-1)){
  Y[i+1] <- f(Y[i])
}

hist(Y, freq=F, col="lightblue")

#############
n <- 100000
Y <- numeric(n)
Y[1] <- 0.62
for(i in 1:(n-1)){
  Y[i+1] <- f(Y[i])
}

hist(Y, freq=F, col="lightblue")

#############
n <- 100000
Y <- numeric(n)
Y[1] <- 0.501
for(i in 1:(n-1)){
  Y[i+1] <- f(Y[i])
}

hist(Y, freq=F, col="lightblue")

#############
n <- 100000
Y <- numeric(n)
Y[1] <- 0.50001
for(i in 1:(n-1)){
  Y[i+1] <- f(Y[i])
}

hist(Y, freq=F, col="lightblue")

# THE second END!