power.determination: Power Determination

In connolr3/eegdatasim: Tool for EEG signals and ERPs

Power Determination

Description

Obtains an estimate of bernoulli power based on monte carlo estimation.

Obtains this power in terms of ERP frequency estimation, and ERP ampltide estimation.

Usage

power.determination(0.1,7,5,200,250,40,50)

Arguments

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

Value

results

dataframe containing ps for ERP frequency estimations in first column, and ps for ERP amplitude estimations in second.

Note

Has a significant run-time.

Author(s)

Rose Connolly [e-mail]

Examples

#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))
}

connolr3/eegdatasim documentation built on April 14, 2022, 10:39 a.m.