Re-updated Ch5 and Ch6 to Python3

Author: mmargenot

Chapter5_LossFunctions/Chapter5.ipynb

@@ -58,7 +58,7 @@ def calc_delta_r(x_predicted,y_predicted,x_true,y_true):
     for  perm in it.permutations(a[num_halos-2],num_halos):
         which_true_halos=[]
         which_predicted_halos=[]
-        for j in xrange(num_halos): #loop through all the true halos with the
+        for j in range(num_halos): #loop through all the true halos with the
 
             distances_perm[count,j]=np.sqrt((x_true[j]-x_predicted[int(perm[j])])**2\
                                       +(y_true[j]-y_predicted[int(perm[j])])**2)
@@ -141,7 +141,7 @@ def convert_to_360(angle, x_in, y_in):
         theta: the angle in the range 0:2pi
     """
     n = len(x_in)
-    for i in xrange(n):
+    for i in range(n):
         if x_in[i] < 0 and y_in[i] > 0:
             angle[i] = angle[i]+mt.pi
         elif x_in[i] < 0 and y_in[i] < 0:
@@ -204,7 +204,7 @@ def main_score( nhalo_all, x_true_all, y_true_all, x_ref_all, y_ref_all, sky_pre
 
         x_predicted=np.array([],dtype=float)
         y_predicted=np.array([],dtype=float)
-        for i in xrange(nhalo):
+        for i in range(nhalo):
             x_predicted=np.append(x_predicted,float(sky[0])) #get the predicted values
             y_predicted=np.append(y_predicted,float(sky[1]))
             #The solution file for the test data provides masses 
@@ -271,9 +271,9 @@ def main_score( nhalo_all, x_true_all, y_true_all, x_ref_all, y_ref_all, sky_pre
     W1=1./1000. #Weight the av_r such that < 1 is a good score > 1 is not so good.
     W2=1.
     metric = W1*av_r + W2*angle_vec #Weighted metric, weights TBD
-    print 'Your average distance in pixels you are away from the true halo is', av_r
-    print 'Your average angular vector is', angle_vec
-    print 'Your score for the training data is', metric
+    print('Your average distance in pixels you are away from the true halo is', av_r)
+    print('Your average angular vector is', angle_vec)
+    print('Your score for the training data is', metric)
     return metric
 
 
@@ -316,10 +316,9 @@ def main(user_fname, fname):
                                                          #first input would be
                                                          #a float, if succeed it
                                                          #is not a header
-        print 'THE INPUT FILE DOES NOT APPEAR TO HAVE A HEADER'
+        print('THE INPUT FILE DOES NOT APPEAR TO HAVE A HEADER')
     except :
-        print 'THE INPUT FILE APPEARS TO HAVE A HEADER, SKIPPING THE FIRST LINE'
-
+        print('THE INPUT FILE APPEARS TO HAVE A HEADER, SKIPPING THE FIRST LINE')
         skip_header = sky_prediction.next()
 
 
@@ -331,7 +330,7 @@ def main(user_fname, fname):
         if does_it_exist > 0: #If it does then find the matching solutions to the sky_id
                             selectskyinsolutions=true_sky_id.index(sky_id)-1
         else: #Otherwise exit
-            print 'Sky_id does not exist, formatting problem: ',sky_id
+            print('Sky_id does not exist, formatting problem: ',sky_id)
             sys.exit(2)
 
 
@@ -342,7 +341,7 @@ def main(user_fname, fname):
 
         x_predicted=np.array([],dtype=float)
         y_predicted=np.array([],dtype=float)
-        for i in xrange(nhalo):
+        for i in range(nhalo):
             x_predicted=np.append(x_predicted,float(sky[2*i+1])) #get the predicted values
             y_predicted=np.append(y_predicted,float(sky[2*i+2]))
             #The solution file for the test data provides masses 
@@ -409,9 +408,9 @@ def main(user_fname, fname):
     W1=1./1000. #Weight the av_r such that < 1 is a good score > 1 is not so good.
     W2=1.
     metric = W1*av_r + W2*angle_vec #Weighted metric, weights TBD
-    print 'Your average distance in pixels you are away from the true halo is', av_r
-    print 'Your average angular vector is', angle_vec
-    print 'Your score for the training data is', metric
+    print('Your average distance in pixels you are away from the true halo is', av_r)
+    print('Your average angular vector is', angle_vec)
+    print('Your score for the training data is', metric)
 
 
 if __name__ == "__main__":

Chapter5_LossFunctions/DarkWorldsMetric.py

@@ -2,14 +2,14 @@
 from matplotlib.patches import Ellipse
 import numpy as np
 
-def draw_sky( galaxies ):
+def draw_sky(galaxies):
     """adapted from Vishal Goklani"""
     size_multiplier = 45
     fig = plt.figure(figsize=(10,10))
     #fig.patch.set_facecolor("blue")
     ax = fig.add_subplot(111, aspect='equal')
     n = galaxies.shape[0]
-    for i in xrange(n):
+    for i in range(n):
         _g = galaxies[i,:]
         x,y = _g[0], _g[1]
         d = np.sqrt( _g[2]**2 + _g[3]**2 )

Chapter5_LossFunctions/draw_sky2.py

@@ -11,7 +11,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "#Chapter 6\n",
+    "# Chapter 6\n",
     "\n",
     "____\n",
     "\n",
@@ -27,7 +27,7 @@
     "\n",
     "Up until now, we have mostly ignored our choice of priors. This is unfortunate as we can be very expressive with our priors, but we also must be careful about choosing them. This is especially true if we want to be objective, that is, not to express any personal beliefs in the priors. \n",
     "\n",
-    "###Subjective vs Objective priors\n",
+    "### Subjective vs Objective priors\n",
     "\n",
     "Bayesian priors can be classified into two classes: *objective* priors, which aim to allow the data to influence the posterior the most, and *subjective* priors, which allow the practitioner to express his or her views into the prior. \n",
     "\n",
@@ -926,7 +926,7 @@
     "figsize(12.0, 8)\n",
     "beta = stats.beta\n",
     "hidden_prob = beta.rvs(1, 13, size=35)\n",
-    "print hidden_prob\n",
+    "print(hidden_prob)\n",
     "bandits = Bandits(hidden_prob)\n",
     "bayesian_strat = BayesianStrategy(bandits)\n",
     "\n",
@@ -1042,8 +1042,8 @@
     "                       \"AMZN\": (0.03, 0.02),\n",
     "                       }\n",
     "\n",
-    "for i, (name, params) in enumerate(expert_prior_params.iteritems()):\n",
-    "    plt.subplot(2, 2, i)\n",
+    "for i, (name, params) in enumerate(expert_prior_params.items()):\n",
+    "    plt.subplot(2, 2, i + 1)\n",
     "    y = normal.pdf(x, params[0], scale=params[1])\n",
     "    #plt.plot( x, y, c = colors[i] )\n",
     "    plt.fill_between(x, 0, y, color=colors[i], linewidth=2,\n",
@@ -1103,7 +1103,7 @@
     "\n",
     "stocks = [\"AAPL\", \"GOOG\", \"TSLA\", \"AMZN\"]\n",
     "\n",
-    "enddate = datetime.datetime.now().strftime(\"%Y-%m-%d\")  # today's date.\n",
+    "enddate = \"2015-04-27\"\n",
     "startdate = \"2012-09-01\"\n",
     "\n",
     "stock_closes = {}\n",
@@ -1120,7 +1120,7 @@
     "    _previous_day = np.roll(stock_closes[stock], -1)\n",
     "    stock_returns[stock] = ((stock_closes[stock] - _previous_day) / _previous_day)[:n_observations]\n",
     "\n",
-    "dates = map(lambda x: datetime.datetime.strptime(x, \"%Y-%m-%d\"), x[1:n_observations + 1, 0])"
+    "dates = list(map(lambda x: datetime.datetime.strptime(x, \"%Y-%m-%d\"), x[1:n_observations + 1, 0]))"
    ]
   },
   {
@@ -1144,12 +1144,12 @@
    "source": [
     "figsize(12.5, 4)\n",
     "\n",
-    "for _stock, _returns in stock_returns.iteritems():\n",
+    "for _stock, _returns in stock_returns.items():\n",
     "    p = plt.plot((1 + _returns)[::-1].cumprod() - 1, '-o', label=\"%s\" % _stock,\n",
     "                 markersize=4, markeredgecolor=\"none\")\n",
     "\n",
     "plt.xticks(np.arange(100)[::-8],\n",
-    "           map(lambda x: datetime.datetime.strftime(x, \"%Y-%m-%d\"), dates[::8]),\n",
+    "           list(map(lambda x: datetime.datetime.strftime(x, \"%Y-%m-%d\"), dates[::8])),\n",
     "           rotation=60);\n",
     "\n",
     "plt.legend(loc=\"upper left\")\n",
@@ -1179,9 +1179,9 @@
     "figsize(11., 5)\n",
     "returns = np.zeros((n_observations, 4))\n",
     "\n",
-    "for i, (_stock, _returns) in enumerate(stock_returns.iteritems()):\n",
+    "for i, (_stock, _returns) in enumerate(stock_returns.items()):\n",
     "    returns[:, i] = _returns\n",
-    "    plt.subplot(2, 2, i)\n",
+    "    plt.subplot(2, 2, i+1)\n",
     "    plt.hist(_returns, bins=20,\n",
     "             normed=True, histtype=\"stepfilled\",\n",
     "             color=colors[i], alpha=0.7)\n",
@@ -1258,7 +1258,7 @@
     "for i in range(4):\n",
     "    plt.hist(mu_samples[:, i], alpha=0.8 - 0.05 * i, bins=30,\n",
     "             histtype=\"stepfilled\", normed=True,\n",
-    "             label=\"%s\" % stock_returns.keys()[i])\n",
+    "             label=\"%s\" % list(stock_returns.keys())[i])\n",
     "\n",
     "plt.vlines(mu_samples.mean(axis=0), 0, 500, linestyle=\"--\", linewidth=.5)\n",
     "\n",
@@ -1302,8 +1302,8 @@
     "    plt.subplot(2, 2, i + 1)\n",
     "    plt.hist(mu_samples[:, i], alpha=0.8 - 0.05 * i, bins=30,\n",
     "             histtype=\"stepfilled\", normed=True, color=colors[i],\n",
-    "             label=\"%s\" % stock_returns.keys()[i])\n",
-    "    plt.title(\"%s\" % stock_returns.keys()[i])\n",
+    "             label=\"%s\" % list(stock_returns.keys())[i])\n",
+    "    plt.title(\"%s\" % list(stock_returns.keys())[i])\n",
     "    plt.xlim(-0.15, 0.15)\n",
     "\n",
     "plt.suptitle(\"Posterior distribution of daily stock returns\")\n",
@@ -2512,21 +2512,21 @@
  "metadata": {
   "anaconda-cloud": {},
   "kernelspec": {
-   "display_name": "Python 2",
+   "display_name": "Python [conda env:bayes]",
    "language": "python",
-   "name": "python2"
+   "name": "conda-env-bayes-py"
   },
   "language_info": {
    "codemirror_mode": {
     "name": "ipython",
-    "version": 2
+    "version": 3
    },
    "file_extension": ".py",
    "mimetype": "text/x-python",
    "name": "python",
    "nbconvert_exporter": "python",
-   "pygments_lexer": "ipython2",
-   "version": "2.7.10"
+   "pygments_lexer": "ipython3",
+   "version": "3.5.2"
   }
  },
  "nbformat": 4,

Chapter6_Priorities/Chapter6.ipynb

@@ -3,13 +3,13 @@
 
 import scipy.stats as stats
 import numpy as np
-from pymc import rbeta
+#from pymc import rbeta
 
 rand = np.random.rand
 beta = stats.beta
 
 
-class GeneralBanditStrat( object ):	
+class GeneralBanditStrat(object):	
 
     """
     Implements a online, learning strategy to solve
