%matplotlib inline from IPython.display import Image import numpy as np import matplotlib.pyplot as plt from scipy.integrate import quad, dblquad, tplquad # defino una función a integrar def f(x): return x x_lower = 0 # límite inferior de x x_upper = 1 # límite superior de x val, abserr = quad(f, x_lower, x_upper) print "integral value =", val, ", absolute error =", abserr def integrand(x, n): return x*n x_lower = 0 # the lower limit of x x_upper = 1 # the upper limit of x val, abserr = quad(integrand, x_lower, x_upper, args=3) print val, abserr val, abserr = quad(lambda x: np.exp(-x ** 2), -np.inf, np.inf) # Inf = infinito print "numerical =", val, abserr analytical = np.sqrt(np.pi) print "analytical =", analytical def integrand(x, y): return np.exp(-x**2-y**2) x_lower = 0 x_upper = 10 y_lower = 0 y_upper = 10 val, abserr = dblquad(integrand, x_lower, x_upper, lambda x : y_lower, lambda x: y_upper) print val, abserr import scipy.linalg as la A = np.array([[1,2,3], [5,6,7], [9,10,11]]) b = np.array([1,2,3]) x = la.solve(A, b) x np.dot(A, x) - b # chequeo A = np.random.rand(3,3) B = np.random.rand(3,3) X = la.solve(A, B) X la.norm(np.dot(A, X) - B) # chequeo evals = la.eigvals(A) evals evals, evecs = la.eig(A) evals evecs n = 1 la.norm(np.dot(A, evecs[:,n]) - evals[n] * evecs[:,n]) la.inv(A) # inversa de una matrix la.det(A) # determinante la.norm(A, ord=2), la.norm(A, ord=np.inf) # normas de varios ordenes from scipy import optimize def f(x): return 4*x**3 + (x-2)**2 + x**4 fig, ax = plt.subplots() x = np.linspace(-5, 3, 100) ax.plot(x, f(x)); x_min = optimize.fmin_bfgs(f, -2) x_min optimize.fmin_bfgs(f, 0.5) optimize.brent(f) optimize.fminbound(f, -4, 2) def f(x): return np.sin(x) fig, ax = plt.subplots(figsize=(10,4)) x = np.linspace(-10, 10, 50) y = f(x) ax.plot(x, y) optimize.fsolve(f, 0.5) optimize.fsolve(f, 4.0) optimize.fsolve(f, -3.0) from scipy import interpolate def f(x): return np.sin(x) n = np.arange(0, 10) x = np.linspace(0, 9, 100) y_meas = f(n) + 0.1 * np.random.randn(len(n)) y_real = f(x) linear_interpolation = interpolate.interp1d(n, y_meas) y_interp1 = linear_interpolation(x) cubic_interpolation = interpolate.interp1d(n, y_meas, kind='cubic') y_interp2 = cubic_interpolation(x) fig, ax = plt.subplots(figsize=(10,4)) ax.plot(n, y_meas, 'bs', label='datos originales') ax.plot(x, y_real, 'k', lw=2, label='funcion real') ax.plot(x, y_interp1, 'r', label='interpolacion lineal') ax.plot(x, y_interp2, 'g', label='interpolacion cubica') ax.legend(loc=3); from scipy import stats X = stats.poisson(3.5) # crea una variable aleatoria discreta con distribución de Poisson n = np.arange(0,15) fig, axes = plt.subplots(3,1, sharex=True) axes[0].step(n, X.pmf(n)) # (PMF) axes[1].step(n, X.cdf(n)) # (CDF) axes[2].hist(X.rvs(size=1000)); #histograma de 1000 realizaciones aleatorias de la variable estocástica X Y = stats.norm() # crea una variable aleatoria continua con distribución normal x = np.linspace(-5,5,100) fig, axes = plt.subplots(3,1, sharex=True) axes[0].plot(x, Y.pdf(x)) axes[1].plot(x, Y.cdf(x)); axes[2].hist(Y.rvs(size=1000), bins=50); X.mean(), X.std(), X.var() # distribución de Poisson Y.mean(), Y.std(), Y.var() # distribución normal t_statistic, p_value = stats.ttest_ind(X.rvs(size=1000), X.rvs(size=1000)) print "t-statistic =", t_statistic print "p-value =", p_value stats.ttest_1samp(Y.rvs(size=1000), 0.1) Y.mean() stats.ttest_1samp(Y.rvs(size=1000), Y.mean()) from scipy import ndimage from scipy import misc lena = misc.lena() shifted_lena = ndimage.shift(lena, (50, 50)) shifted_lena2 = ndimage.shift(lena, (50, 50), mode='nearest') rotated_lena = ndimage.rotate(lena, 30) cropped_lena = lena[50:-50, 50:-50] zoomed_lena = ndimage.zoom(lena, 2) fig, axes = plt.subplots(1,5, figsize=(16, 5)) axes[0].imshow(shifted_lena, cmap=plt.cm.gray) axes[1].imshow(shifted_lena2, cmap=plt.cm.gray) axes[2].imshow(rotated_lena, cmap=plt.cm.gray) axes[3].imshow(cropped_lena, cmap=plt.cm.gray) axes[4].imshow(zoomed_lena, cmap=plt.cm.gray) for ax in range(5): axes[ax].axis('off') import numpy as np noisy_lena = np.copy(lena).astype(np.float) noisy_lena += lena.std()*0.5*np.random.standard_normal(lena.shape) blurred_lena = ndimage.gaussian_filter(noisy_lena, sigma=3) median_lena = ndimage.median_filter(blurred_lena, size=5) from scipy import signal wiener_lena = signal.wiener(blurred_lena, (5,5)) fig, axes = plt.subplots(1,4, figsize=(16, 5)) axes[0].imshow(noisy_lena, cmap=plt.cm.gray) axes[1].imshow(blurred_lena, cmap=plt.cm.gray) axes[2].imshow(median_lena, cmap=plt.cm.gray) axes[3].imshow(wiener_lena, cmap=plt.cm.gray) for ax in range(4): axes[ax].axis('off') im = np.zeros((256, 256)) im[64:-64, 64:-64] = 1 im = ndimage.rotate(im, 15, mode='constant') im = ndimage.gaussian_filter(im, 8) sx = ndimage.sobel(im, axis=0, mode='constant') sy = ndimage.sobel(im, axis=1, mode='constant') sob = np.hypot(sx, sy) plt.figure(figsize=(16, 5)) plt.subplot(141) plt.imshow(im, cmap=plt.cm.gray) plt.axis('off') plt.title('Cuadrado', fontsize=16) plt.subplot(142) plt.imshow(sx) plt.axis('off') plt.title('Sobel (direccion x)', fontsize=16) plt.subplot(143) plt.imshow(sob) plt.axis('off') plt.title('Sobel', fontsize=16) im += 0.07*np.random.random(im.shape) sx = ndimage.sobel(im, axis=0, mode='constant') sy = ndimage.sobel(im, axis=1, mode='constant') sob = np.hypot(sx, sy) plt.subplot(144) plt.imshow(sob) plt.axis('off') plt.title('Sobel con ruido', fontsize=20) plt.subplots_adjust(wspace=0.02, hspace=0.02, top=1, bottom=0, left=0, right=0.9) plt.show() np.random.seed(1) n = 10 l = 256 im = np.zeros((l, l)) points = l*np.random.random((2, n**2)) im[(points[0]).astype(np.int), (points[1]).astype(np.int)] = 1 im = ndimage.gaussian_filter(im, sigma=l/(4.*n)) mask = (im > im.mean()).astype(np.float) img = mask + 0.3*np.random.randn(*mask.shape) binary_img = img > 0.5 # Eliminamos los pequeños puntos blancos open_img = ndimage.binary_opening(binary_img) # Eliminamos los pequeños puntos negros close_img = ndimage.binary_closing(open_img) plt.figure(figsize=(12, 3)) l = 128 plt.subplot(141) plt.imshow(binary_img[:l, :l], cmap=plt.cm.gray) plt.axis('off') plt.subplot(142) plt.imshow(open_img[:l, :l], cmap=plt.cm.gray) plt.axis('off') plt.subplot(143) plt.imshow(close_img[:l, :l], cmap=plt.cm.gray) plt.axis('off') plt.subplot(144) plt.imshow(mask[:l, :l], cmap=plt.cm.gray) plt.contour(close_img[:l, :l], [0.5], linewidths=2, colors='r') plt.axis('off') plt.subplots_adjust(wspace=0.02, hspace=0.3, top=1, bottom=0.1, left=0, right=1) plt.show()