<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 State assumptions about your data
Let X<sub>D</sub> denote features of depressed people, X<sub>ND</sub> denote features of people that are not depressed. Let μ<sub>D</sub> be the mean of X<sub>D</sub> and μ<sub>ND</sub> be the mean of X<sub>ND</sub>.
Assume the real X<sub>D</sub> and X<sub>ND</sub> are both from a normal distribution.
Formally define statistical test
The null and alternative hypotheses are
Step1:  Note
Step2:  Sample data from null
Step3:  Sample data from alternate
Step4:  Plot power vs n on null set
Step5:  Compute p-value on your real data
Step6:  Sample from real data and plot power
<ASSISTANT_TASK:>
Python Code:
from scipy.stats import pearsonr
from scipy.stats import chisquare
from scipy.stats import ks_2samp
from scipy.stats import anderson_ksamp
pearsonr_test = lambda x: pearsonr(x[0], x[1])[1]
chi_test = lambda x: chisquare(x[0], x[1])[1]
ks_test = lambda x: ks_2samp(x[0], x[1])[1]
anderson_ksamp_test = lambda x: anderson_ksamp(x)[2]
Explanation: State assumptions about your data
Let X<sub>D</sub> denote features of depressed people, X<sub>ND</sub> denote features of people that are not depressed. Let μ<sub>D</sub> be the mean of X<sub>D</sub> and μ<sub>ND</sub> be the mean of X<sub>ND</sub>.
Assume the real X<sub>D</sub> and X<sub>ND</sub> are both from a normal distribution.
Formally define statistical test
The null and alternative hypotheses are:
H<sub>0</sub>: μ<sub>D</sub> = μ<sub>ND</sub>
H<sub>A</sub>: μ<sub>D</sub> != μ<sub>ND</sub>
Provide algorithm for implementing test
We use Kolmogorov-Smirnov Goodness-of-Fit Test, which is a two-sided test for the null hypothesis that 2 independent samples are drawn from the same continuous distribution.
The original 754 features was reduced to 27 features, and on each of the 27 features, K-S test is applied below.
End of explanation
import itertools
np.random.seed(123456789) # for reproducibility, set random seed
alpha = 0.05 
r = 20  # number of rois
N = 100 # number of samples at each iteration
# define number of subjects per class
S = np.array((4, 6, 8, 10, 14, 18, 20, 26, 30, 40,
              50, 60, 70, 80, 100, 120, 150, 200, 250,
              300, 400, 500, 750, 1000, 1500, 2000,
              3000, 5000))
Explanation: Note: Generate random samples and plot power vs n on null set based on Greg's code.
Random Sampling Setup
End of explanation
pow_null = np.array((), dtype=np.dtype('float64'))
# compute this statistic for various sizes of datasets
for s in S:
    s0 = s/2
    s1 = s - s0
    # compute this many times for each operating point to get average
    pval = np.array((), dtype=np.dtype('float64'))    
    for _ in itertools.repeat(None,N):
        g0 = 1 * (np.random.rand(r, r, s0) > 0.5) # (null), 0.52 (classes)
        g1 = 1 * (np.random.rand(r, r, s1) > 0.5) # (null), 0.48 (classes)
        # compute feature of data
        pbar0 = 1.0*np.sum(g0, axis=(0,1))/(r**2 * s0)
        pbar1 = 1.0*np.sum(g1, axis=(0,1))/(r**2 * s1)
        # compute K-S test on feature
        pval = np.append(pval, ks_2samp(pbar0, pbar1)[1])
    
    # record average p value at operating point
    pow_null = np.append(pow_null, np.sum(1.0*(pval < alpha))/N)
Explanation: Sample data from null
End of explanation
pow_alt = np.array((), dtype=np.dtype('float64'))
# compute this statistic for various sizes of datasets
for s in S:
    s0 = s/2
    s1 = s - s0
    # compute this many times for each operating point to get average
    pval = np.array((), dtype=np.dtype('float64'))    
    for _ in itertools.repeat(None,N):
        g0 = 1 * (np.random.rand(r, r, s0) > 0.52) # (null), 0.52 (classes)
        g1 = 1 * (np.random.rand(r, r, s1) > 0.48) # (null), 0.48 (classes)
        # compute feature of data
        pbar0 = 1.0*np.sum(g0, axis=(0,1))/(r**2 * s0)
        pbar1 = 1.0*np.sum(g1, axis=(0,1))/(r**2 * s0)
        # compute K-S test on feature
        pval = np.append(pval, ks_2samp(pbar0, pbar1)[1])
    
    # record average p value at operating point
    pow_alt = np.append(pow_alt, np.sum(1.0*(pval < alpha))/N)
Explanation: Sample data from alternate
End of explanation
plt.scatter(S, pow_null, hold=True, label='null')
plt.scatter(S, pow_alt, color='green', hold=True, label='alt')
plt.xscale('log')
plt.xlabel('number of samples')
plt.ylabel('power')
plt.title('Strength of depression classification under null model')
plt.axhline(alpha, color='red', linestyle='--', label='alpha')
plt.legend(loc=5)
plt.show()
Explanation: Plot power vs n on null set
End of explanation
p_vals = list()
for a, b in zip(depr, not_depr):
    p_vals.append(round(ks_2samp(a, b)[1], 5))
print p_vals
Explanation: Compute p-value on your real data
End of explanation
from skbio.stats.power import subsample_power
# Computer power of a sub sample set
def compute_sub_power(test, samples):
    pwr, counts = subsample_power(test=test,
                                   samples=samples,
                                   max_counts=1205,
                                   min_counts=100,
                                   counts_interval=100,
                                   draw_mode="ind",
                                   alpha_pwr=0.05)
    return pwr, counts
from mpl_toolkits.axes_grid1 import Grid
plt.close('all')
fig = plt.figure()
fig.set_size_inches(18, 9)
grid = Grid(fig, rect=111, nrows_ncols=(4, 7),
            axes_pad=0.25, label_mode='L',
            )
    
def plot_power(i, ax):
    a, b = depr[i], not_depr[i]
    samples = [np.array(a), np.array(b)]
    pwr, counts = compute_sub_power(ks_test, samples) 
    ax.plot(counts, pwr.mean(0))
    
for i, ax in enumerate(grid):
    if i < 27:
        plot_power(i, ax)
        title = 'p = ' + str(p_vals[i])
        ax.set_title(title)
plt.tight_layout()
Explanation: Sample from real data and plot power
End of explanation
<END_TASK>