<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Constant in space
Possibilities include 
* Bayesian inference using 2 zones which leads to no siginificant difference.
* fitting a line, should give gradient not different from 0
    * problem here is that sigma is unknown
* Test for normality (should give high GOF)
Step1:  Test two zones
Step2:  Test for gradient
Step3:  Calculate the relative change in density per milimeter
Step4:  Constant
Step5:  Compare AICc
The AICc of 2 zone model is probably wrong (switchpoint changes Likelihood only every 50 $\mu m$
Step6:  Normtest
Step7:  Minimal p-Value = 0.041 -> no strong indication for non-Gaussuian distribution of density in space
Constant in time
Step8:  p-Value of 0.08 does not indicate significant difference between the samples from different time points.
Step9:  Calculate mean density
<ASSISTANT_TASK:>
Python Code:
burn = 100000
iter_ = 1100000
thin = 1000
Explanation: Constant in space
Possibilities include 
* Bayesian inference using 2 zones which leads to no siginificant difference.
* fitting a line, should give gradient not different from 0
    * problem here is that sigma is unknown
* Test for normality (should give high GOF)
End of explanation
def make_2zone_model(data, SOX2_min, SOX2_max, sigma_min, sigma_max):
    
    positions = sp.array(data['pos'])
    values_SOX2 = sp.array(data['SOX2'])
   
    def step_function(x, switchpoint, left_value, right_value):
        ''' This function should return something in the same format as the passed array 
          Specifically, it produces an output that has an array of the same size of the experimental data
        but whose contents are the lower average until the switchpoint, and the upper average past the switchpoint.
        For all purposes, this builds the model to which we want to compare the data.
        '''
        return sp.where(x<=switchpoint, left_value, right_value)
    
    SOX2_left_mean = pymc.Uniform('SOX2_left_mean', lower = SOX2_min, upper = SOX2_max, value = 30.0)
    SOX2_right_mean = pymc.Uniform('SOX2_right_mean', lower = SOX2_min, upper = SOX2_max, value = 30.0)
    switchpoint = pymc.Uniform('switchpoint', lower = positions.min(), upper = positions.max(),\
                               value = 0.5*(positions.max() - positions.min()) + positions.min())
    sigma = pymc.Uniform('sigma', lower = sigma_min, upper = sigma_max, value = 2.0)
    
    @pymc.deterministic(name='SOX2_step')
    def SOX2_step(SOX2_left_mean = SOX2_left_mean, SOX2_right_mean = SOX2_right_mean):
        return SOX2_left_mean - SOX2_right_mean
    
    @pymc.deterministic(name='SOX2_mean')
    def SOX2_mean(x = positions, switchpoint = switchpoint, left_value = SOX2_left_mean, right_value = SOX2_right_mean):
        return step_function(x, switchpoint, left_value, right_value)
    
    cells_SOX2 = pymc.Normal('cells_SOX2', mu=SOX2_mean, tau = 1/sigma**2, value = values_SOX2,\
                                           plot = False, observed = True)
    return locals()
SOX2_min, SOX2_max = 0.0, 50.0
sigma_min, sigma_max = 0.1, 10.0
%%cache step_df.pkl step_df --cachedir ../cache/density
%%time 
step_df = pd.DataFrame()
for ID, data in SOX2_numbers.set_index('ID').dropna().groupby(level = 'ID'):
    map_ = pymc.MAP(make_2zone_model(data, SOX2_min, SOX2_max, sigma_min, sigma_max))
    map_.fit()
    M = pymc.MCMC(map_)
    M.sample(iter=iter_, burn=burn, thin=thin, progress_bar=False)
#     pymc.Matplot.plot(M.SOX2_step)
    trace = M.SOX2_step.trace()
    step_df.loc[ID, 'mean'] = trace.mean()
    step_df.loc[ID, 'minus3sigma'] = sp.percentile(trace, 0.5*(100.0-99.73))
    step_df.loc[ID, 'plus3sigma'] = sp.percentile(trace, 100.0 - 0.5*(100.0-99.73))
    step_df.loc[ID, 'minus2sigma'] = sp.percentile(trace, 0.5*(100.0-95.4))
    step_df.loc[ID, 'plus2sigma'] = sp.percentile(trace, 100.0 - 0.5*(100.0-95.4))
    step_df.loc[ID, 'minus_p05'] = sp.percentile(trace, 0.5*(100.0-97.5))
    step_df.loc[ID, 'plus_p05'] = sp.percentile(trace, 100.0 - 0.5*(100.0-97.5))
    
    step_df.loc[ID, 'AICc'] = map_.AICc
step_df.query('minus_p05 > 0 or plus_p05 < 0')
Explanation: Test two zones
End of explanation
def make_line_model(data, m_min, m_max, n_min, n_max, sigma_min, sigma_max):
    
    positions = sp.array(data['pos'])
    values_SOX2 = sp.array(data['SOX2'])
   
    def line(x, m, n):
        return m * x + n
    
    m = pymc.Uniform('m', lower = m_min, upper = m_max, value = 0.0)
    n = pymc.Uniform('n', lower = n_min, upper = n_max, value = 30.0)
    sigma = pymc.Uniform('sigma', lower = sigma_min, upper = sigma_max, value = 2.0)
    
    @pymc.deterministic(name='SOX2_mean')
    def SOX2_mean(x = positions, m = m, n = n):
        return line(x, m, n)
    
    cells_SOX2 = pymc.Normal('cells_SOX2', mu=SOX2_mean, tau = 1/sigma**2, value = values_SOX2,\
                                           plot = False, observed = True)
    return locals()
m_min, m_max = -1.0, 1.0
n_min, n_max = 0.0, 50.0
sigma_min, sigma_max = 0.1, 10.0
%%cache m_df.pkl m_df --cachedir ../cache/density
%%time 
m_df = pd.DataFrame()
for ID, data in SOX2_numbers.set_index('ID').dropna().groupby(level = 'ID'):
    map_ = pymc.MAP(make_line_model(data, m_min, m_max, n_min, n_max, sigma_min, sigma_max))
    map_.fit()
    M = pymc.MCMC(map_)
    M.sample(iter=iter_, burn=burn, thin=thin, progress_bar=False)
#     pymc.Matplot.plot(M.m)
    trace =  M.m.trace()
    m_df.loc[ID, 'mean'] = trace.mean()
    m_df.loc[ID, 'minus3sigma'] = sp.percentile(trace, 0.5*(100.0-99.73))
    m_df.loc[ID, 'plus3sigma'] = sp.percentile(trace, 100.0 - 0.5*(100.0-99.73))
    m_df.loc[ID, 'minus2sigma'] = sp.percentile(trace, 0.5*(100.0-95.4))
    m_df.loc[ID, 'plus2sigma'] = sp.percentile(trace, 100.0 - 0.5*(100.0-95.4))
    m_df.loc[ID, 'minus_p05'] = sp.percentile(trace, 0.5*(100.0-97.5))
    m_df.loc[ID, 'plus_05'] = sp.percentile(trace, 100.0 - 0.5*(100.0-97.5))
    m_df.loc[ID, 'AICc'] = map_.AICc
m_df.query('minus_p05 > 0 or plus_05 < 0')
Explanation: Test for gradient
End of explanation
rel_dens_change = 1000 * sp.absolute(m_df.query('minus_p05 > 0 or plus_05 < 0')['mean']) / SOX2_numbers.set_index('ID')['SOX2'].groupby(level = 'ID').agg('mean')
rel_dens_change
Explanation: Calculate the relative change in density per milimeter
End of explanation
def make_constant_model(data, SOX2_min, SOX2_max, sigma_min, sigma_max):
    
    positions = sp.array(data['pos'])
    values_SOX2 = sp.array(data['SOX2'])
    
    SOX2_mean = pymc.Uniform('SOX2_mean', lower = SOX2_min, upper = SOX2_max, value = 30.0)
    sigma = pymc.Uniform('sigma', lower = sigma_min, upper = sigma_max, value = 2.0)
    
    cells_SOX2 = pymc.Normal('cells_SOX2', mu=SOX2_mean, tau = 1/sigma**2, value = values_SOX2,\
                                           plot = False, observed = True)
    return locals()
SOX2_min, SOX2_max = 0.0, 50.0
sigma_min, sigma_max = 0.1, 10.0
%%cache c_df.pkl c_df --cachedir ../cache/density
%%time 
c_df = pd.DataFrame()
for ID, data in SOX2_numbers.set_index('ID').dropna().groupby(level = 'ID'):
    map_ = pymc.MAP(make_constant_model(data, SOX2_min, SOX2_max, sigma_min, sigma_max))
    map_.fit()
    c_df.loc[ID, 'AICc'] = map_.AICc
Explanation: Constant
End of explanation
AICc = pd.DataFrame()
AICc['constant'] = c_df['AICc']
AICc['line'] = m_df['AICc']
AICc['2zone'] = step_df['AICc']
AICc['constant<line'] = AICc['constant'] < AICc['line']
AICc
Explanation: Compare AICc
The AICc of 2 zone model is probably wrong (switchpoint changes Likelihood only every 50 $\mu m$
End of explanation
pvalues = SOX2_numbers[['ID', 'SOX2']].groupby('ID').agg(lambda x: sp.stats.normaltest(x).pvalue)
pvalues
Explanation: Normtest
End of explanation
SOX2_space_mean = SOX2_numbers.groupby('ID').mean()
grouped = SOX2_space_mean.groupby('time')
samples = grouped.apply(lambda group: sp.array(group['SOX2'])).tolist()
sp.stats.f_oneway(*samples)
Explanation: Minimal p-Value = 0.041 -> no strong indication for non-Gaussuian distribution of density in space
Constant in time
End of explanation
SOX2_space_time_mean = grouped.agg(['mean', 'sem']).reset_index()
fig, ax = plt.subplots(figsize=(25.0 / 25.4, 25.0 / 25.4))
# plt.errorbar(sp.array(SOX2_space_time_mean['time']),\
#              sp.array(SOX2_space_time_mean['SOX2', 'mean']),\
#              sp.array(SOX2_space_time_mean['SOX2', 'sem']),\
#              fmt = 'o', color = 'black')
ax.scatter(SOX2_space_mean['time'], SOX2_space_mean['SOX2'], color = 'black', s = 10)
ax.plot([-0.5, 8.5], 2 * [SOX2_space_time_mean['SOX2', 'mean'].mean()], color = 'black')
ax.fill_between([-0.5, 8.5], 2 * [SOX2_space_time_mean['SOX2', 'mean'].mean() + SOX2_space_time_mean['SOX2', 'mean'].std()], 2 * [SOX2_space_time_mean['SOX2', 'mean'].mean() - SOX2_space_time_mean['SOX2', 'mean'].std()], alpha = 0.2, color = 'black', label = 'conf.belt')
ax.set_xlim(-0.5, 8.5)
ax.set_xlabel('Time (days)')
ax.set_ylabel('# SOX2+ cells')
ax.set_yticks(sp.arange(26, 35, 2))
ax.set_yticks(sp.arange(26, 35, 2)+1, minor = True)
ax.set_ylim(26, 34)
plt.savefig('../../figure_plots/Fig2_stem_cells_time.svg')
plt.show()
Explanation: p-Value of 0.08 does not indicate significant difference between the samples from different time points.
End of explanation
Nsection = SOX2_space_time_mean['SOX2', 'mean'].mean()
Nsection
Nsection_sem = SOX2_space_time_mean['SOX2', 'mean'].sem()
Nsection_sem
# see Rodrigo Albos et al., 2015
lcell = 13.2
lcell_sem = 0.1
rho = Nsection / lcell
rho
rho_ci = rho * (lcell_sem / lcell + Nsection_sem / Nsection)
rho_ci
Explanation: Calculate mean density
End of explanation
<END_TASK>