<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 SciPy
SciPy es paquete que incluye una colección de algoritmos matemáticos y funciones construidas sobre el paquete NumPy. 
Dirigida al mismo tipo de usuarios que los de aplicaciones como MATLAB, GNU Octave, y Scilab, SciPy contiene varios módulos, entre ellos
Step1:  File input/output
Step2:  Así de sencillo, ya hemos "exportado" los valores de la variable 'a' a un archivo en formato *.mat.
Veamos si podemos localizarlo
Step3:  Imaginemos la situación inversa. ¿Cómo importar archivos de datos de MATLAB a Python
Step4:  Podemos listar los variables importadas
Step5:  Existen otras formas de cargar datos
Step6:  Vaya, parece que para estas muestras los histogramas no nos muestran información clara tampoco. Asumiendo que las muestras son representativas, podemos realizar la prueba de t de Student para ver si existen diferencias significativas entre ellas.
Step7:  Podemos ver que el valor de $\frac{p}{2} > 0.05 $ por lo que la hipótesis nula se cumple (no existen diferencias significativas entre ellas).
Ejercicio
Step8:  Ahora que hemos realizado nuestras simulaciones, podemos calcular nuestro valor $p$, que es simplemente la proporción de simulaciones que resultaron en una diferencia mayor o igual a 6.6 (la diferencia original)
Step9:  $$ p = \frac{N_{>6.6}}{N_{total}} = \frac{1512}{10000} = 0.15 $$
Step10:  Este tutorial de scipy.stats muestra más ejemplos que podemos realizar. Por el momento, vamos a continuar explorando SciPy. Retomaremos más trabajo estadísitco cuando lleguemos a  pandas, 
Interpolación
Step11:  Para interpolar utilizaremos el paquete interpolate de SciPy
Step12:  Para crear una función interpolante utilizaremos el objeto InterpolatedUnivariateSpline del paquete interpolate. A este objeto solo hay que pasarle los puntos de interpolación y el grado, y generará un spline.
Step13:  ¿Cómo obtengo los puntos desde aquí? El resultado que hemos obtenido es una función y admite como argumento la $x$.
Step14:  Vamos a representar esta función junto con los puntos de interpolación. Fíjate en que, ahora que tenemos una función interpolante, podemos representarla en un dominio
Step15:  Retrocede ahora y comprueba lo que pasa si cambias el grado del spline. Dale a un vistazo a todas las opciones que SciPy ofrece para interpolar datos.
Ajuste
El ajuste funciona de manera totalmente distinta a la interpolación
Step16:  Generaremos unos datos para ver cómo funcionaría, del tipo
Step17:  Utilicemos ahora la función polynomial.polyfit, que recibe los puntos de interpolación y el grado del polinomio. El resultado serán los coeficientes del mismo, en orden de potencias crecientes.
Step18:  ¡Muy similares a lo que esperábamos! Para evaluar un polinomio con estos coeficientes, o bien construimos la función nosotros mismos o usamos la función polynomial.polyval
Step19:  Si la función que queremos ajustar es más compleja, necesitaremos ajustar los datos a una curva mediante un algoritmo de optimización.
Step20:  Generemos una vez más los datos añadiendo un poco de ruído. ¿Puedes leer ya funciones descritas con NumPy?
Step21:  Como en este ejemplo sintético conocemos los valores exactos, podemos ver que existe variación respecto a los valores originales debido al ruído añadido.
Normalmente es mejor ayudar al solver definiendo unas restricciones
Step22:  Veamos los resultados en una gráfica
Step23:  Otra forma de ajustar una función a datos experimentales es minimizando el error por mínimos cuadrados. Para hacer este ejemplo más interesante, añadamos además de valores atípicos (outliers, en inglés). Este ejemplo está tomado del Cookbook de Scipy, robust regression
Step24:  Una vez creado el modelo estamos listos para generar los datos.
Step25:  La función que calcula los residuos se puede definir como
Step26:  Ya tenemos todo lo que necesitamos para realizar el ajuste por mínimos cuadrados
Step27:  Veamos los resultados
Step28:  El paquete optimize incluye multitud de métodos para optimización, ajuste de curvas y búsqueda de raíces. La ayuda de este paquete es bastante larga (puedes consultarla también en http
Step29:  Lo más sencillo en estos casos es aplicar un filtro por ventana
Step30:  Como podemos ver, si la frecuencia de muestreo no es muy alta o el tamaño de nuestra ventana no es el adecuado, el resultado puede ser no satisfactorio.
El filtro de Savitzky–Golay suele ser interesante en estos casos.
Step31:  Si estamos interesados en obtener una señal sinuisidal y dado que el ruído ocurre a una frecuencia más alta, otra opción es generar filtro de paso bajo.
Step32:  Por último, si la señal tiene un deriva (drifting) podemos corregirla fácilmente con
<ASSISTANT_TASK:>
Python Code:
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
Explanation: SciPy
SciPy es paquete que incluye una colección de algoritmos matemáticos y funciones construidas sobre el paquete NumPy. 
Dirigida al mismo tipo de usuarios que los de aplicaciones como MATLAB, GNU Octave, y Scilab, SciPy contiene varios módulos, entre ellos:
Funciones especiales (scipy.special)
Integración (scipy.integrate)
Optimización (scipy.optimize)
Interpolación (scipy.interpolate)
Transformadas de Fourier (scipy.fftpack)
Procesado de señales (scipy.signal)
Algebra lineal (scipy.linalg)
Problemas de valores propios con matrices dispersas (ARPACK)
Grafos (scipy.sparse.csgraph)
Análisis espacial: (scipy.spatial)
Estadística (scipy.stats)
Tratamiento de imágenes multidimensionales (scipy.ndimage)
Entrada y salida de archivos (scipy.io)
En esta clase nos vamos repasar brevemente algunas de las funcionalidades más relevantes para el análisis de datos. Pero antes, carguemos las librerías de NumPy y matplotlib que ya conocemos.
End of explanation
from scipy import io as spio
a = np.ones((3, 3)) # creamos una matriz de 3x3 
spio.savemat('archivo.mat', # nombre del archivo
             {'a': a}) # asignamos y referenciamos el nombre con un diccionario 
Explanation: File input/output: scipy.io
Este módulo de SciPy permite leer y escribir archivos de datos en una gran variedad de formatos (ver documentación).
Por ejemplo, para trabajar con datos de MATLAB (*.dat):
End of explanation
%ls *.mat
Explanation: Así de sencillo, ya hemos "exportado" los valores de la variable 'a' a un archivo en formato *.mat.
Veamos si podemos localizarlo:
End of explanation
data = spio.loadmat('archivo.mat')
data['a']
Explanation: Imaginemos la situación inversa. ¿Cómo importar archivos de datos de MATLAB a Python
End of explanation
spio.whosmat('archivo.mat')
Explanation: Podemos listar los variables importadas:
End of explanation
a = np.array([84, 72, 57, 46, 63, 76, 99, 91])
b = np.array([81, 69, 74, 61, 56, 87, 69, 65, 66, 44, 62, 69])
plt.hist(b, bins=5, alpha=0.5)  
plt.hist(a, bins=5, alpha=0.5)  
plt.plot(a,np.zeros(len(a)),'^')
plt.plot(b,np.zeros(len(b)),'^')
plt.title("Histogram")
plt.show()
print("La media de 'a' es {0:.1f}, y su desviación estándar es {1:.1f}".format(a.mean(), a.std()))
print("La media de 'b' es {0:.1f}, y su desviación estándar es {0:.1f}".format(b.mean(), b.std()))
print("La diferencia entre las medias es de {0:.1f}".format(a.mean()- b.mean()))
Explanation: Existen otras formas de cargar datos:
* Datos de archivos de texto: numpy.loadtxt() | numpy.savetxt()
* Datos de archivos de texto incompletos: numpy.genfromtxt() | numpy.recfromcsv()
* Datos con formato NumPy-binario (eficiente): numpy.save() | numpy.load
* pandas.read_csv: que cubriremos en detalle más tarde
Estadística: scipy.stats
Este módulo contiene un gran número de distribuciones de probabilidad, tanto continuas como discretas, así como un creciente número de funciones estadísticas. 
Veamos un ejemplo simple tomado de Jake VanderPlas y su charla (Statistics for Hackers). Queremos saber si dos sets de muestras son diferentes (a y b).
End of explanation
from scipy import stats
stats.ttest_ind(a,b)
Explanation: Vaya, parece que para estas muestras los histogramas no nos muestran información clara tampoco. Asumiendo que las muestras son representativas, podemos realizar la prueba de t de Student para ver si existen diferencias significativas entre ellas.
End of explanation
samples = np.concatenate([a,b])
num_simulations = 10000
differences = np.zeros(num_simulations)
for i in range(num_simulations):
    np.random.shuffle(samples)
    a_new = samples[0:len(a)]
    b_new = samples[len(a):len(samples)]
    a_mean = a_new.mean()
    b_mean = b_new.mean()
    differences[i]= (a_mean-b_mean)
Explanation: Podemos ver que el valor de $\frac{p}{2} > 0.05 $ por lo que la hipótesis nula se cumple (no existen diferencias significativas entre ellas).
Ejercicio:
Si las muestras no son diferentes, no debería importar si intercambiamos los valores de forma aleatoria, ¿verdad?.
¡Comprobémoslo!
End of explanation
p = np.sum(differences>(a.mean()-b.mean()))/num_simulations
p
Explanation: Ahora que hemos realizado nuestras simulaciones, podemos calcular nuestro valor $p$, que es simplemente la proporción de simulaciones que resultaron en una diferencia mayor o igual a 6.6 (la diferencia original)
End of explanation
plt.hist(differences, bins=50)
plt.axvline((a.mean()-b.mean()),color='r')
plt.xlabel('mean difference')
plt.ylabel('number')
Explanation: $$ p = \frac{N_{>6.6}}{N_{total}} = \frac{1512}{10000} = 0.15 $$
End of explanation
x_i = [0.0, 0.9, 1.8, 2.7, 3.6, 4.4, 5.3, 6.2, 7.1, 8.0]
y_i = [0.0, 0.8, 1.0, 0.5, -0.4, -1.0, -0.8, -0.1, 0.7, 1.0]
plt.plot(x_i, y_i, 'x', mew=2)
Explanation: Este tutorial de scipy.stats muestra más ejemplos que podemos realizar. Por el momento, vamos a continuar explorando SciPy. Retomaremos más trabajo estadísitco cuando lleguemos a  pandas, 
Interpolación: scipy.interpolate
Una de las tareas fundamentales a las que nos enfrentaremos a diario en el análisis de datos es la de interpolar.
Supongamos que tenemos una serie de puntos que representan los datos de un cierto experimento. Por simplicidad, vamos a generar unos puntos de una función $\sin{x}$ de ejemplo para explicar el funcionamiento.
End of explanation
from scipy import interpolate
Explanation: Para interpolar utilizaremos el paquete interpolate de SciPy:
End of explanation
f_interp = interpolate.InterpolatedUnivariateSpline(x_i, y_i, k=1)
f_interp
Explanation: Para crear una función interpolante utilizaremos el objeto InterpolatedUnivariateSpline del paquete interpolate. A este objeto solo hay que pasarle los puntos de interpolación y el grado, y generará un spline.
End of explanation
f_interp(np.pi / 2)
Explanation: ¿Cómo obtengo los puntos desde aquí? El resultado que hemos obtenido es una función y admite como argumento la $x$.
End of explanation
x = np.linspace(0, 8)
y_interp = f_interp(x)
plt.plot(x_i, y_i, 'x', mew=2)
plt.plot(x, y_interp)
Explanation: Vamos a representar esta función junto con los puntos de interpolación. Fíjate en que, ahora que tenemos una función interpolante, podemos representarla en un dominio:
End of explanation
from numpy.polynomial import polynomial
Explanation: Retrocede ahora y comprueba lo que pasa si cambias el grado del spline. Dale a un vistazo a todas las opciones que SciPy ofrece para interpolar datos.
Ajuste
El ajuste funciona de manera totalmente distinta a la interpolación: obtendremos una curva que no pasará por ninguno de los puntos originales, pero que a cambio tendrá una expresión analítica simple que conocemos a priori.
Veamos un ejemplo simple para realizar ajustes polinómicos vamos a utilizar el paquete np.polynomial.polynomial (sí, está dos veces).
End of explanation
x_i = np.linspace(-2, 3, num=10)
y_i = x_i ** 2 - x_i + 1 + 0.5 * np.random.randn(10)
plt.plot(x_i, y_i, 'x', mew=2)
Explanation: Generaremos unos datos para ver cómo funcionaría, del tipo:
$$y(x) = x^2 - x + 1 + \text{ruido}$$
End of explanation
a, b, c = polynomial.polyfit(x_i, y_i, deg=2)
a, b, c
Explanation: Utilicemos ahora la función polynomial.polyfit, que recibe los puntos de interpolación y el grado del polinomio. El resultado serán los coeficientes del mismo, en orden de potencias crecientes.
End of explanation
x = np.linspace(-2, 3)
#y_fit = a + b * x + c * x ** 2
y_fit = polynomial.polyval(x, (a, b, c))
l, = plt.plot(x, y_fit)
plt.plot(x_i, y_i, 'x', mew=2, c=l.get_color())
Explanation: ¡Muy similares a lo que esperábamos! Para evaluar un polinomio con estos coeficientes, o bien construimos la función nosotros mismos o usamos la función polynomial.polyval:
End of explanation
from scipy.optimize import curve_fit
Explanation: Si la función que queremos ajustar es más compleja, necesitaremos ajustar los datos a una curva mediante un algoritmo de optimización.
End of explanation
def func(x, a, b, c):
    return a * np.exp(-b * x) + c
a, b, c = 2.5, 1.3, 0.5
xdata = np.linspace(0, 4, 50)
y = func(xdata, a, b, c)
y_noise = 1.5 * np.random.normal(size=xdata.size)
ydata = y + y_noise
plt.plot(xdata, ydata, 'x',mew=2, label='exp. data')
plt.plot(xdata, func(xdata, a, b, c), '-', label='true function')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.show()
popt, pcov = curve_fit(func, xdata, ydata)
popt
Explanation: Generemos una vez más los datos añadiendo un poco de ruído. ¿Puedes leer ya funciones descritas con NumPy?
End of explanation
popt_bounds, pcov_bounds = curve_fit(func, xdata, ydata,
                                     bounds=([1, 1, 0], [3., 2., 1.]))
popt_bounds
Explanation: Como en este ejemplo sintético conocemos los valores exactos, podemos ver que existe variación respecto a los valores originales debido al ruído añadido.
Normalmente es mejor ayudar al solver definiendo unas restricciones:
End of explanation
plt.plot(xdata, ydata, 'x',mew=2, label='exp. data')
plt.plot(xdata, func(xdata, a, b, c), '-', label='true function')
plt.plot(xdata, func(xdata, *popt), 'r-', label='fit')
plt.plot(xdata, func(xdata, *popt_bounds), 'g--', label='fit-with-bounds')
         
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.show()
Explanation: Veamos los resultados en una gráfica:
End of explanation
def generate_data(t, A, sigma, omega, noise=0, n_outliers=0, random_state=0):
    y = A * np.exp(-sigma*t) * np.sin(omega*t)
    rnd = np.random.RandomState(random_state)
    error = noise * rnd.randn(t.size)
    outliers = rnd.randint(0, t.size, n_outliers)
    error[outliers] = error[outliers] * 35
    return y + error
# Parametros del modelo
A = 2
sigma = 0.1
omega = 0.1 * 2 * np.pi
x_true = np.array([A, sigma, omega])
noise = 0.1
t_min = 0
t_max = 30
Explanation: Otra forma de ajustar una función a datos experimentales es minimizando el error por mínimos cuadrados. Para hacer este ejemplo más interesante, añadamos además de valores atípicos (outliers, en inglés). Este ejemplo está tomado del Cookbook de Scipy, robust regression
End of explanation
t= np.linspace(t_min, t_max, 30)
y_exp = generate_data(t, A, sigma, omega, noise=noise, n_outliers=4)
y_true = generate_data(t, A, sigma, omega) # ¿por qué no necesito indicar nada más?
plt.plot(t, y_exp, 'o', label='exp data')
plt.plot(t, y_true, label='true')
plt.xlabel('$t$')
plt.ylabel('$y$')
plt.legend()
Explanation: Una vez creado el modelo estamos listos para generar los datos.
End of explanation
def fun_res(x, t, y):
    A, sigma, omega = x # parámetros 
    return (A * np.exp(-sigma * t) * np.sin(omega * t)) - y
x0 = np.ones(3) # valores inciales de A, sigma y omega
Explanation: La función que calcula los residuos se puede definir como:
End of explanation
from scipy.optimize import least_squares
res_lsq = least_squares(fun_res, x0, args=(t, y_exp))
res_lsq
res_robust = least_squares(fun_res, x0,
                           loss='soft_l1', # Norma del tipo L1 (más robusta)
                           f_scale=0.1, # restringe los errores
                           args=(t, y_exp))
res_robust
Explanation: Ya tenemos todo lo que necesitamos para realizar el ajuste por mínimos cuadrados
End of explanation
y_lsq = generate_data(t, *res_lsq.x)
y_robust = generate_data(t, *res_robust.x)
plt.plot(t, y_exp, 'o', label='exp data')
plt.plot(t, y_true, label='true')
plt.plot(t, y_lsq, label='lsq')
plt.plot(t, y_robust, label='robust lsq')
plt.xlabel('$t$')
plt.ylabel('$y$')
plt.legend()
Explanation: Veamos los resultados:
End of explanation
N = 100 # number of samples
T = 1./N # sample spacing
t = np.linspace(-1, N*T, N)
y = (np.sin(
    2*np.pi*0.75*t*(1-t) + 2.1) +
    0.1*np.sin(2*np.pi*1.25*t + 1) +
    0.18*np.cos(2*np.pi*3.85*t)
    )
t_exp = (t + 1)
y_exp = y + np.random.randn(len(t)) * 0.30 # ruído
plt.plot(t_exp, y_exp, label='exp data', alpha=0.75)
plt.xlabel('$t$')
plt.ylabel('$y$')
plt.legend()
Explanation: El paquete optimize incluye multitud de métodos para optimización, ajuste de curvas y búsqueda de raíces. La ayuda de este paquete es bastante larga (puedes consultarla también en http://docs.scipy.org/doc/scipy/reference/tutorial/optimize.html).
Procesado de señales (scipy.signal)
Afrontémoslo, lo habitual es que tus datos tengan siempre un ruído asociado y no dispongamos de un modelo al que ajustarlos a priori. Bajo esta situación se suelen optar por diferentes técnicas para filtrar la señal.
Primero generaremos una señal con ruído:
End of explanation
from scipy.signal import medfilt
n_elements = 11 # nº de elementos de en el que se aplica el filtro 
y_exp_filt = medfilt(y_exp, n_elements)
plt.plot(t_exp, y_exp, label='exp data', alpha=0.55)
plt.plot(t_exp, y_exp_filt, label='filt. (median)')
plt.plot(t_exp, y, '-k', label='true', )
plt.xlabel('$t$')
plt.ylabel('$y$')
plt.legend()
Explanation: Lo más sencillo en estos casos es aplicar un filtro por ventana:
End of explanation
from scipy.signal import savgol_filter
n_elements = 11 # nº de elementos de en el que se aplica el filtro 
polyorder = 3
y_exp_filt = savgol_filter(y_exp, n_elements, polyorder)
plt.plot(t_exp, y_exp, label='exp data', alpha=0.55)
plt.plot(t_exp, y_exp_filt, label='filt. (Savitzky-Golay)')
plt.plot(t_exp, y, '-k', label='true', )
plt.xlabel('$t$')
plt.ylabel('$y$')
plt.legend()
Explanation: Como podemos ver, si la frecuencia de muestreo no es muy alta o el tamaño de nuestra ventana no es el adecuado, el resultado puede ser no satisfactorio.
El filtro de Savitzky–Golay suele ser interesante en estos casos.
End of explanation
from scipy import signal
filt_order = 3 # Filter order
Wn = 0.2 # # Cutoff frequency [rad·s^-1]
# Create the filter
b, a = signal.butter(filt_order, Wn, btype='low')
w, h = signal.freqs(b, a)
# Second, apply the filter
y_exp_lowpass = signal.filtfilt(b,a, y_exp)
plt.plot(t_exp, y_exp, label='exp data', alpha=0.55)
plt.plot(t_exp, y_exp_lowpass, label='filt. (low pass)')
plt.plot(t_exp, y, '-k', label='true', )
plt.xlabel('$t$')
plt.ylabel('$y$')
plt.legend()
Explanation: Si estamos interesados en obtener una señal sinuisidal y dado que el ruído ocurre a una frecuencia más alta, otra opción es generar filtro de paso bajo.
End of explanation
from scipy import signal
N = 100
t0 = 0
tf = 100
m = 1
b = 0
t = np.linspace(t0, tf, N)
y = m*t+b
y_exp = y + 100* np.random.rand(N)
plt.plot(t, y_exp, linewidth=2, label='exp data') 
plt.plot(t, signal.detrend(y_exp), linewidth=2, label='drift correction') 
plt.xlabel('$t$')
plt.ylabel('$y$')
plt.legend()
Explanation: Por último, si la señal tiene un deriva (drifting) podemos corregirla fácilmente con:
End of explanation
<END_TASK>