@@ -32,72 +32,72 @@ class GeneralBanditStrat( object ):
     def __init__(self, bandits, choice_function):
 
         self.bandits = bandits
-        n_bandits = len( self.bandits )
-        self.wins = np.zeros( n_bandits )
-        self.trials = np.zeros(n_bandits )
+        n_bandits = len(self.bandits)
+        self.wins = np.zeros(n_bandits)
+        self.trials = np.zeros(n_bandits)
         self.N = 0
         self.choices = []
         self.score = []
         self.choice_function = choice_function
 
-    def sample_bandits( self, n=1 ):
+    def sample_bandits(self, n=1):
 
-        score = np.zeros( n )
-        choices = np.zeros( n )
+        score = np.zeros(n)
+        choices = np.zeros(n)
 
         for k in range(n):
             #sample from the bandits's priors, and select the largest sample
             choice = self.choice_function(self)
 
             #sample the chosen bandit
-            result = self.bandits.pull( choice )
+            result = self.bandits.pull(choice)
 
             #update priors and score
-            self.wins[ choice ] += result
-            self.trials[ choice ] += 1
-            score[ k ] = result 
+            self.wins[choice] += result
+            self.trials[choice] += 1
+            score[k] = result 
             self.N += 1
-            choices[ k ] = choice
+            choices[k] = choice
 
-        self.score = np.r_[ self.score, score ]
-        self.choices = np.r_[ self.choices, choices ]
+        self.score = np.r_[self.score, score]
+        self.choices = np.r_[self.choices, choices]
         return 
 
 
 def bayesian_bandit_choice(self):
-	return np.argmax( rbeta( 1 + self.wins, 1 + self.trials - self.wins) )
+	return np.argmax(np.random.beta(1 + self.wins, 1 + self.trials - self.wins))
 
-def max_mean( self ):
+def max_mean(self):
     """pick the bandit with the current best observed proportion of winning """
-    return np.argmax( self.wins / ( self.trials +1 ) )
+    return np.argmax(self.wins / (self.trials +1))
 
 def lower_credible_choice( self ):
     """pick the bandit with the best LOWER BOUND. See chapter 5"""
     def lb(a,b):
-        return a/(a+b) - 1.65*np.sqrt( (a*b)/( (a+b)**2*(a+b+1) ) )
+        return a/(a+b) - 1.65*np.sqrt((a*b)/( (a+b)**2*(a+b+1)))
     a = self.wins + 1
     b = self.trials - self.wins + 1
-    return np.argmax( lb(a,b) )
+    return np.argmax(lb(a,b))
 
-def upper_credible_choice( self ):
+def upper_credible_choice(self):
     """pick the bandit with the best LOWER BOUND. See chapter 5"""
     def lb(a,b):
-        return a/(a+b) + 1.65*np.sqrt( (a*b)/( (a+b)**2*(a+b+1) ) )
+        return a/(a+b) + 1.65*np.sqrt((a*b)/((a+b)**2*(a+b+1)))
     a = self.wins + 1
     b = self.trials - self.wins + 1
-    return np.argmax( lb(a,b) )
+    return np.argmax(lb(a,b))
 
-def random_choice( self):
-    return np.random.randint( 0, len( self.wins ) )
+def random_choice(self):
+    return np.random.randint(0, len(self.wins))
 
 
-def ucb_bayes( self ):
+def ucb_bayes(self):
 	C = 0
 	n = 10000
-	alpha =1 - 1./( (self.N+1) )
-	return np.argmax( beta.ppf( alpha,
+	alpha =1 - 1./((self.N+1))
+	return np.argmax(beta.ppf(alpha,
 							   1 + self.wins, 
-							   1 + self.trials - self.wins ) )
+							   1 + self.trials - self.wins))
 
 
 
@@ -117,7 +117,7 @@ def __init__(self, p_array):
         self.p = p_array
         self.optimal = np.argmax(p_array)
 
-    def pull( self, i ):
+    def pull(self, i):
         #i is which arm to pull
         return rand() < self.p[i]