Monte Carlo Simulations
In [1]:
## importing the libraries
import numpy as np
import math as m
import scipy
import matplotlib.pyplot as plt
from scipy import special
from scipy.optimize import curve_fit
In [2]:
def Erfc(x,sigma):
y = special.erfc(x/(sigma*np.sqrt(2)))
return y
## Introducing the model to be used later for the fitting
def SL_fit(x,w1,w2,w3,s1,s2,s3):
f = 0.5*(w1*Erfc(x,s1)+w2*Erfc(x,s2)+w3*Erfc(x,s3)+ (1-w1-w2-w3))
return f
In [3]:
## Jacobian matrix
def DF(x,w1,w2,w3,s1,s2,s3):
dfw1 = 0.5*(Erfc(x,s1)-1)
dfw2 = 0.5*(Erfc(x,s2)-1)
dfw3 = 0.5*(Erfc(x,s3)-1)
dfs1 = w1*x*np.exp(-x**2/(2*s1**2))/ (np.sqrt(2*np.pi)*s1**2)
dfs2 = w2*x*np.exp(-x**2/(2*s2**2))/ (np.sqrt(2*np.pi)*s2**2)
dfs3 = w3*x*np.exp(-x**2/(2*s3**2))/ (np.sqrt(2*np.pi)*s3**2)
return np.transpose(np.array([dfw1,dfw2,dfw3,dfs1,dfs2,dfs3]))
In [4]:
## Loading the data:
path = "/home/fatima/Desktop/project_3/"
file = np.loadtxt(path+'limb_profile_av_norm_shifted')
x = file[:,0]
y = file[:,1]
In [5]:
ind = np.where(x>=0)
x = x[ind]
y = y[ind]
weights = np.sqrt(np.abs(y)) ## Poisson weighting
In [6]:
p0=[0.3, 0.3, 0.2, 1, 2, 3] ## initial guess best-fit parameters
popt, pcov = curve_fit(SL_fit,x,y,p0,method='lm',sigma=weights,jac=DF,ftol=1e-8,xtol=1e-8,maxfev=5000)
chi_sq_w = np.sum((1/weights**2)*(SL_fit(x,*popt)-y)**2)
red_chi_sq = chi_sq_w/(len(y)-len(popt))
print popt # to print the best-fit parameters
In [94]:
perr = np.sqrt(np.diag(pcov))
print perr ## to print the error
In [9]:
mc_pars = []
total_iterations = 10
for _ in range(total_iterations): # i was replaced by _ here, because it was not used
x_trial = x
y_trial = y + np.random.normal(loc=y, scale=weights, size=y.shape)
try:
iteration_identifiers, _ = curve_fit(SL_fit,x_trial,y_trial,p0,method='lm')
except (ValueError, RuntimeError):
continue
mc_pars.append(iteration_identifiers)
print iteration_identifiers
plt.plot(x_trial, y_trial)
plt.plot(x,y)
plt.show()
In [87]:
mc_pars = np.array(mc_pars)
print(mc_pars.shape)
In [90]:
print(mc_pars[:, 2].mean())
In [11]:
X = np.random.normal(loc=y, scale=weights, size=y.shape)
In [12]:
X.shape
Out[12]:
In [13]:
plt.plot(X)
Out[13]:
In [17]:
plt.plot(y+X)
Out[17]:
In [3]:
import random
def f( x, p0, p1, p2):
return p0*x + 0.4*np.sin(p1*x) + p2
def ff(x, p):
return f(x, *p)
# These are the true parameters
p0 = 1.0
p1 = 40
p2 = 2.0
# These are initial guesses for fits:
pstart = [
p0 + random.random(),
p1 + 5.*random.random(),
p2 + random.random()
]
%matplotlib inline
import matplotlib.pyplot as plt
xvals = np.linspace(0., 1, 120)
yvals = f(xvals, p0, p1, p2)
In [4]:
pstart
Out[4]:
In [40]:
xdata = np.array(xvals)
np.random.seed(42)
err_stdev = 10
yvals_err = np.random.normal(f(xdata, p0, p1, p2), err_stdev, len(xdata))
ydata = f(xdata, p0, p1, p2) + yvals_err
GAUSSIAN = (1/ (np.sqrt(2.0*np.pi)*err_stdev))* np.exp(-(ydata)*(ydata)/(2.0*err_stdev**2))
In [41]:
plt.hist(yvals_err)
#plt.plot(yvals_err, GAUSSIAN)
#plt.plot()
Out[41]:
In [19]:
plt.hist(f(xdata, p0, p1, p2))
Out[19]:
In [20]:
plt.plot(xdata,ydata)
Out[20]:
In [21]:
f(xdata, p0, p1, p2)
Out[21]:
In [22]:
yvals_err
Out[22]:
In [23]:
plt.hist(ydata)
Out[23]:
In [ ]: