#M480 R functions



bmot<-function(n=100){
#approximates a Brownian motion sample path on [0,100] as n gets large
h <- 10/sqrt(n)
tem<-rnorm(n,mean=0,sd=h)
B<-c(0,cumsum(tem))
times <- 100*0:n/n  
plot(times,B,type="l") 
if(n < 1000) 
list(B=B) 
}



#From Olive (2014, 8.1)

cltsim <- function(n=100, nruns=100){
ybar <- 1:nruns
for(i in 1:nruns){
  ybar[i] <- mean(rexp(n))}
list(ybar=ybar)}


cltsim2 <- function(n=100, nruns=100){
ybar <- 1:nruns
for(i in 1:nruns){
  ybar[i] <- mean(rnorm(n))}
list(ybar=ybar)}


poisprocess <- function(n=100, rate=1){
#generates a Poisson process with n events
#out<-poisprocess(n=100,rate=1/20) #rate = lambda
wtimes <- c(0,cumsum(rexp(n,rate=rate)))
Nt <- 0:n
plot(wtimes,Nt,type="s")  #should have slope near rate
abline(0,rate)
list(wtimes=wtimes,Nt=Nt) 
}


rng<-function(n=100,seed=12345){
#uniform(0,1) random number generator
u<-1:(n+1)
u[1] <- seed
for(i in 2:(n+1)){
u[i] <- (69069*u[i-1] + 1) %% (2^32) 
}
u<-u/(2^32)
u <- u[-1]
list(u=u)
}


rwalk <- function(n=100, type = 1, y0 = 1){
#generates a random walk
#type 1) N(1,1), 2) Cauchy(1,1), 3) EXP(1), 4) uniform(0,2)
times <- 1:n
if(type==1)
tem<-rnorm(n,mean=1,sd=1)
if(type==2) 
tem<-rcauchy(n,location=1,scale=1) #generator might be bad
if(type==3) 
tem<-rexp(n)
if(type==4) 
tem<-runif(n,min=0,max=2)
tem<-c(y0,tem)
yrw<-cumsum(tem)[-1]
plot(times,yrw,type="l")  
abline(y0,1)
list(yrw=yrw) 
}


shorth3<-function(y, alpha = 0.05){
# computes the shorth(c) interval [Ln,Un] for c = cc.
#shorth lists Ln and Un using Frey's correction
n <- length(y)
cc <- ceiling(n * (1 - alpha + 1.12*sqrt(alpha/n)))
cc <- min(n,cc)
  sy <- sort(y)
  rup <- sy[cc]
  rlow <- sy[1]
  olen <- rup - rlow
  if(cc < n){
    for(j in (cc + 1):n){
	zlen <- sy[j] - sy[j - cc + 1]
	if(zlen < olen) {
	   olen <- zlen
	   rup <- sy[j]
	   rlow <- sy[j - cc + 1]
	}
    }
} 
Ln <- rlow; Un <- rup
list(shorth=c(Ln, Un))
}