# 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/"
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

[ 0.52750103  0.28882569  0.10191755  0.25905336  0.76540583  2.83343005]

In [94]:
perr = np.sqrt(np.diag(pcov))
print perr ## to print the error

[ 0.00349415  0.0028457   0.00119396  0.00135734  0.00715294  0.04738134]

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()

[ 0.31340342  0.31363243  0.21164169  0.52898839  1.80269277  1.80271224]

[ 0.31885237  0.31594929  0.20553026  1.14083605  1.14115389  1.14081575]

[ 0.34207102  0.34599831  0.1371758   0.62452273  2.04366996  0.6243128 ]

[ 0.20582384  0.48394574  0.25346475  0.2518755   2.73310171  0.96326215]

[ 0.30655759  0.29281329  0.15200812  0.74734774  0.7473279   0.74721024]

[ 0.6408618   0.16169021  0.09391276  1.04628039  4.13653741  4.1316929 ]

In [87]:
mc_pars = np.array(mc_pars)
print(mc_pars.shape)

(666, 6)

In [90]:
print(mc_pars[:, 2].mean())

1.33754544667

In [11]:
X = np.random.normal(loc=y, scale=weights, size=y.shape)

In [12]:
X.shape

Out[12]:
(279,)
In [13]:
plt.plot(X)

Out[13]:
[<matplotlib.lines.Line2D at 0x7f2d6291acd0>]
In [17]:
plt.plot(y+X)

Out[17]:
[<matplotlib.lines.Line2D at 0x7f2d628f5ed0>]
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]:
[1.316042402611915, 43.27641975825924, 2.72916543068528]
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]:
(array([  1.,   6.,  12.,  12.,  32.,  29.,  13.,   7.,   7.,   1.]),
array([-23.67796735, -18.54075949, -13.40355162,  -8.26634376,
-3.1291359 ,   2.00807197,   7.14527983,  12.28248769,
17.41969556,  22.55690342,  27.69411128]),
<a list of 10 Patch objects>)
In [19]:
plt.hist(f(xdata, p0, p1, p2))

Out[19]:
(array([  8.,   9.,  13.,  16.,  20.,  14.,  15.,  10.,   9.,   6.]),
array([ 1.71765553,  1.88414245,  2.05062937,  2.2171163 ,  2.38360322,
2.55009015,  2.71657707,  2.88306399,  3.04955092,  3.21603784,
3.38252477]),
<a list of 10 Patch objects>)
In [20]:
plt.plot(xdata,ydata)

Out[20]:
[<matplotlib.lines.Line2D at 0x7f2249f26d10>]
In [21]:
f(xdata, p0, p1, p2)

Out[21]:
array([ 2.        ,  2.1403395 ,  2.26591185,  2.36360271,  2.42341847,
2.43960469,  2.41129028,  2.34258492,  2.24211913,  2.12207826,
1.99683862,  1.88135842,  1.78950351,  1.73249545,  1.71765553,
1.74758527,  1.81987531,  1.92737503,  2.05899292,  2.200938  ,
2.33826341,  2.45653935,  2.54346816,  2.59026076,  2.59262038,
2.55122346,  2.47164399,  2.36372958,  2.24049927,  2.11668637,
2.00708938,  1.92491568,  1.88030323,  1.87918593,  1.92262939,
2.00671168,  2.12296232,  2.25931035,  2.40143536,  2.53437035,
2.64417691,  2.71950535,  2.75286497,  2.74146251,  2.68751475,
2.59800044,  2.48387918,  2.35886472,  2.23789005,  2.13543601,
2.06391049,  2.03225964,  2.04496661,  2.10154969,  2.19661631,
2.32046653,  2.46017881,  2.60105617,  2.72827127,  2.82852591,
2.89153951,  2.91119972,  2.88624661,  2.82041365,  2.72200985,
2.60298977,  2.47761548,  2.36086025,  2.26673264,  2.2067086 ,
2.18844696,  2.21493224,  2.28414061,  2.38926636,  2.51948369,
2.66115839,  2.79937385,  2.91960066,  3.00932281,  3.05943856,
3.06527919,  3.02713153,  2.95020589,  2.84405285,  2.72149429,
2.59718835,  2.48598871,  2.40128212,  2.35349006,  2.34890229,
2.38897286,  2.47015737,  2.58430967,  2.71959366,  2.86180802,
2.99597575,  3.10802049,  3.186342  ,  3.22311459,  3.21516299,
3.16431777,  3.07721039,  2.96453106,  2.83983214,  2.71801131,
2.61364411,  2.53935256,  2.50439241,  2.5136172 ,  2.56693499,
2.65931865,  2.78136856,  2.92036466,  3.06169016,  3.19046753,
3.29322374,  3.35939819,  3.38252477,  3.36095555,  3.29804526])
In [22]:
yvals_err

Out[22]:
array([ 0.09934283, -0.02765286,  0.12953771,  0.30460597, -0.04683067,
-0.04682739,  0.31584256,  0.15348695, -0.09389488,  0.10851201,
-0.09268354, -0.09314595,  0.04839245, -0.38265605, -0.34498357,
-0.11245751, -0.20256622,  0.06284947, -0.18160482, -0.28246074,
0.29312975, -0.04515526,  0.01350564, -0.28494964, -0.10887654,
0.02218452, -0.23019872,  0.0751396 , -0.12012774, -0.05833875,
-0.12034132,  0.37045564, -0.00269944, -0.21154219,  0.16450898,
-0.24416873,  0.04177272, -0.39193402, -0.26563721,  0.03937225,
0.14769332,  0.03427366, -0.02312966, -0.06022074, -0.2957044 ,
-0.14396884, -0.09212775,  0.21142445,  0.06872366, -0.35260803,
0.06481679, -0.07701646, -0.1353844 ,  0.12233526,  0.2061999 ,
0.18625602, -0.1678435 , -0.06184248,  0.06625269,  0.19510903,
-0.09583485, -0.0371318 , -0.22126699, -0.23924132,  0.16250516,
0.27124801, -0.01440202,  0.20070658,  0.07232721, -0.12902395,
0.07227912,  0.30760731, -0.00716521,  0.31292873, -0.52394902,
0.1643805 ,  0.01740941, -0.05980147,  0.01835216, -0.39751378,
-0.04393438,  0.07142251,  0.29557881, -0.10365404, -0.16169872,
-0.10035141,  0.18308042,  0.06575022, -0.10595204,  0.10265349,
0.01941551,  0.193729  , -0.14041062, -0.06553243, -0.07842163,
-0.29270299,  0.05922406,  0.05221105,  0.00102269, -0.04691743,
-0.28307415, -0.08412906, -0.0685429 , -0.16045545, -0.03225714,
0.08081017,  0.37723718,  0.03491556,  0.05151008, -0.01488918,
-0.38375424, -0.00530278,  0.01204604,  0.49264842, -0.03847219,
0.06030947, -0.00694235, -0.23373561,  0.22856456,  0.15038661])
In [23]:
plt.hist(ydata)

Out[23]:
(array([  2.,   6.,  11.,  13.,  23.,  26.,  21.,   6.,   9.,   3.]),
array([ 1.3498394 ,  1.57380747,  1.79777555,  2.02174362,  2.24571169,
2.46967976,  2.69364783,  2.9176159 ,  3.14158397,  3.36555204,
3.58952011]),
<a list of 10 Patch objects>)
In [ ]: