In [1]:

    

%pylab


    

    




Welcome to pylab, a matplotlib-based Python environment [backend: TkAgg].
For more information, type 'help(pylab)'.

In [3]:

    

def function(x):
    y = exp(-((x-30.0)**2)/100.0)
    return y


    

In [63]:

    

X = linspace(0, 60, 100)
Y = function(X)


    

In [62]:

    

plot(x, y)


    

    
Out[62]:





[<matplotlib.lines.Line2D at 0x5447410>]






    

In [11]:

    

#1. Make a random number between 0 and 60


    

In [55]:

    

walk = []
x = random.random()*60
walk.append(x)


    

In [56]:

    

for i in range(200000):
    x = random.random()*2-1
    alpha = function(x + walk[-1])/function(walk[-1])

    if alpha>=1.0:
        walk.append(walk[-1]+x)
    else:
        beta = random.random()
        if(beta<=alpha):
            walk.append(walk[-1]+x)
        else:
            walk.append(walk[-1])


    

In [69]:

    

histo = hist(walk, bins=40, normed=True)
f = function(X)
norm = sum(f*(x[1]-x[0]))
plot(X,f/norm, linewidth=1, color='r')


    

    
Out[69]:





[<matplotlib.lines.Line2D at 0x6ca45d0>]






    

In [50]:

    

size(walk)


    

    
Out[50]:





200001

In [2]:

    

def function(x, y):
    f = exp((-(100*(y-x**2)**2)+(1-x)**2)/20.)
    return f


    

In [3]:

    

from mpl_toolkits.mplot3d import Axes3D


    

In [33]:

    

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
X = np.linspace(-20, 20, 300)
Y = np.linspace(-20, 200, 300)
X, Y = np.meshgrid(X,Y)
Z = function(X, Y)
ax.plot_surface(X, Y, log10(Z))
show()


    

    




-c:7: RuntimeWarning: divide by zero encountered in log10
/usr/lib/pymodules/python2.7/mpl_toolkits/mplot3d/axes3d.py:1618: RuntimeWarning: invalid value encountered in subtract
  v1[which_pt] = np.array(ps2[i1]) - np.array(ps2[i2])
/usr/lib/pymodules/python2.7/mpl_toolkits/mplot3d/axes3d.py:1619: RuntimeWarning: invalid value encountered in subtract
  v2[which_pt] = np.array(ps2[i2]) - np.array(ps2[i3])
/usr/lib/python2.7/dist-packages/numpy/core/numeric.py:1342: RuntimeWarning: invalid value encountered in subtract
  y = a[2]*b[0] - a[0]*b[2]
/usr/lib/pymodules/python2.7/mpl_toolkits/mplot3d/axes3d.py:1673: RuntimeWarning: invalid value encountered in divide
  for n in normals])

In [115]:

    

print amin(Y)


    

    




-100.0

In [31]:

    

hist = hist2d(X, Y, weights=log10(Z), bins=(60, 60), normed=True)
colorbar()


    

    




/usr/lib/pymodules/python2.7/matplotlib/colorbar.py:808: RuntimeWarning: invalid value encountered in divide
  z = np.take(y, i0) + (xn-np.take(b,i0))*dy/db






    
Out[31]:





<matplotlib.colorbar.Colorbar instance at 0x41c4758>

In [6]:

    

%pylab inline
#load the list of points
obs_data = loadtxt("movimiento.dat")
t_obs = obs_data[:,0]
x_obs = obs_data[:,1]


    

    




Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline].
For more information, type 'help(pylab)'.

In [7]:

    

scatter(t_obs,x_obs)


    

    
Out[7]:





<matplotlib.collections.PathCollection at 0x40afd10>






    

In [4]:

    

def likelihood(y_obs, y_model):
    chi_squared = (1.0/2.0)*sum((y_obs-y_model)**2)
    return exp(-chi_squared)


    

In [5]:

    

def my_model(t_obs, x_obs, a, b):
    return a*t_obs + b*t_obs**2


    

In [ ]:

    

a_walk = []
b_walk = []
l_walk = []

m_walk = append(random.random())
b_walk = append(random.random())

y_init = my_model(t_obs, m_walk[0], b_walk[0])
l_walk = append(l_walk, likelihood(y_obs, y_init))
print m_walk
print b_walk
print l_walk