#Benjamin Preis
#MAS.864, Spring 2020
#Problem Set 3, Random Systems

library(ggplot2)

###6.2c
rannum<-function(x){
  x<-(8121*x+28411)%%134456
  return(x)
}

x1<-c()
x2<-c()
x01<-1
x02<-2
for(i in 1:1000){
  if(i==1){
    x1[i]<-rannum(x01)
    x2[i]<-rannum(x02)
  }
  else{
    x1[i]<-rannum(x1[i-1])
    x2[i]<-rannum(x2[i-1])
  }
}
x1<-x1/134456
x2<-x2/134456

y1<-sqrt(-2*log(x1))*sin(x2)
y2<-sqrt(-2*log(x1))*cos(x2)

#C1
mean(c(y1,y2))
#C2
sd(c(y1,y2))^2


###6.4
rannum<-function(x){
  x<-(8121*x+28411)%%134456
  return(x)
  
seeds<-rannum(c(1:10))

position<-list()

temp<-c()
for(j in 1:10){
  for(i in 1:1000){
    if(i ==1){
      x<-rannum(seeds[j])
      temp[i]<-0
    }
    else{
      x<-rannum(x)
      
      if(x/134456>.5){
        temp[i]<-temp[i-1]+1
      }
      else{
        temp[i]<-temp[i-1]-1
      }
    }
  }
  position[[j]]<-temp
}

position<-lapply(position, function(x) x<-data.frame("x"=c(1:1000),"y"=x))

g<-ggplot()

for(i in 1:10){
  g<-g+geom_line(data=position[[i]],aes(x=x,y=y))
}
g<-g+xlab("t")+ylab("Position")

df2<-data.frame("x"=100*c(1:10))
df2$ymin<-(-1)*3*sqrt(df2$x)
df2$ymax<-3*sqrt(df2$x)
g<-g+geom_errorbar(data=df2,aes(x=x,ymin=ymin,ymax=ymax))