| power.determination | R Documentation |
Obtains an estimate of bernoulli power based on monte carlo estimation.
Obtains this power in terms of ERP frequency estimation, and ERP ampltide estimation.
power.determination(0.1,7,5,200,250,40,50)
accuracy_window |
percentage used to caluclate range acceptable estimate is found in |
freq |
expected peak frequency |
amp |
expected peak amplitude |
frames |
Frames in each trial |
srate |
Maximum number of trials you wish to calculate p for, default = 40 |
srate |
number of trials for wish to avaergae over to estimate p for each trial, default = 100 |
results |
dataframe containing ps for ERP frequency estimations in first column, and ps for ERP amplitude estimations in second. |
Has a significant run-time.
Rose Connolly connolr3@tcd.ie
#NOT RUN
# my.ps<-power.determination(0.1,7,5,200,250,30,50)
#plot(my.ps$my_frequency_ps,main="p estimates for 1 to 30 trials averaged over 50 simulations")
## The function is currently defined as
power.determination <-function(accuracy_window,freq,amp,frames,srate,maxtrial = 40,averagingN=100)
{
set.seed("123")
my_frequency_ps <- numeric(length=maxtrial)
my_amplitude_ps <- numeric(length=maxtrial)
meanpower <- R.matlab::readMat("meanpower.mat")$meanpower
for (N in 1:maxtrial)
{
freq_ones <- 0
amp_ones <- 0
for (j in 1:averagingN)
{
#create sample data
mysignal <-noise(frames, N, 250, meanpower=meanpower) + amp * peak(frames, N, srate, freq)
#prepare signal: average and standardise
my_averaged_signal <- signal_averaging(mysignal, frames, N)
hats <- est_sig_hat(my_averaged_signal, frames / 2)
standata <- (my_averaged_signal - hats[2]) / hats[1]
#estimate peak range
mypeak_range <- find_ERP_range(standata, 1.7)
#prepare data and starting values for optim
yis <- my_averaged_signal[mypeak_range]
peak_center_estimate = which(my_averaged_signal == max(my_averaged_signal))
pars <-optimise_ERP(mypeak_range, yis ,mysr = srate,pkcntr = peak_center_estimate)
if (abs(freq - pars$par[1]) <= freq * accuracy_window)
freq_ones=freq_ones+ 1
if (abs(amp - pars$par[2]) <= amp * accuracy_window)
amp_ones=amp_ones+ 1
}
cat("\n Completed (no. trials): ",N)
freq_p <- freq_ones / averagingN
amp_p <- amp_ones / averagingN
my_frequency_ps[N]<-freq_p
my_amplitude_ps[N]<-amp_p
}
results<-data.frame(my_frequency_ps,my_amplitude_ps)
return((results))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.