How_many.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Determine how many FRBs are needed to distinguish mutliple populations of FRBs based on the trends observed with few objects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "import numpy as np\n",
    "import sklearn.cluster as cl\n",
    "from sklearn import metrics\n",
    "import scipy.stats as stats\n",
    "from sklearn import mixture\n",
    "from sklearn.decomposition import PCA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "params = np.loadtxt('quadratic_fit_parameters_20150908.txt')\n",
    "name, W, DM = np.loadtxt(\"FRBs.txt\", comments=\"#\", unpack=True, usecols = [0,1,2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def quad(x, a, b, c):\n",
    "    return (a*x*x + b * x + c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def norm_data(N,SC):\n",
    "    pop_A = quad(np.linspace(2,10,N/2), params[0],params[1],params[2])\n",
    "    pop_B = quad(np.linspace(1,5,N/2), params[3],params[4],params[5])\n",
    "    scatter_A = pop_A + SC*np.random.randn(len(pop_A))\n",
    "    scatter_B = pop_B + SC*np.random.randn(len(pop_B))\n",
    "    \n",
    "    total_DM = np.append(scatter_A, scatter_B)\n",
    "    #total_DM = np.append(total_DM, DM)\n",
    "    \n",
    "    total_W = np.append(np.linspace(2,10,N/2), np.linspace(1,5,N/2))\n",
    "    #total_W = np.append(total_W, W)\n",
    "    \n",
    "    total_W = np.log10(total_W)\n",
    "    total_DM = np.log10(total_DM)\n",
    "    \n",
    "    total_DM = np.nan_to_num(total_DM)\n",
    "    \n",
    "    #print np.linspace(2,10,N), scatter_A\n",
    "    #print np.linspace(2,10,N), scatter_B\n",
    "    \n",
    "    y_mean = np.mean(total_DM)\n",
    "    x_mean = np.mean(total_W)\n",
    "    y_diff = total_DM - y_mean\n",
    "    x_diff = total_W - x_mean\n",
    "    y_norm = y_diff/np.std(total_DM,ddof=1)\n",
    "    x_norm = x_diff/np.std(total_W,ddof=1)\n",
    "    X_norm = np.array(zip(x_norm,y_norm))\n",
    "\n",
    "    return X_norm\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def K_mean(X):\n",
    "    kmeans_model = cl.KMeans(n_clusters=2).fit(X)\n",
    "    labels = kmeans_model.labels_\n",
    "    sil = metrics.silhouette_score(X,labels)\n",
    "    return labels, sil"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def Hierarchical(X):\n",
    "    Ward_model =cl.AgglomerativeClustering(n_clusters = 2).fit(X)\n",
    "    labels = Ward_model.labels_\n",
    "    sil = metrics.silhouette_score(X,labels)\n",
    "    return labels, sil"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def GMM(X):\n",
    "    GMM_model = mixture.GMM(n_components=2, covariance_type='full')\n",
    "    GMM_model.fit(X)\n",
    "    labels = GMM_model.predict(X)\n",
    "    sil = metrics.silhouette_score(X,labels)\n",
    "    return labels, sil"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def checkEqual2(iterator):\n",
    "       return len(set(iterator)) <= 1\n",
    "# eps=1.3,min_samples=3\n",
    "def DBScan(X):\n",
    "    #DBSCAN_model = cl.DBSCAN(eps=1.3).fit(X)\n",
    "    tmp = cl.DBSCAN(eps=1.3*np.sqrt(10./len(X)))#, min_samples=3.*len(X)/10.)\n",
    "    DBSCAN_model = tmp.fit(X)\n",
    "    #print tmp.get_params(\"eps\")\n",
    "    labels = DBSCAN_model.labels_ \n",
    "    n_clusters = len(set(labels))-(1 if -1 in labels else 0) \n",
    "    if n_clusters > 1:\n",
    "    #print checkEqual2(labels)\n",
    "    #if checkEqual2(labels) == False:#len(np.where(labels == 0 )[0]) > 0:\n",
    "    #    try:            \n",
    "            sil = metrics.silhouette_score(X,labels)\n",
    "        #print sil\n",
    "            return labels, sil\n",
    "        #except ValueError:\n",
    "            #print 'labels for bad sil: ', labels\n",
    "        #    return labels, -1.0\n",
    "    else:\n",
    "        return labels, -1.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# choose number of bootstrap samples\n",
    "def Permutation_test(data,nsims, other):\n",
    "    ndata_1 = len(data)\n",
    "    fakeidx1=np.floor(random.rand(ndata_1,nsims)*ndata_1)\n",
    "    #print 'len(fakeidx)', np.shape(fakeidx1)\n",
    "    fakeidx1 = fakeidx1.astype(int64)\n",
    "    fake1=data[fakeidx1]\n",
    "    #print 'fakedata', np.shape(fake1.T)\n",
    "    other = np.array([other]*len(fake1.T))\n",
    "    #print 'otherdata', np.shape(other)\n",
    "    return zip(fake1.T, other)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def Gaussian_Permutation_test(Xdata, nsims):    \n",
    "    pca = PCA(n_components=2)\n",
    "    pca.fit(Xdata)\n",
    "    C = pca.components_\n",
    "    ys = np.dot(C, Xdata.T)\n",
    "    \n",
    "    x = ys[0]\n",
    "    y = ys[1]\n",
    "    x_mean = np.mean(x)\n",
    "    x_std = np.std(x, ddof =1)\n",
    "    y_mean = np.mean(y)\n",
    "    y_std = np.std(y, ddof =1)\n",
    "    \n",
    "    x_ran = np.random.normal(x_mean, x_std, (len(x),nsims))\n",
    "    y_ran = np.random.normal(y_mean, y_std, (len(x), nsims))\n",
    "    \n",
    "    new = np.array(zip(x_ran,y_ran))\n",
    "\n",
    "    new_trans = []\n",
    "    for i in range(0,nsims):\n",
    "        new_trans.append(np.dot(np.linalg.inv(C), new[:,:,i].T))\n",
    "    return new_trans "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def p_value_K(sil, X):\n",
    "    sil_K = np.zeros(len(X))\n",
    "    for i in range(0,len(X)):\n",
    "        permed = np.array(zip(X[i][0], X[i][1]))\n",
    "        l, sil1 = K_mean(permed)\n",
    "        sil_K[i] = sil1\n",
    "    pval = (100. - stats.percentileofscore(sil_K,sil))/100.\n",
    "    #plt.hist(sil_K)\n",
    "    return pval\n",
    "\n",
    "def p_value_H(sil, X):\n",
    "    sil_K = np.zeros(len(X))\n",
    "    for i in range(0,len(X)):\n",
    "        permed = np.array(zip(X[i][0], X[i][1]))\n",
    "        l, sil1 = Hierarchical(permed)\n",
    "        sil_K[i] = sil1\n",
    "    pval = (100. - stats.percentileofscore(sil_K,sil))/100.\n",
    "    return pval\n",
    "\n",
    "def p_value_GMM(sil, X):\n",
    "    sil_GMM = np.zeros(len(X))\n",
    "    for i in range(0,len(X)):\n",
    "        permed = np.array(zip(X[i][0], X[i][1]))\n",
    "        l, sil1 = GMM(permed)\n",
    "        sil_GMM[i] = sil1\n",
    "    pval = (100. - stats.percentileofscore(sil_GMM,sil))/100.\n",
    "    return pval\n",
    "\n",
    "def p_value_DB(sil, X):\n",
    "    sil_K = np.zeros(len(X))\n",
    "    #print np.shape(X)\n",
    "    for i in range(0,len(X)):\n",
    "        permed = np.array(zip(X[i][0], X[i][1]))\n",
    "        l, sil1 = DBScan(permed)\n",
    "        sil_K[i] = sil1\n",
    "    pval = (100. - stats.percentileofscore(sil_K,sil))/100.\n",
    "    #plt.hist(sil_K)\n",
    "    #plt.vlines(sil, 0, 20)\n",
    "    return pval\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def beta(p):\n",
    "    b = len(np.array(np.where(p < 0.05))[0])\n",
    "    return (b+0.0)/len(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def dat_to_pval_K(N, SC, nsims):\n",
    "    dat = norm_data(N,SC)\n",
    "    labels, sil = K_mean(dat)\n",
    "    X2 = Permutation_test(dat[:,1], nsims, dat[:,0])\n",
    "    pval = p_value_K(sil, X2)\n",
    "    return pval\n",
    "\n",
    "def dat_to_pval_H(N, SC, nsims):\n",
    "    dat = norm_data(N,SC)\n",
    "    labels, sil = Hierarchical(dat)\n",
    "    X2 = Permutation_test(dat[:,1], nsims, dat[:,0])\n",
    "    pval = p_value_H(sil, X2)\n",
    "    return pval\n",
    "\n",
    "def dat_to_pval_GMM(N, SC, nsims):\n",
    "    dat = norm_data(N,SC)\n",
    "    labels, sil = GMM(dat)\n",
    "    X2 = Permutation_test(dat[:,1], nsims, dat[:,0])\n",
    "    pval = p_value_GMM(sil, X2)\n",
    "    return pval\n",
    "\n",
    "def dat_to_pval_DB(N, SC, nsims):\n",
    "    dat = norm_data(N,SC)\n",
    "    labels, sil = DBScan(dat)\n",
    "    #print 'sil',sil\n",
    "    #print 'labels', labels\n",
    "    X2 = Permutation_test(dat[:,1], nsims, dat[:,0])\n",
    "    pval = p_value_DB(sil, X2)\n",
    "    return pval\n",
    "\n",
    "def dat_to_pval_H_G(N, SC, nsims):\n",
    "    dat = norm_data(N,SC)\n",
    "    labels, sil = Hierarchical(dat)\n",
    "    #print 'sil',sil\n",
    "    #print 'labels', labels\n",
    "    X2 = Gaussian_Permutation_test(dat, nsims)\n",
    "    pval = p_value_H(sil, X2)\n",
    "    return pval"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def rev_dat_to_pval_K(N, SC, nsims):\n",
    "    dat = norm_data(N,SC)\n",
    "    labels, sil = K_mean(dat)\n",
    "    X2 = Permutation_test(dat[:,0], nsims, dat[:,1])\n",
    "    pval = p_value_K(sil, X2)\n",
    "    return pval\n",
    "\n",
    "def rev_dat_to_pval_H(N, SC, nsims):\n",
    "    dat = norm_data(N,SC)\n",
    "    labels, sil = Hierarchical(dat)\n",
    "    X2 = Permutation_test(dat[:,0], nsims, dat[:,1])\n",
    "    pval = p_value_H(sil, X2)\n",
    "    return pval\n",
    "\n",
    "def rev_dat_to_pval_GMM(N, SC, nsims):\n",
    "    dat = norm_data(N,SC)\n",
    "    labels, sil = GMM(dat)\n",
    "    X2 = Permutation_test(dat[:,0], nsims, dat[:,1])\n",
    "    pval = p_value_GMM(sil, X2)\n",
    "    return pval\n",
    "\n",
    "def rev_dat_to_pval_DB(N, SC, nsims):\n",
    "    dat = norm_data(N,SC)\n",
    "    labels, sil = DBScan(dat)\n",
    "    #print 'sil',sil\n",
    "    #print 'labels', labels\n",
    "    X2 = Permutation_test(dat[:,0], nsims, dat[:,1])\n",
    "    pval = p_value_DB(sil, X2)\n",
    "    return pval"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 92  96 100 104 108 112 116 120 124 128 132 136]\n"
     ]
    }
   ],
   "source": [
    "N = np.arange(92, 140, 4)# 120, 160, 10)\n",
    "print N"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "big_ps_GMM = []\n",
    "N = [100]\n",
    "bet = np.zeros(len(N))\n",
    "for j in range(0, len(N)):\n",
    "    ps = np.zeros(1000)#1000)\n",
    "    for i in range(0, 1000):#00):\n",
    "        ps[i] = dat_to_pval_GMM(N[j], 100, 1000)#00)\n",
    "    bet[j] = beta(ps)\n",
    "    big_ps_GMM.append(ps) \n",
    "    print bet\n",
    "    \n",
    "#np.savetxt('pval_GMM.txt', big_ps_GMM)\n",
    "#np.savetxt('beta_GMM.txt', bet)\n",
    "\n",
    "big_ps_GMM = np.zeros(len(N))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "big_ps_K = []\n",
    "bet = np.zeros(len(N))\n",
    "for j in range(0, len(N)):\n",
    "    ps = np.zeros(1000)\n",
    "    for i in range(0, 1000):\n",
    "        ps[i] = dat_to_pval_K(N[j], 100, 1000)\n",
    "    bet[j] = beta(ps)\n",
    "    big_ps_K.append(ps) \n",
    "    print bet\n",
    "    \n",
    "np.savetxt('pval_K.txt', big_ps_K)\n",
    "np.savetxt('beta_K.txt', bet)\n",
    "\n",
    "big_ps_K = np.zeros(len(N))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "rev_big_ps_K = []\n",
    "rev_bet = np.zeros(len(N))\n",
    "for j in range(0, len(N)):\n",
    "    rev_ps = np.zeros(1000)\n",
    "    for i in range(0, 1000):\n",
    "        rev_ps[i] = rev_dat_to_pval_K(N[j], 100, 1000)\n",
    "    rev_bet[j] = beta(rev_ps)\n",
    "    rev_big_ps_K.append(rev_ps) \n",
    "    print rev_bet\n",
    "    \n",
    "np.savetxt('rev_pval_K.txt', rev_big_ps_K)\n",
    "np.savetxt('rev_beta_K.txt', rev_bet)\n",
    "\n",
    "rev_big_ps_K = np.zeros(len(N))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#N = [160]\n",
    "#N = [70]\n",
    "big_ps_H = []\n",
    "bet_H = np.zeros(len(N))\n",
    "for j in range(0, len(N)):\n",
    "    if N[j] ==0:\n",
    "        bet_H[j] = 0.0\n",
    "    else:\n",
    "        ps1 = np.zeros(1000)\n",
    "        for i in range(0, 1000):\n",
    "            ps1[i] = dat_to_pval_H_G(N[j], 100, 1000)\n",
    "        bet_H[j] = beta(ps1)\n",
    "        big_ps_H.append(ps1) \n",
    "        print bet_H\n",
    "\n",
    "#np.savetxt('pval_Hier.txt', big_ps_H)\n",
    "#np.savetxt('beta_Hier.txt', bet_H)\n",
    "\n",
    "big_ps_H = np.zeros(len(N))\n",
    "#[ 0.89]   [ 0.857  0.873  0.875  0.896, [ 0.907]]  N = np.arange(120, 170, 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/karaponder/anaconda/lib/python2.7/site-packages/numpy/core/_methods.py:59: RuntimeWarning: Mean of empty slice.\n",
      "  warnings.warn(\"Mean of empty slice.\", RuntimeWarning)\n",
      "/Users/karaponder/anaconda/lib/python2.7/site-packages/numpy/core/_methods.py:71: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  ret = ret.dtype.type(ret / rcount)\n",
      "/Users/karaponder/anaconda/lib/python2.7/site-packages/IPython/kernel/__main__.py:14: RuntimeWarning: invalid value encountered in log10\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.896  0.     0.     0.     0.     0.     0.     0.     0.     0.     0.\n",
      "  0.   ]\n",
      "[ 0.896  0.897  0.     0.     0.     0.     0.     0.     0.     0.     0.\n",
      "  0.   ]\n",
      "[ 0.896  0.897  0.932  0.     0.     0.     0.     0.     0.     0.     0.\n",
      "  0.   ]\n",
      "[ 0.896  0.897  0.932  0.898  0.     0.     0.     0.     0.     0.     0.\n",
      "  0.   ]\n",
      "[ 0.896  0.897  0.932  0.898  0.912  0.     0.     0.     0.     0.     0.\n",
      "  0.   ]\n",
      "[ 0.896  0.897  0.932  0.898  0.912  0.921  0.     0.     0.     0.     0.\n",
      "  0.   ]\n",
      "[ 0.896  0.897  0.932  0.898  0.912  0.921  0.93   0.     0.     0.     0.\n",
      "  0.   ]\n",
      "[ 0.896  0.897  0.932  0.898  0.912  0.921  0.93   0.931  0.     0.     0.\n",
      "  0.   ]\n",
      "[ 0.896  0.897  0.932  0.898  0.912  0.921  0.93   0.931  0.934  0.     0.\n",
      "  0.   ]"
     ]
    }
   ],
   "source": [
    "rev_big_ps_H = []\n",
    "rev_bet_H = np.zeros(len(N))\n",
    "for j in range(0, len(N)):\n",
    "    rev_ps1 = np.zeros(1000)\n",
    "    for i in range(0, 1000):\n",
    "        rev_ps1[i] = rev_dat_to_pval_H(N[j], 100, 1000)\n",
    "    rev_bet_H[j] = beta(rev_ps1)\n",
    "    rev_big_ps_H.append(rev_ps1) \n",
    "    print rev_bet_H\n",
    "\n",
    "np.savetxt('rev_pval_Hier_150.txt', rev_big_ps_H)\n",
    "np.savetxt('rev_beta_Hier_150.txt', rev_bet_H)\n",
    "\n",
    "rev_big_ps_H = np.zeros(len(N))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "N=[70]\n",
    "big_ps_DB = []\n",
    "bet_DB = np.zeros(len(N))\n",
    "for j in range(0, len(N)):\n",
    "    ps2 = np.zeros(500)\n",
    "    for i in range(0, 500):\n",
    "        ps2[i] = dat_to_pval_DB(N[j], 100, 500)\n",
    "    bet_DB[j] = beta(ps2)\n",
    "    big_ps_DB.append(ps2) \n",
    "    print bet_DB\n",
    "    #print big_ps_DB\n",
    "    \n",
    "#np.savetxt('beta_DB.txt', bet_DB)\n",
    "#np.savetxt('pval_DB.txt', big_ps_DB)\n",
    "\n",
    "big_ps_DB = np.zeros(len(N))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "big_ps_DB = []\n",
    "bet_DB = np.zeros(len(N))\n",
    "for j in range(0, len(N)):\n",
    "    ps2 = np.zeros(1000)\n",
    "    for i in range(0, 1000):\n",
    "        ps2[i] = dat_to_pval_DB(N[j], 100, 1000)\n",
    "    bet_DB[j] = beta(ps2)\n",
    "    big_ps_DB.append(ps2) \n",
    "    print bet_DB\n",
    "    \n",
    "np.savetxt('beta_DB1.txt', bet_DB)\n",
    "np.savetxt('pval_DB1.txt', big_ps_DB)\n",
    "\n",
    "big_ps_DB = np.zeros(len(N))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "rev_big_ps_DB = []\n",
    "rev_bet_DB = np.zeros(len(N))\n",
    "for j in range(0, len(N)):\n",
    "    rev_ps2 = np.zeros(1000)\n",
    "    for i in range(0, 1000):\n",
    "        rev_ps2[i] = rev_dat_to_pval_DB(N[j], 100, 1000)\n",
    "    rev_bet_DB[j] = beta(rev_ps2)\n",
    "    rev_big_ps_DB.append(rev_ps2) \n",
    "    print rev_bet_DB\n",
    "    \n",
    "np.savetxt('rev_beta_DB.txt', rev_bet_DB)\n",
    "np.savetxt('rev_pval_DB.txt', rev_big_ps_DB)\n",
    "\n",
    "rev_big_ps_DB = np.zeros(len(N))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "rev_big_ps_DB = []\n",
    "rev_bet_DB = np.zeros(len(N))\n",
    "for j in range(0, len(N)):\n",
    "    rev_ps2 = np.zeros(1000)\n",
    "    for i in range(0, 1000):\n",
    "        rev_ps2[i] = rev_dat_to_pval_DB(N[j], 100, 1000)\n",
    "    rev_bet_DB[j] = beta(rev_ps2)\n",
    "    rev_big_ps_DB.append(rev_ps2) \n",
    "    print rev_bet_DB\n",
    "    \n",
    "np.savetxt('rev_beta_DB1.txt', rev_bet_DB)\n",
    "np.savetxt('rev_pval_DB1.txt', rev_big_ps_DB)\n",
    "\n",
    "rev_big_ps_DB = np.zeros(len(N))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\n",
    "## to plot\n",
    "\n",
    "# plt.scatter(X[:, 0], X[:, 1], c=labels, cmap='autumn_r',s=50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#one, two = Permutation_test(dat[:,0], 20)\n",
    "#print np.shape(X2)#np.array([dat[:,1]]*len(X2)))\n",
    "#print len(X2)\n",
    "#print np.shape(one), np.shape(two), np.shape(dat)\n",
    "#print np.shape(one[:,0])\n",
    "#print np.shape(np.row_stack((one,two)))\n",
    "#comb = np.row_stack((one,two))\n",
    "#print np.shape(comb[0])\n",
    "#print np.shape(comb.T)\n",
    "#print comb.T\n",
    "\n",
    "#print 'pvalue:', (100.-stats.percentileofscore(sil_K,sil))/100.\n",
    "#print sil\n",
    "#sil = K_mean(zip(rrr[14][0], rrr[14][1]))\n",
    "\n",
    "#print sil_K\n",
    "#plt.hist(sil_K)\n",
    "#plt.vlines(sil, 0, 300)\n",
    "\n",
    "#print len(dat[:,0])#len(np.array(np.where(labels == 2))[0])\n",
    "\n",
    "\n",
    "#dat_to_pval_K(10, 100, 1000)\n",
    "#X = norm_data(10,100)\n",
    "#print X\n",
    "#sil = DBScan(X)\n",
    "#X2 = Permutation_test(X[:,1],10, X[:,0])\n",
    "#p_value_DB(sil, X2)\n",
    "#print np.array(zip(X2[1][0], X2[1][1]))\n",
    "\n",
    "#dat_to_pval_DB(100, 100, 1000)\n",
    "#dat = norm_data(0,100)\n",
    "#labels, sil = DBScan(dat)\n",
    "#print sil , '\\n'\n",
    "#print labels, '\\n'\n",
    "#plt.figure()\n",
    "#plt.scatter(dat[:, 0], dat[:, 1], c=labels, cmap='autumn_r',s=50)\n",
    "#plt.figure()\n",
    "#dat1 = norm_data(100,100)\n",
    "#labels1, sil1 = DBScan(dat1)\n",
    "#plt.scatter(dat1[:, 0], dat1[:, 1], c=labels1, cmap='autumn_r',s=50)\n",
    "#print labels1, '\\n', sil1"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}