#globalVariables(c(‘FLCohort’,‘m’,‘FLPar’,‘dims’))
##############################################################
#' noise
#'
#' A noise generator
#'
#' @param n number of iterations
#' @param len an \code{FLQuant}
#' @param sd standard deviation
#' @param b autocorrelation parameter, a real number in [0,1]
#' @param trunc get rid of first values equal to trunc, i.e. to allow burn in
#'
#' @export
#' @docType methods
#' @rdname noise
#'
#' @return A vector with autocorrelation equal to B.
#'
#' @references Ranta and Kaitala 2001 Proc. R. Soc.
#' vt = b * vt-1 + s * sqrt(1 - b^2)
#' s is normally distributed random variable with mean = 0
#' b is the autocorrelation parameter
#' @export
#'
#' @examples
#' \dontrun{
#' white <- noise(1000,sd=.3,b=0)
#' plot(white)
#' acf(white)
#'
#' red <- noise(1000,sd=.3,b=0.8)
#' plot(red)
#' acf(white)
#'
#' blue <- noise(1000,sd=.3,b=-0.8)
#' plot(blue)
#' acf(blue)
#' }
setGeneric('noise', function(n,len,...) standardGeneric('noise'))
# n =100
# len=FLQuant(0,dimnames=list(year=1:55))
# b =0
# sd =0.3
# trunc=0
noiseFn<- function(len,sd=1,b=0,trunc=0){
mn=0
x <- rep(0, len+1) # going to hack off the first value at the end
s <- rnorm(len,mean=mn,sd=sd)
for(i in 1:len){
x[i+1] <- b * x[i] + s[i] * sqrt(1 - b^2)
if(trunc>0){
if (x[i+1] > (1-trunc)) x[i+1] <- ( 1-trunc)
if (x[i+1] < (-1+trunc)) x[i+1] <- (-1+trunc)}
}
x<-x[-1]
return(x)}
# setMethod("noise", signature(n='numeric', len="missing"),
# function(n,len))
setMethod("noise", signature(n='numeric', len="FLQuant"),
function(n=n,len=len,sd=0.3,b=0) {
if (dims(len)$iter!=1) stop("len can not have iter>1")
res=aaply(len,c(1,3:5), function(x)
maply(seq(n), function(x) noiseFn(dims(len)$year,b=b,sd=sd,trunc=0)))
ln=length(dim(res))
names(dimnames(res))[[ln]]="year"
names(dimnames(res))[[ln-1]]="iter"
if (ln==2)
res=aperm(res,c(ln,ln-1))
else
res=aperm(res,c(seq(ln-2),ln,ln-1))
dmns=dimnames(len)
dmns$iter=seq(n)
FLQuant(c(res),dimnames=dmns)})
## cohort effects
coEff=function(x,sd=1,B=0){
dev =log(rlnorm(length(dmns$cohort),0,cv))
for(i in 2:(length(dev)))
dev[i]=dev[i]+dev[i-1]*rho
parC[var]=parC[var]*exp(dev)
tmp=len2wt(parC,vonB(ages(FLCohort(m(stk))),parC[c("linf","t0","k")]))
res=window(as(tmp,"FLQuant"),start=1,end=dims(m(stk))$year)
res}
cEff=function(stk,par,sigma,rho=0,var="k"){
cv=sigma #getS(sigma,rho)
dmns =dimnames(par)
dmns$cohort=dimnames(FLCohort(m(stk)))$cohort
parC =FLPar(rep(c(par),length(dmns$cohort)),dimnames=dmns[c(1,3,2)])
units(parC)=""
dev =log(rlnorm(length(dmns$cohort),0,cv))
for(i in 2:(length(dev)))
dev[i]=dev[i]+dev[i-1]*rho
parC[var]=parC[var]*exp(dev)
tmp=len2wt(parC,vonB(ages(FLCohort(m(stk))),parC[c("linf","t0","k")]))
res=window(as(tmp,"FLQuant"),start=1,end=dims(m(stk))$year)
res}
#Spectral analysis function
spectra <- function(x,fs=1,norm = FALSE, pl = TRUE,omit=-(1:5))
{
# Pad x with zeroes to make it's length a power of 2, i.e. length should be 2^something
# This makes the fft faster
oldx <- x # keep for later
if(norm == TRUE) x <- x - mean(x)
nfft <- (2^ceiling(log2(length(x))))
x[(length(x)+1): nfft] <- 0
fftx <- fft(x)
# It's symmetrical so throw away second half. Only first 1 + nfft points are unique
NoUnPo <- ceiling((nfft+1)/2) # Number of unique points
fftx <- fftx[1:NoUnPo]
# First element is DC component, last is the Nyquist component
# Take magnitude of fft of x and scale the fft so that it is not a function of length
mx <- abs(fftx) / length(x)
# Take square of magnitude
mx <- mx^2
# As we dropped the first half of fft so multiply by 2 to keep same energy
# DC component (first element) and Nyquist component (last element if even
# number points (should be)) should not be multiplied by 2
mx[2:(length(mx)-1)] <- mx[2:(length(mx)-1)] * 2
f <- seq(from =0, to = NoUnPo-1) * fs/nfft # frequency axis for plot
return(as.data.frame(list(mx = mx, f = f)))}
#ggplot(spectra(x))+geom_line(aes(f,mx))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.