######################################################
## PS-239 - Third R class - Fall 2006
######################################################

# 1. If statements

if(condition) {
  statement or block of statements
}

# If the condition is true, the statements are executed
# If the condition is false, the statements are not executed

matrix1<-matrix(c(3,0,0,0,3,0,0,0,3),nrow=3,ncol=3)
dim(matrix1)

# Function to check whether the number of cols equals the number of rows
agrees<-function(x) {
  error=(ncol(x)!=nrow(x))
  if (error==TRUE)  return("matrix is not square")
}
agrees(matrix1) # Does nothing because the condition is false

matrix2<-matrix(c(3,0,0,0,1,5),nrow=2,ncol=3)
dim(matrix2)
matrix2
agrees(matrix2)

# We can also add an "else" condition ==> whenever the condition in the if is false,
# it executes the statements following the else

agrees2<-function(x) {
  error=(ncol(x)!=nrow(x))
  if (error==TRUE)  return("matrix is not square")
  else return("matrix is square")
}
agrees(matrix1) # Does nothing because there is no else condition
agrees2(matrix1) # Does something


# Now let's use the if statement for something useful
# We'll invert a matrix but only if it's square

invert<-function(x) {
  error=(ncol(x)!=nrow(x))
  if (error==TRUE)  return("matrix is not square")
  else inver1<-solve(x)
  list(inver1=inver1)
}
invert(matrix1)
invert(matrix2)

# If you want to do more than one statement ==> use ""{}"

invert2<-function(x) {
  error=(ncol(x)!=nrow(x))
  if (error==TRUE)  {
    return("matrix is not square")
  }
  else {
    inver1<-solve(x)
    inver2<-ginv(x)
  }
  list(inver1=inver1,inver2=inver2)
}

invert2(matrix1)
invert2(matrix2)

# 2. Loops

# 2.1 For loops

for (i in start.number:end.number) {
  statement or block statements
}

x<-rnorm(25,mean=0,sd=1)
cum.sum=0
for (i in 1:25) {
  cum.sum=x[i]+cum.sum
}
cum.sum
sum(x)  # Gives the same number!!
mean.x<-cum.sum/length(x)
mean.x
mean(x)

# 2.2 While loops
x<-1
while(x<4) {   # Tests at the beginning
  x<-x+0.5
}
x

# 2.2 Do -While loops

do {

while }

######## MONTECARLO ######################

# Example of Monte Carlo simulations
# Example taken from Tom Huber, Physics Department, Gustavus Adolphus College. (http://physics.gac.edu/~huber/envision/instruct/montecar.htm)
# Square that has one corner at the origin of a coordinate system and has sides of length 1 ==>area=1
# Inscribe a quarter of a circle with radius inside this square ==> area=pi/4.
# We can do a Monte Carlo simulation to find the relative area of the circle and square and then multiply the circle's area by 4 to find pi.
# For a point (X,Y) to be inside of a circle of radius 1, its distance from the origin (X2+Y2) will be less than or equal to 1.
# We can generate thousands of random (X,Y) positions and determine whether each of them are inside of the circle.
# Each time it is inside of the circle, we will add one to a counter.
# After generating a large number of points, the ratio of the number of points inside the circle to the total number of points
# generated will approach the ratio of the area of the circle to the area of the square.
# Therefore pi is approximately equal to 4*(Number of Points inside circle)/ (Total Points Generated)

nsimul<-100
cum<-0
for (i in 1:nsimul) {
  x<-runif(1)
  y<-runif(1)
  distan<-sqrt(x^2+y^2)     # Distance of point (x,y) from the origin
  if (distan<=1) cum<-cum+1
}
4*cum/nsimul

# Using matrices
nsimul<-100000
cum<-rep(FALSE,nsimul)
x<-runif(nsimul)
y<-runif(nsimul)
distan<-sqrt(x^2+y^2)
cum[(distan<=1)]<-TRUE
4*sum(as.numeric(cum))/nsimul