Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- # see more: http://reakkt.cm
- n <- 1000
- mean.eq <- 0.001
- sd.eq <- 0.02
- contamination.ratio <- 0.01
- contaminated.sd <- 6
- cf1 <- 1
- cf2 <- 0.8
- ###
- x <- numeric(n)
- contaminatedN <- round(n*contamination.ratio)
- contaminated.idx <- sample(1:n,contaminatedN,replace=FALSE)
- x[contaminated.idx] <- rnorm(contaminatedN,mean.eq,contaminated.sd)/100
- sd.d <- numeric(n)
- for (i in 1:n) {
- if (is.na(match(i,contaminated.idx))) {
- sd.d[i] <- sum( (cf1/abs(i-contaminated.idx)^cf2) * abs(x[contaminated.idx]) )
- x[i] <- rnorm( 1,0,sd.eq + ifelse(sd.d[i]>sd.eq,sd.d[i],0) )
- } else sd.d[i] <- sd.d[i-1]
- }
- par(mfrow=c(3,1))
- plot(cumsum(x),type="l")
- plot(x,type="l")
- plot(sd.d,type="l")
- abline(h=sd.eq,col="Red",lty="dotted")
- par(mfrow=c(1,1))
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement