<jupyter_start><jupyter_text>Numerical Methods  -- Assignment 6   Problem 1 -- Minimisation  (a) Brent's method  In numerical analysis, Brent's method is a root-finding algorithm combining the bisection method, the secant method and inverse quadratic interpolation. It has the reliability of bisection but it can be as quick as some of the less-reliable methods. Now given the following function:\begin{equation}f(x) = \frac{(x+3)^2}{2}+\frac{x-5}{3}-\frac{10}{4}\end{equation}<jupyter_code>%config InlineBackend.figure_format = 'retina' 
from scipy.optimize import brent
import numpy as np
import matplotlib.pyplot as plt

f = lambda x: 0.5*np.power(x+3,2)+(x-5)/3-10./4
minimum = brent(f,brack=(-100,100))
x = np.linspace(-100,100,100)
y = f(x)
print "The minimum of the function is found to be:",minimum
#Plot to check the shape of the function
plt.plot(x,y,'-')
plt.legend('f(x)')
plt.show()<jupyter_output>The minimum of the function is found to be: -3.33333338268<jupyter_text>(b) Conjugate gradient method $$f(x,y) = 0.26((x+0.5)^2+(y-1)^2)-0.48xy$$<jupyter_code>from scipy.optimize import fmin_cg,minimize

def func(x):
    return 0.26*((x[0]+0.5)**2+(x[1]-1)**2)-0.48*x[0]*x[1]

def fprime(x):
    return np.array((0.26*2*(x[0]+0.5)-0.48*x[1],0.26*2*(x[1]-1)-0.48*x[0]))
fc = lambda x: 0.26*(np.power(x[0]+0.5,2)+np.power(x[1]-1,2))-0.48*x[0]*x[1]
fp = lambda x: np.array((0.26*2*(x[0]+0.5)-0.48*x[1],0.26*2*(x[1]-1)-0.48*x[0]))
x0 = [2,5]
fmin = minimize(fc,x0,jac=fp,method='CG')
print "The minimum of the function is found at (x,y)=:",fmin.x
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.pyplot import figure, show
import numpy as np

x = np.linspace(-100,100,100)
y = x
xx,yy = np.meshgrid(x,y)
X = [xx,yy]
zz = func(X)
#plot
fig = plt.figure(figsize=(10,6))
frame = fig.add_subplot(1,1,1,projection ='3d')
surf = frame.plot_surface(xx,yy,zz,cmap='Dark2')
frame.set_xlabel('x')
frame.set_ylabel('y')
frame.set_zlabel('z')
frame.set_title('f(x,y)')
fig.colorbar(surf,ax=frame)
show()<jupyter_output><empty_output><jupyter_text>Problem 2 -- x,y data The $\textit{Maximum Likelihood}$ is done by minimizing the log-likelihood function.$$ln\,L = -\frac{1}{2}\sum_{i=1}^{N}\left(\frac{(data-model)^2}{2\sigma^2}+ln\,(2\pi\sigma^2)\right)$$,where the model is described by $y = kx+m$, and $\sigma^2 = yerr^2+f^2(kx+m)^2$. First use $\chi^2$ fitting and reduced $\chi^2$ to get a reference of the values of parameters.<jupyter_code>from kapteyn import kmpfit

data = np.genfromtxt('Downloads/table_1.dat')
x = data[:,0]
y = data[:,1]
ey = data[:,2]
#define a model function, since assuming linear correlation
def model(params, x):
    a, b = params 
    return a + b*x

# A residuals function
def residuals(p, data):  # Residuals function needed by kmpfit
    x, y, err= data
    return (y - model(p,x))/err


p0 = (1, 1)
fitobj = kmpfit.Fitter(residuals=residuals, data=(x,y,ey))
status = fitobj.fit(params0=p0)
print("Was fit successful?", bool(status))

a, b = fitobj.params
da, db = fitobj.stderr
s = u"Model parameters a={} \u00B1 {} and b={} \u00B1 {}".format(a, da, b, db)
print(s)
print("Chi squared, reduced Chi squared:", fitobj.chi2_min, fitobj.rchi2_min)
print("The errors derived from the covariance matrix:", fitobj.xerror)


from matplotlib.pyplot import figure, show

fig = figure()
frame = fig.add_subplot(1,1,1)
frame.errorbar(x, y, yerr=ey,marker='o', ls='')
xmin = x.min(); xmax = x.max()
p = (a,b)
ymin = model(p,xmin); ymax = model(p,xmax)
frame.plot([xmin, xmax], [ymin, ymax])
frame.set_xlabel('Distance (Mpc)')
frame.set_ylabel('Velocity (km/s)')
show()<jupyter_output><empty_output><jupyter_text>(a) -- Optimal parameters<jupyter_code>from math import *
import time

def log_like(params,x,y,ey):
    k,m,f = params
    sig = np.sqrt(ey**2+(f*(k*x+m))**2)
    L = ((k*x+m)-y)**2/(sig**2)+np.log(2*pi*sig**2)
    return 0.5*np.sum(L)
#minimize via three methods
t1 = time.time()
l_BFGS = minimize(log_like,x0=[0,0,0],args=(x,y,ey),method='BFGS')
t2 = time.time()
t3 = time.time()
l_NM = minimize(log_like,x0=[0,0,0],args=(x,y,ey),method='Nelder-Mead')
t4 = time.time()
t5 = time.time()
l_CG = minimize(log_like,x0=[0,0,0],args=(x,y,ey),method='CG')
t6 = time.time()
print '%6s %6s %6s %6s %12s %12s' % ("Methods","k","m","f","time","Success")
print '%6s %.7f %.7f %.7f %.6f %8s' % ("BFGS",l_BFGS.x[0],l_BFGS.x[1],l_BFGS.x[2],t2-t1,l_BFGS.success)
print '%6s %.7f %.7f %.7f %.6f %8s' % ("NM",l_NM.x[0],l_NM.x[1],np.abs(l_NM.x[2]),t4-t3,l_NM.success)
print '%6s %.7f %.7f %.7f %.6f %8s' % ("CG",l_CG.x[0],l_CG.x[1],l_CG.x[2],t6-t5,l_CG.success)<jupyter_output>Methods      k      m      f         time      Success
  BFGS 0.6850430 3.3064525 0.2998975 0.014389     True
    NM 0.6850535 3.3064455 0.2998962 0.021051     True
    CG 0.6850429 3.3064530 0.2998975 0.013389     True<jupyter_text>The table printed above listed the optimization results via the three methods. In this case CG method performs the fastest. All three methods uses x,y and error in y to minimize the log-likelihood function.  (b) -- Plot fitting<jupyter_code>fig = figure()
frame = fig.add_subplot(1,1,1)
frame.errorbar(x, y, yerr=ey,marker='*', fillstyle='none',ls='',color='indigo')
xmin = x.min(); xmax = x.max()
p1 = [l_BFGS.x[1],l_BFGS.x[0]]
p2 = [l_NM.x[1],l_NM.x[0]]
p3 = [l_CG.x[1],l_CG.x[0]]
ymin1 = model(p1,xmin); ymax1 = model(p1,xmax)
ymin2 = model(p2,xmin); ymax2 = model(p2,xmax)
ymin3 = model(p3,xmin); ymax3 = model(p3,xmax)
frame.plot([xmin, xmax], [ymin1, ymax1])
frame.plot([xmin, xmax], [ymin1, ymax1])
frame.plot([xmin, xmax], [ymin1, ymax1])
frame.set_xlabel('x')
frame.set_ylabel('y')
frame.legend(['BFGS','NM','CG','data'])
show()<jupyter_output><empty_output><jupyter_text>Notice that because of small deviations the three fits overlay on each other and is not distinguashiable on the graph above.  Problem 3 -- Dispersion relation<jupyter_code>data = np.genfromtxt('./Downloads/table_2.dat')
x = data[:,0]
y = data[:,1]
v = data[:,2]
ve = data[:,3]<jupyter_output><empty_output><jupyter_text>(a) -- Minimization<jupyter_code>def L(params, x,y,v,ve):
    e0, vsys = params               
    v_model = vsys    #model without A                    
    sigma = np.sqrt(e0**2 + ve**2)                          
    L = ((v - v_model)**2/ (sigma**2) + np.log(2*pi*sigma**2))
    return 0.5 * np.sum(L)        
log_BFGS = minimize(L,x0=[1,1],args=(x,y,v,ve))
print "Was minimization successful?",log_BFGS.success
print '%6s %6s' % ('sig0(rad)','v_sys(km/s)')
print '%.3f %.3f' % (log_BFGS.x[0],log_BFGS.x[1])<jupyter_output>Was minimization successful? True
sig0(rad) v_sys(km/s)
127.074 134.081<jupyter_text>(b) -- Minimization (with $A$ and $\theta_0$)<jupyter_code>def logL(params,x,y,v,ve):
    e0, vsys, A, theta0 = params 
    theta = np.arctan2(y,x)
    v_model = vsys + A*np.sin(theta-theta0)    #model with A                    
    sigma = np.sqrt(e0**2 + ve**2)                          
    L = ((v - v_model)**2/ (sigma**2) + np.log(2*pi*sigma**2))
    return 0.5 * np.sum(L)  

ln_BFGS = minimize(logL,x0=[1,1,1,1],args=(x,y,v,ve))
print "Was minimization successful?",ln_BFGS.success
print '%6s %6s %4s %8s' % ('sig0(rad)','v_sys(km/s)','A(km/s)','theta0(rad)')
print '%.6f %.6f %.6f %.6f' % (ln_BFGS.x[0],ln_BFGS.x[1],ln_BFGS.x[2],ln_BFGS.x[3])<jupyter_output>Was minimization successful? True
sig0(rad) v_sys(km/s) A(km/s) theta0(rad)
120.918200 132.610064 56.638726 0.812978<jupyter_text>(c) -- MCMC<jupyter_code>import emcee

def LMC(params,x,y,v,ve):                         
    e0, vsys, A, theta0 = params 
    theta = np.arctan2(y,x)
    v_model = vsys + A*np.sin(theta-theta0)                                               
    sigma = np.sqrt(e0**2 + ve**2)                          
    L = ((v - v_model)**2/ (sigma**2) + np.log(2*pi*sigma**2))
    if A > 0: #A stays positive
        return -0.5 * np.sum(L)
    else:
        return -np.inf
nwalkers, ndim = 500, 4
p0 = np.zeros(shape=(500,4))
for i in range(500):
    p0[i,0] = ln_BFGS.x[0]+1e-4*np.random.uniform(-ln_BFGS.x[0]/10, ln_BFGS.x[0]/10)
    p0[i,1] = ln_BFGS.x[1]+1e-4*np.random.uniform(-ln_BFGS.x[1]/10, ln_BFGS.x[1]/10)
    p0[i,2] = ln_BFGS.x[2]+1e-4*np.random.uniform(-ln_BFGS.x[2]/10, ln_BFGS.x[2]/10)
    p0[i,3] = ln_BFGS.x[3]+1e-6*np.random.uniform(-ln_BFGS.x[3]/10, ln_BFGS.x[3]/10)    
    
l = emcee.EnsembleSampler(nwalkers, ndim, LMC, args=[x,y,v,ve])

MCMC = l.run_mcmc(p0, 500)
e0 = [i[0] for i in MCMC[0]]
vsys = [i[1] for i in MCMC[0]]     
A = [i[2] for i in MCMC[0]]
theta0 = [i[3] for i in MCMC[0]]

#print the result
print '%6s %6s %4s %8s' % ('sig0(rad)','v_sys(km/s)','A(km/s)','theta0(rad)')
print '%.6f %.6f %.6f %.6f' % (np.mean(e0),np.mean(vsys),np.mean(A),np.mean(theta0))<jupyter_output>sig0(rad) v_sys(km/s) A(km/s) theta0(rad)
124.025509 131.732557 54.471844 0.806364<jupyter_text>(d) -- Histogram<jupyter_code>plt.hist(e0, bins=100)
plt.xlabel(r"$\sigma_0(rad)$")
plt.ylabel('occurance')
plt.axvline(x=ln_BFGS.x[0],color='#9370db',ls='--')
plt.axvline(x=np.mean(e0),color='#b22222',ls='--')
plt.legend(['BFGS','MCMC'])
plt.show()
plt.hist(vsys, bins=100)
plt.xlabel(r"$v_{sys}(km/s)$")
plt.ylabel('occurance')
plt.axvline(x=ln_BFGS.x[1],color='#9370db',ls='--')
plt.axvline(x=np.mean(vsys),color='#b22222',ls='--')
plt.legend(['BFGS','MCMC'])
plt.show()
plt.hist(A, bins=100)
plt.xlabel(r"$A(km/s)$")
plt.ylabel('occurance')
plt.axvline(x=ln_BFGS.x[2],color='#9370db',ls='--')
plt.axvline(x=np.mean(A),color='#b22222',ls='--')
plt.legend(['BFGS','MCMC'])
plt.show()
plt.hist(theta0, bins=100)
plt.xlabel(r"$\theta_0(rad)$")
plt.ylabel('occurance')
plt.axvline(x=ln_BFGS.x[3],color='#9370db',ls='--')
plt.axvline(x=np.mean(theta0),color='#b22222',ls='--')
plt.legend(['BFGS','MCMC'])
plt.show()<jupyter_output><empty_output><jupyter_text>(e) -- $v_{rot} \ vs \theta$<jupyter_code>from matplotlib.pyplot import figure,show

params1 = [ln_BFGS.x[1],ln_BFGS.x[2],ln_BFGS.x[3]]
params2 = [np.mean(vsys),np.mean(A),np.mean(theta0)]

def v_model(params,x,y):
    vsys, A, theta0 = params
    theta = np.arctan2(y,x)
    v_rot = vsys + A*np.sin(theta-theta0)
    return theta, v_rot
Theta1, V1 = v_model(params1,x,y)
Theta2, V2 = v_model(params2,x,y)

#plot data and fits
fig = figure(figsize=(12,9))
frame = fig.add_subplot(1,1,1)
frame.errorbar(np.arctan2(y,x),v,ve,marker = 'd',ls='',color='indigo',fillstyle='none')
frame.plot(Theta1,V1,'.',color = 'darkorange')
frame.plot(Theta2,V2,'.',color='#b22222')
frame.set_xlabel(r'$\theta(radian)$')
frame.set_ylabel(r'$v_{rot}(km/s)$')
frame.legend(['data','BFGS','MCMC'])
show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.21202979164702554, 'fraction_numerical': 0.0663712642594513, 'mean_word_length': 4.615669666490207, 'pattern_counts': {'":': 0, '<': 40, '<?xml version=': 0, '>': 41, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '511732', 'license': 'MIT', 'n_tokens_mistral': 5054, 'n_tokens_neox': 4715, 'n_words': 845, 'path': 'numerical6.ipynb', 'repo_name': 'fatginger1024/NumericalMethods', 'size': '10607'}