Update 2D uncertainties and intrinsic variance exercises

straightline/src/data_allerr.dat

@@ -0,0 +1,21 @@
+# index & x & y & sigm_y & sigma_x & corr 
+1 & 201 & 592 & 61 & 9 & -0.84
+2 & 244 & 401 & 25 & 4 & 0.31
+3 & 47 & 583 & 38 & 11 & 0.64
+4 & 287 & 402 & 15 & 7 & -0.27
+5 & 203 & 495 & 21 & 5 & -0.33
+6 & 58 & 173 & 15 & 9 & 0.67
+7 & 210 & 479 & 27 & 4 & -0.02
+8 & 202 & 504 & 14 & 4 & -0.05
+9 & 198 & 510 & 30 & 11 & -0.84
+10 & 158 & 416 & 16 & 7 & -0.69
+11 & 165 & 393 & 14 & 5 & 0.30
+12 & 201 & 442 & 25 & 5 & -0.46
+13 & 157 & 317 & 52 & 5 & -0.03
+14 & 131 & 311 & 16 & 6 & 0.50
+15 & 166 & 400 & 34 & 6 & 0.73
+16 & 160 & 337 & 31 & 5 & -0.52
+17 & 186 & 423 & 42 & 9 & 0.90
+18 & 125 & 334 & 26 & 8 & 0.40
+19 & 218 & 533 & 16 & 6 & -0.78
+20 & 146 & 344 & 22 & 5 & -0.56

straightline/src/data_yerr.dat

@@ -0,0 +1,21 @@
+# index & x & y & sigm_y
+1 & 201 & 592 & 61
+2 & 244 & 401 & 25
+3 & 47 & 583 & 38
+4 & 287 & 402 & 15
+5 & 203 & 495 & 21
+6 & 58 & 173 & 15
+7 & 210 & 479 & 27
+8 & 202 & 504 & 14
+9 & 198 & 510 & 30
+10 & 158 & 416 & 16
+11 & 165 & 393 & 14
+12 & 201 & 442 & 25
+13 & 157 & 317 & 52
+14 & 131 & 311 & 16
+15 & 166 & 400 & 34
+16 & 160 & 337 & 31
+17 & 186 & 423 & 42
+18 & 125 & 334 & 26
+19 & 218 & 533 & 16
+20 & 146 & 344 & 22

straightline/src/ex12.py

@@ -97,17 +97,17 @@ def ex12(exclude=sc.array([1,2,3,4]),plotfilename='ex12.png',
                                        (Z,ycovar),disp=False)
     if linalg.norm(bestfit2d-bestfit2d1) > 10**-12:
         if linalg.norm(bestfit2d-bestfit2d1) < 10**-6:
-            print "Different optimizers give slightly different results..."
+            print("Different optimizers give slightly different results...")
         else:
-            print "Different optimizers give rather different results..."
-        print "The norm of the results differs by %g" % linalg.norm(bestfit2d-bestfit2d1)
+            print("Different optimizers give rather different results...")
+        print("The norm of the results differs by %g" % linalg.norm(bestfit2d-bestfit2d1))
         try:
             x=raw_input('continue to plot? [yn]\n')
         except EOFError:
-            print "Since you are in non-interactive mode I will assume 'y'"
+            print("Since you are in non-interactive mode I will assume 'y'")
             x='y'
         if x == 'n':
-            print "returning..."
+            print("returning...")
             return -1
 
     #Plot result
@@ -175,5 +175,5 @@ def objective(mb,Z,ycovar):
     cost= v[1]
     delta= sc.dot(v,Z.T)-mb[0]*cost
     sigma2= sc.dot(v,sc.dot(ycovar,v))
-    return 0.5*sc.sum(delta**2./sigma2+sc.log(sigma2))
+    return 0.5*sc.sum(delta**2./sigma2+sc.log(sigma2)+sc.log(1.+mb[1]**2.))
 

straightline/src/ex13.py

@@ -168,7 +168,7 @@ def ex13(exclude=sc.array([1,2,3,4]),plotfilename='ex13.png',
                    xlabel=r'$x$',ylabel=r'$y$',zorder=2)
     for ii in range(10):
         #Random sample
-        ransample= sc.floor((stats.uniform.rvs()*nsamples))
+        ransample= sc.floor((stats.uniform.rvs()*nsamples)).astype('int')
         ransample= samples.T[ransample,0:2]
         mf= m.sqrt(1./m.cos(ransample[1])**2.-1.)
         b= ransample[0]/m.cos(ransample[1])
@@ -263,7 +263,7 @@ def objective(pars,Z,ycovar):
     for ii in range(ndata):
         detVycovar[ii]= m.sqrt(linalg.det(V+ycovar[:,ii,:]))
         deltaOUT[ii]= sc.dot(Z[ii,:]-Zb,sc.dot(linalg.inv(V+ycovar[:,ii,:]),Z[ii,:]-Zb))
-    return sc.sum(sc.log((1.-Pb)/sc.sqrt(2.*m.pi*sigma2)*
+    return sc.sum(sc.log((1.-Pb)/sc.sqrt(2.*m.pi*sigma2/sc.cos(t)**2.)*
                          sc.exp(-0.5*delta**2./sigma2)
                          +Pb/2./m.pi/detVycovar
                          *sc.exp(-0.5*deltaOUT)))
@@ -287,4 +287,4 @@ def Pbad(pars,Z,ycovar):
     detVycovar= m.sqrt(linalg.det(V+ycovar[:,:]))
     deltaOUT= sc.dot(Z-Zb,sc.dot(linalg.inv(V+ycovar[:,:]),Z-Zb))
     Pbad= Pb/2./m.pi/detVycovar*sc.exp(-0.5*deltaOUT)
-    return Pbad/(Pbad+(1.-Pb)/sc.sqrt(2.*m.pi*sigma2)*sc.exp(-0.5*delta**2./sigma2))
+    return Pbad/(Pbad+(1.-Pb)/sc.sqrt(2.*m.pi*sigma2)*sc.exp(-0.5*delta**2./sigma2/sc.cos(t)**2.))

straightline/src/ex16.py

@@ -90,7 +90,7 @@ def ex16(exclude=sc.array([3]),plotfilename='ex16.png',
     bestfitvar= linalg.inv(bestfitvar)
     bestfit= sc.dot(bestfitvar,bestfit)
     #Now optimize
-    initial= sc.array([bestfit[0],bestfit[1],sc.log(100)])
+    initial= sc.array([bestfit[0],bestfit[1],sc.log(1000)])
     bestfit2d1= optimize.fmin(objective,initial,(Z,ycovar),disp=False)
     #Restart the optimization once using a different method
     bestfit2d= optimize.fmin_powell(objective,initial,
@@ -113,6 +113,7 @@ def ex16(exclude=sc.array([3]),plotfilename='ex16.png',
     b= bestfit2d[0]
     mf= bestfit2d[1]
     V=sc.exp(bestfit2d[2]/2.)
+    print(b,mf,V)
     cost= 1./sc.sqrt(1+mf**2.)
     bcost= b*cost
 
@@ -182,4 +183,4 @@ def objective(mbV,Z,ycovar):
     cost= v[1]
     delta= sc.dot(v,Z.T)-mbV[0]*cost
     sigma2= sc.dot(v,sc.dot(ycovar,v))+V
-    return 0.5*sc.sum(delta**2./sigma2+sc.log(sigma2))
+    return 0.5*sc.sum(delta**2./sigma2+sc.log(sigma2)+sc.log(1.+mbV[1]**2.))

straightline/src/ex17.py

@@ -156,7 +156,7 @@ def ex17(exclude=sc.array([3]),plotfilename='ex17.png',
     #Plot
     plot.bovy_print(**bovyprintargs)
     hist, bins, patchess= plot.bovy_hist(sc.exp(samples.T[:,2]/2.),edgecolor='k',
-                                      bins=round(sc.sqrt(nsamples)/2.),
+                                      bins=int(round(sc.sqrt(nsamples)/2.)),
                                       xlabel=r'$\sqrt{V}$',normed=True,
                                       histtype='step')
     cumhist= sc.cumsum(hist)/sc.sum(hist)/(bins[1]-bins[0])
@@ -265,5 +265,5 @@ def objective(pars,Z,ycovar):
     v= sc.array([-sc.sin(t),sc.cos(t)])
     delta= sc.dot(v,Z.T)-bcost
     sigma2= sc.dot(v,sc.dot(ycovar,v))+V
-    return -0.5*sc.sum(delta**2./sigma2+sc.log(sigma2))#+pars[2]
+    return -0.5*sc.sum(delta**2./sigma2+sc.log(sigma2/sc.cos(t)**2.))
 

straightline/src/generate_data.py

@@ -359,7 +359,7 @@ def read_data(datafilename='data_yerr.dat',allerr=False):
             rawdata.append(float(i.group()))
             nvalue+= 1
         if nvalue != ncol:
-            print "Warning, number of columns for this record does not match the expected number"
+            print("Warning, number of columns for this record does not match the expected number")
     #Now process the raw data
     out=[]
     for ii in range(nline):