Robust standard deviation and median of outliers

rstd.py

import pandas as pd
import collections
import numpy as np
import matplotlib.pyplot as plt

from allensdk.core.brain_observatory_nwb_data_set import BrainObservatoryNwbDataSet

df = pd.read_csv('/allen/aibs/technology/nicholasc/brain_observatory_cell_df.csv', sep=' ').drop('Unnamed: 0', axis=1)
# print df.columns

def robust_std(x):
    '''
    Robust estimate of std of noise in df/f
    '''

    # first pass removing big pos peaks
    x = x[x< 1.5*np.abs(x.min())]
    MAD = np.median(np.abs(x - np.median(x)))
    rstd = 1.4826*MAD

    # second pass removing remaining pos and neg peaks
    x = x[abs(x) < 2.5*rstd]
    MAD = np.median(np.abs(x - np.median(x)))
    return 1.4826*MAD

def get_median_deviation(dff, rstd):
    
    return np.median(dff[dff>rstd])

data_dict = collections.defaultdict(list)
for file_name in df['file_name'].unique()[:10]:
    data_set = BrainObservatoryNwbDataSet(file_name)
    session_type = data_set.get_session_type()
    timestamps, dff_array = data_set.get_dff_traces()
    csids = data_set.get_cell_specimen_ids()
    
    for curr_dff, curr_csid in zip(dff_array, csids):
        curr_rstd = robust_std(curr_dff)
        median_deviation = get_median_deviation(curr_dff, curr_rstd)

        data_dict['csid'].append(curr_csid)
        data_dict['rstd'].append(curr_rstd)
        data_dict['median_deviation'].append(median_deviation)
        data_dict['session_type'].append(session_type)


rstd_analysis_df = pd.merge(pd.DataFrame(data_dict), df, on=['csid', 'session_type'])
rstd_analysis_df.to_csv('/home/nicholasc/projects/brain_observatory_analysis/ephys_receptive_field/rstd_analysis_df.csv', sep=' ', index=False)