function ds_stat=cosmo_montecarlo_phase_stat(ds,varargin)
% compute phase statistics based on Monte Carlo simulation
%
% ds_stat=cosmo_montecarlo_phase_stat(ds,...)
%
% Inputs:
% ds dataset struct with fields:
% .samples PxQ complex matrix for P samples (trials,
% observations) and Q features (e.g. combinations
% of time points, frequencies and channels)
% .sa.targets Px1 array with trial conditions.
% There must be exactly two conditions, thus
% .sa.targets must have exactly two unique
% values. A balanced number of samples is
% requires, i.e. each of the two unique values in
% .sa.targets must occur equally often.
% .sa.chunks Px1 array indicating which samples can be
% considered to be independent. It is required
% that all samples are independent, therefore
% all values in .sa.chunks must be different from
% each other
% .fa } optional feature attributes
% .a } optional sample attributes
% 'output',p Return statistic, one of the following:
% - 'pbi': phase bifurcation index
% - 'pos': phase opposition sum
% - 'pop': phase opposition product
% 'niter',niter Generate niter null datasets by random
% shuffling the targets.
% If you have no idea what value to use, consider
% using niter=1000.
% 'zscore',z (optional, default='non_parametric')
% Compute z-score using either:
% 'non_parametric': non-parametric approach based
% on how many values in the
% original dataset show an
% output greater than the
% null data.
% 'parametric' : parametric approach based on
% mean and standard deviation
% of the null data. This would
% assume normality of the
% computed output statistic.
% 'extreme_tail_set_nan',n (optional, default=true)
% If n==true and all the output value in the
% original dataset for a particular feature i is
% less than or greater than all output values for
% that feature in the null datasets, set the
% output in stat.samples(i) to NaN.
% If n==false, the p value corresponding to the
% output is limited to the range
% [1/niter,1-1/niter]
% 'progress',p (optional,default=1)
% Show progress every p null datasets.
% 'permuter_func',f (optional, default is function defined in
% this function's body)
% Function handle with signature
% idxs=f(iter)
% which returns permuted indices in the range
% 1:nsamples for the iter-th iteration with that
% seed value. The targets are resamples using
% these permuted indices.
% 'seed',s (optional, default=1)
% Use seed s when generating pseudo-random
% permutations for null distribution.
%
% Output:
% stat_ds Dataset with field
% .samples 1xQ z-scores indicating the probability of the
% observed data in ds.samples, under the null
% hypothesis of no phase difference. z-scores are
% not corrected for multiple comparisons.
%
% Notes:
% - this function computes phase statistics for each feature separately;
% it does not correct for multiple comparisons
% - p-values are computed by dividing as (r+1) / (niter+1), with r the
% number of times that the original data as less then the null
% distributions. This follows the recommendation of North et al (2002).
%
% Reference
% - North, Bernard V., David Curtis, and Pak C. Sham. "A note on the
% calculation of empirical P values from Monte Carlo procedures." The
% American Journal of Human Genetics 71.2 (2002): 439-441.
%
% See also: cosmo_phase_stat, cosmo_phase_itc
