Tarea 2 Valeria Alvarez Navarro

Preguntas Teóricas

1.- ¿Por qué una matriz equivale a una transformación lineal entre espacios vectoriales?

Una matriz vive en un espacio R mxn (pudiéndose descomponer en n vectores de dimensión Rm), al multiplicar una matriz por un vector de dimensión n (x elemento de Rn), se le está aplicando una transformación lineal que va a llevar a este vector al espacio Rm con lo cual se está llevando una transformación lineal entre espacios vectoriales (Rn a Rm).

2.- ¿ Cuál es el efecto de una transformación lineal de una matriz diagonal y el de una matriz ortogonal?

En el caso de la matriz diagonal, el efecto de una transformación lineal es que los vectores transformados pueden expandirse o contraerse. En el caso de la matriz ortogonal, el efecto de una transformación lineal es que todos los ángulos y longitudes de los vectores transformados permanece sin cambio. Esto es, al multiplicar por una escalar c, los vectores van a conservar el mismo ángulo entre ellos y van a preservar la longitud original cuando la transformación los haya rotado.

3.- ¿Qué es la descomposición en valores singulares de una matriz?

Es una manera de transformar los datos contenidos en una matriz (sobre todo cuando estos tienen una dimensión grande). Esta descomposición existe para todas las matrices sin importar su dimensión (al contrario de la descomposición en valores propios).Para llevar a cabo la descomposición en valores singulares (SVD), se utliza la factorización de la matriz A de mxn en tres matrices: dos matrices ortogonales que contienen los vectores singulares (U y V) y una matriz cuasi diagonal que contiene los valores singulares de A, denominada Sigma (S). Para lograr que la matriz S sea diagonal, se le añaden ceros. Los valores de la diagonial de S son siempre no negativos y las entradas se encuentran ordenadas de mayor a menor. La SVD es un método para obtener la aproximación más cercana de oreden k al rango de la matriz A. La elección de k es un trade off entre precisión y velocidad de procesamiento de los datos. A se reconstruye multiplicando las tres matrices: U, S, VT (V transpuesta)

4.- ¿Qué es diagonalizar una matriz y qué representan los eigenvalores?

Diagonalizar una matriz cuadrada (A), se refiere al proceso por el cual se encuentra una matriz (P) tal que al multiplicar P inversa por A por P, se obtiene una matriz que únicamente tiene elementos en la diagonal y en todas sus demás entradas tiene ceros, la cual se llama matriz Diagonal (D). Las entradas de la matriz diagonal D, son los eigenvalores de A, los cuáles representan las raíces de su polinomio característico det|(A-LamdaI)|=0. Para que la matriz A sea pueda ser diagonalizable sus n vectores deben ser linealmente independientes. Más aun, P debe ser una matriz invertible en la cual sus columnas son los n eigenvectores linealmente independientes de A. Los eigenvalores representan el factor escalar por el cual un eigenvector ha sido multiplicado al llevar a cabo una transformación lineal; los cuales a su vez no son afectados por transformaciones lineales en cuanto a su dirección (permanecen como los originales).

5.- Intuitivamente, ¿Qué son los eigenvectores?

Son los vectores que al aplicar una transformación lineal, no son alterados en su dirección o la transformación sólo hace que se multipliquen por un escalar, lo cual tampoco cambia su dirección. Son los vectores que caracterizan las transformación lineal.

6.- ¿Cómo interpretas la descomposición en valores singulares como una composición de tres tipos de transformaciones lineales simples?

La primera matriz al ser ortogonal, debe llevar a cabo una rotación. Posteriormente al aplicar la matriz diagonal se lleva a cabo una expansión o contracción de ejes. Por último se vuelve a llevar a cabo una rotación ya que la tercera matriz también es ortogonal. Preg. 2 y 3.

7.- ¿Qué relación hay entre la descomposición en valores singulares y la diagonalización?

La diagonalización sólo se puede llevar a cabo en matrices cuadradas y debe existir la matriz inversa (P, pregunta 4), ese proceso se llama también eigen descomposición. En cuando a la SVD, no importa el tamaño de las matrices y no se requiere encontrar P inversa, ya que se puede construir una matriz pseudo.inversa. Adicionalmente, la SVD es más estable y con menos limitaciones para manejar grandes volúmenes de datos.

8.- ¿Cómo se usa la descomposición en valores singulares para dar una aproximación de rango menor a una matriz?

La descomposiciòn SVD de A, proporciona una representación de una matriz en la forma de combinaciones lineales de vectores los cuales están ordenados por importancia (matriz S), quedarse con los "k" más importantes o los primeros es equivalente a dar una aproximaciòn de rango menor (orden k) de la matriz A; donde k<=R. Para llevar a cabo este proceso, se deben quedar las primeras k columnas de la matriz U, los primeros k valores de la matriz S, y los primeros k renglones de la matriz Vt).

9.- Describe el método de minimización por descenso gradiente

El método de minimización por descenso gradiente, es algoritmo de aproximación de mínimos locales de una función. El descenso gradiente inicia con un cojunto de parámetros iniciales y de manera iterativa se mueve hacia aquellos que minimizan la función. La manera en que lleva a cabo las iteraciones es a través de cálculo diferencial, siguiendo la dirección de la pendientes de la superficie creada por la función objetivo en la dirección opuesta a la del gradiente.

10.- Menciona 4 ejemplos de problemas de optimización (dos con restricciones y dos sin restricciones) que te parezcan interesantes como Científico de Datos

  • El algoritmo de asignación de registros que pretenecen a una misma clase (Data Quality)
  • La tasa de interés vs el riesgo de crédito
  • La partición (course classing) de variables tal que la pendiente de la variable dependiente sea máxima
  • El punto de corte de un score de originación de tarjetas de crédito sujeto a mantener una tasa de rechazo menor al 10% y una tasa de incumplimiento menor al 15%.

Parte 2: Aplicaciones de Python

Ejercicio 1

- Importar imagen


In [121]:
import numpy as np
import PIL.Image

In [122]:
im = PIL.Image.open("/Users/valeriaalvarez/Documents/rh.jpeg")
col,row =  im.size
A = np.zeros((row*col, 5))
pixels = im.load()
print(pixels[187,250])
for i in range(col):
    for j in range(row):
        #print("i=%d, j=%d" % (i,j))        
        r,g,b =  pixels[i,j]
        radiohead[i*col + j,:] = r,g,b,i,j


(146, 145, 151)

In [123]:
im


Out[123]:

In [228]:
im.size


Out[228]:
(188, 251)

In [124]:
print(radiohead)


[[ 1.  1.  1.  0.  0.]
 [ 1.  1.  1.  0.  1.]
 [ 1.  1.  1.  0.  2.]
 ..., 
 [ 0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.]]

In [125]:
radiohead.shape


Out[125]:
(47188, 5)

- Realizar la descomposición SVD


In [11]:
import numpy as np
RH = np.linalg

In [199]:
U, s, Vh = RH.svd(radiohead, full_matrices=False)
assert np.allclose(radiohead, np.dot(U, np.dot(np.diag(s), Vh)))

In [200]:
R_H = np.dot(np.dot(U, np.diag(s)), Vh)
print(np.std(radiohead), np.std(R_H), np.std(radiohead - R_H))


70.3299615217 70.3299615217 2.92360175958e-13

In [201]:
U


Out[201]:
array([[ -3.26374927e-05,   6.39795402e-05,   3.05808469e-05,
         -4.61751752e-05,   1.29611575e-05],
       [ -4.10402353e-05,   9.67648940e-06,   8.41018955e-05,
         -2.76703868e-05,   2.01277414e-05],
       [ -4.94429779e-05,  -4.46265614e-05,   1.37622944e-04,
         -9.16559835e-06,   2.72943255e-05],
       ..., 
       [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00],
       [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00],
       [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00]])

In [202]:
s


Out[202]:
array([ 44326.45026186,  14134.75339743,   9745.52184027,    435.22285058,
           96.94693191])

In [203]:
s.shape


Out[203]:
(5,)

In [204]:
Vh


Out[204]:
array([[ -4.80112687e-01,  -4.78211716e-01,  -4.88379795e-01,
         -4.04427117e-01,  -3.72463751e-01],
       [  2.96387116e-01,   2.97135870e-01,   3.10812037e-01,
         -3.71634390e-01,  -7.67560232e-01],
       [  9.79235304e-02,   1.03857283e-01,   9.62454984e-02,
         -8.35646307e-01,   5.21590548e-01],
       [ -3.62141636e-01,  -4.65485891e-01,   8.07531036e-01,
         -2.25492904e-03,   8.05370678e-03],
       [  7.35466453e-01,  -6.74952131e-01,  -5.92577769e-02,
         -4.09279314e-03,   6.94778332e-04]])

In [205]:
import scipy.linalg as sc

In [206]:
S = sc.diagsvd(s, 5, 5)
S


Out[206]:
array([[ 44326.45026186,      0.        ,      0.        ,      0.        ,
             0.        ],
       [     0.        ,  14134.75339743,      0.        ,      0.        ,
             0.        ],
       [     0.        ,      0.        ,   9745.52184027,      0.        ,
             0.        ],
       [     0.        ,      0.        ,      0.        ,    435.22285058,
             0.        ],
       [     0.        ,      0.        ,      0.        ,      0.        ,
            96.94693191]])

- Verificar la descomposición SVD


In [207]:
U @ S @ Vh


Out[207]:
array([[  1.00000000e+00,   1.00000000e+00,   1.00000000e+00,
         -1.98847759e-12,   1.43523071e-12],
       [  1.00000000e+00,   1.00000000e+00,   1.00000000e+00,
         -1.77517576e-12,   1.00000000e+00],
       [  1.00000000e+00,   1.00000000e+00,   1.00000000e+00,
         -3.33564562e-13,   2.00000000e+00],
       ..., 
       [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00],
       [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00],
       [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00]])

In [208]:
# otra forma de hacer la verificación
rd_head=dot(U, np.dot(np.diag(s), Vh))
rd_head


Out[208]:
array([[  1.00000000e+00,   1.00000000e+00,   1.00000000e+00,
         -1.98853310e-12,   1.43534174e-12],
       [  1.00000000e+00,   1.00000000e+00,   1.00000000e+00,
         -1.77514800e-12,   1.00000000e+00],
       [  1.00000000e+00,   1.00000000e+00,   1.00000000e+00,
         -3.33398029e-13,   2.00000000e+00],
       ..., 
       [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00],
       [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00],
       [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00]])

con el método anterior, no logré imprimir la imagen con el imshow

- Hacer la aproximación de orden k


In [195]:
from scipy import *
from pylab import *
img = imread("/Users/valeriaalvarez/Documents/rh.jpeg")[:,:,0]
gray()
figure(1)
imshow(img)


Out[195]:
<matplotlib.image.AxesImage at 0x1162df940>

In [220]:
m,n = img.shape
m,n


Out[220]:
(251, 188)

In [235]:
U,s,Vt = svd(img)

In [226]:
Vt.shape


Out[226]:
(188, 188)

In [224]:
s.shape


Out[224]:
(188,)

In [222]:
S = resize(s,[m,1])*eye(m,n)
S


Out[222]:
array([[ 21793.67836849,      0.        ,      0.        , ...,
             0.        ,      0.        ,      0.        ],
       [     0.        ,   6015.72361847,      0.        , ...,
             0.        ,      0.        ,      0.        ],
       [     0.        ,      0.        ,   3552.98001426, ...,
             0.        ,      0.        ,      0.        ],
       ..., 
       [     0.        ,      0.        ,      0.        , ...,
             0.        ,      0.        ,      0.        ],
       [     0.        ,      0.        ,      0.        , ...,
             0.        ,      0.        ,      0.        ],
       [     0.        ,      0.        ,      0.        , ...,
             0.        ,      0.        ,      0.        ]])

In [223]:
S.shape


Out[223]:
(251, 188)

In [232]:
m,n = img.shape
U,s,Vt = svd(img)
S = resize(s,[m,1])*eye(m,n)
k = 5
imshow(dot(U[:,1:k],dot(S[1:k,1:k],Vt[1:k,:])))
show()



In [233]:
k = 10
imshow(dot(U[:,1:k],dot(S[1:k,1:k],Vt[1:k,:])))
show()



In [234]:
k = 25
imshow(dot(U[:,1:k],dot(S[1:k,1:k],Vt[1:k,:])))
show()


¿Qué tiene que ver este proyecto con compresión de imágenes?

La descomposición de valores singulares tiene una de sus mayores aplicaciones en la compresión de imágenes. Lo que entiendo es que esta aplicación es importante conocerla pues la SVD es el mecanismo para que estas imágenes puedan ser trasmitidas de manera eficiente a través de internet con la cantidad mínima de información necesaria para que la imagen no pierda sus elemetos escenciales y con utilizando menor capacidad de almacenamiento.

En la compresión de imágenes es probablemente, donde se puede ver con mayor claridad el potencial del SVD para el análisis de datos.

Ejercicio 2

Aplicación de pseudoinversa y sistema de ecuaciones

  • Programar una función que dada cualquier matriz devuelva la pseudainversa usando la descomposición SVD.
  • Hacer otra función que resuelva cualquier sistema de ecuaciones de la forma Ax=b usando esta pseudoinversa.

In [241]:
from numpy import *
import numpy as np

def gen_mat(i,j):
    A = floor(random.rand(i,j)*20-0) # se está haciendo una matriz aleatoria de 4x4
    b = floor(random.rand(j,1)*20-0) # este es el vector de resultados b
    return A,b

In [242]:
A,b= gen_mat(4,4)

In [243]:
A


Out[243]:
array([[  3.,   9.,  12.,  17.],
       [  0.,  14.,  13.,  12.],
       [  4.,   1.,   6.,   3.],
       [ 15.,   0.,  13.,   8.]])

In [244]:
b


Out[244]:
array([[ 17.],
       [  8.],
       [  5.],
       [  7.]])

In [245]:
#Esto sólo sirve para matrices cuadradas.
def Inversa(A):
    if((A.shape[0] == A.shape[1])):
        U,s,V=np.linalg.svd(A)
        Inversa = np.dot(np.dot(V.T,linalg.inv(diag(s))),U.T)
        return Inversa
    else:
        return "La Matriz no es cuadrada, calcula la pseudoinversa"

In [246]:
w=Inversa(A)
w


Out[246]:
array([[-0.06877898,  0.07959815, -0.49536321,  0.21251932],
       [-0.11089645,  0.17890263, -0.50656878,  0.1572643 ],
       [-0.0007728 , -0.04404946,  0.6236476 , -0.16615147],
       [ 0.13021638, -0.07766615, -0.08462133, -0.00347759]])

In [193]:
m,n=A.shape
m,n


Out[193]:
(5, 4)

In [209]:
U, s, V = np.linalg.svd(A)

In [211]:
s


Out[211]:
array([ 49.52981067,  13.87115693,   7.16708431,   3.16571681])

In [201]:
U.shape


Out[201]:
(5, 5)

In [202]:
V.shape


Out[202]:
(4, 4)

In [199]:
S=eye(m,n)
S


Out[199]:
array([[ 1.,  0.,  0.,  0.],
       [ 0.,  1.,  0.,  0.],
       [ 0.,  0.,  1.,  0.],
       [ 0.,  0.,  0.,  1.],
       [ 0.,  0.,  0.,  0.]])

In [210]:
S = resize(s,[m,1])*eye(m,n)
S


Out[210]:
array([[ 49.52981067,   0.        ,   0.        ,   0.        ],
       [  0.        ,  13.87115693,   0.        ,   0.        ],
       [  0.        ,   0.        ,   7.16708431,   0.        ],
       [  0.        ,   0.        ,   0.        ,   3.16571681],
       [  0.        ,   0.        ,   0.        ,   0.        ]])

In [200]:
Sigma = np.zeros([U.shape[1],V.shape[0]])
Sigma


Out[200]:
array([[ 0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.]])

In [216]:
def P_Inversa(A):
    import numpy as np
    if isinstance(A, np.ndarray): 
        U, s, V = np.linalg.svd(A)
        m,n=A.shape
        S = resize(s,[m,1])*eye(m,n) 
        for i in range(S.shape[0]):
            for j in range(S.shape[1]):
                if (i == j):
                    if (s[i] == 0):
                        S[i,j] == 0
                    else:
                        S[i,j] = 1/s[i]
                        
        P_Inversa = np.dot(np.dot(V.T,S.T),U.T)
        return(P_Inversa)
    else:
        return "Errores en la especificación"

In [217]:
A_I=P_Inversa(A)
A_I


Out[217]:
array([[-0.10203517, -0.0511652 ,  0.03658153,  0.15545005, -0.06976987],
       [-0.01682481, -0.08537225,  0.05294657,  0.04191888,  0.01896219],
       [ 0.0310797 ,  0.16506386, -0.02903122, -0.15212095,  0.05136191],
       [ 0.10515633, -0.01732195, -0.04166221, -0.02419105,  0.02124916]])

Comprobación de que el proceso anterior sale correctamente.


In [228]:
B = np.linalg.pinv(A)
np.allclose(A, np.dot(A, np.dot(B, A)))
B


Out[228]:
array([[-0.10203517, -0.0511652 ,  0.03658153,  0.15545005, -0.06976987],
       [-0.01682481, -0.08537225,  0.05294657,  0.04191888,  0.01896219],
       [ 0.0310797 ,  0.16506386, -0.02903122, -0.15212095,  0.05136191],
       [ 0.10515633, -0.01732195, -0.04166221, -0.02419105,  0.02124916]])

Para obtener la solución para cualquier matriz, hay que juntar las 2 funciones anteriores


In [254]:
def Solucion(A,b):
    import numpy as np
    if isinstance(A, np.ndarray):
        if isinstance(b, np.ndarray):
            if((A.shape[1] == b.shape[0])):#la matriz y el vector son compatibles
                if((A.shape[0] == A.shape[1])):
                    A_inv=Inversa(A)
                    x_sol = np.dot(A_inv,b)
                else:
                    A_inv=P_Inversa(A)
                    x_sol = np.dot(A_inv.T,b)
                return(x_sol)   
            else:
                return "A y b son incompatibles"
        else:
            return "Problemas con b"
    else:
        return "Problemas con A"

In [255]:
x = Solucion(A,b)
x


Out[255]:
array([[-1.52163833],
       [-1.88601236],
       [ 1.58964451],
       [ 1.14489954]])
  • Jugar con el sistema Ax=b donde A = [[1,1],[0,0]] y b puede tomar distintos valores. a) Observar que pasa si b está en la imagen de A (contestar cuál es la imagen) y si no está (ej. b = [1, 1]). b) Contestar, ¿la solución resultante es única? Si hay más de una solución, investigar que caracteriza a la solución directa. c) Repetir cambiando A = [[1,1],[0, 1e-32]], ¿en este caso, la solución es única? ¿Cambia el valor devuelto de x en cada posible valor de b en el punto anterior?

In [256]:
import numpy as np
A1 =np.array([[1,1],[0,0]])
A1


Out[256]:
array([[1, 1],
       [0, 0]])

In [257]:
b1=np.array([[1],[1]])
b1


Out[257]:
array([[1],
       [1]])

In [258]:
x1 = Solucion(A1,b1)
x1


---------------------------------------------------------------------------
LinAlgError                               Traceback (most recent call last)
<ipython-input-258-2caa02fe4478> in <module>()
----> 1 x1 = Solucion(A1,b1)
      2 x1

<ipython-input-254-33b1b306ee80> in Solucion(A, b)
      5             if((A.shape[1] == b.shape[0])):#la matriz y el vector son compatibles
      6                 if((A.shape[0] == A.shape[1])):
----> 7                     A_inv=Inversa(A)
      8                     x_sol = np.dot(A_inv,b)
      9                 else:

<ipython-input-245-db56ccd9c698> in Inversa(A)
      3     if((A.shape[0] == A.shape[1])):
      4         U,s,V=np.linalg.svd(A)
----> 5         Inversa = np.dot(np.dot(V.T,linalg.inv(diag(s))),U.T)
      6         return Inversa
      7     else:

/Users/valeriaalvarez/anaconda/lib/python3.6/site-packages/numpy/linalg/linalg.py in inv(a)
    524     signature = 'D->D' if isComplexType(t) else 'd->d'
    525     extobj = get_linalg_error_extobj(_raise_linalgerror_singular)
--> 526     ainv = _umath_linalg.inv(a, signature=signature, extobj=extobj)
    527     return wrap(ainv.astype(result_t, copy=False))
    528 

/Users/valeriaalvarez/anaconda/lib/python3.6/site-packages/numpy/linalg/linalg.py in _raise_linalgerror_singular(err, flag)
     88 
     89 def _raise_linalgerror_singular(err, flag):
---> 90     raise LinAlgError("Singular matrix")
     91 
     92 def _raise_linalgerror_nonposdef(err, flag):

LinAlgError: Singular matrix

Al hacer el ejercicio con A1 = [[1,1],[0,0]] y b=([[1],[1]]), se obtiene un error de que el sistema no se puede resolver porque A es una matriz singular, las columnas de A son iguales por lo que su determinante es igual a 0. Para evitar este error y que se resuelva por la pseudoinversa, se mete la condición del determinante en la función.


In [108]:
d=np.linalg.det(A1) 
d


Out[108]:
0.0

In [296]:
def Solucion(A,b):
    import numpy as np
    import warnings
    if isinstance(A, np.ndarray):
        if isinstance(b, np.ndarray):
            if((A.shape[1] == b.shape[0])):#la matriz y el vector son compatibles
                if((A.shape[0] == A.shape[1]) and np.linalg.det(A) != 0):
                    A_inv=Inversa(A)
                    x_sol = np.dot(A_inv,b)
                    return(x_sol)
                elif (np.linalg.det(A)==0):
                    A_inv=P_Inversa(A)
                    x_sol = np.dot(A_inv.T,b)
                    print(x_sol) 
                    print("La Matriz A es singular")
                else:
                    A_inv=P_Inversa(A)
                    x_sol = np.dot(A_inv.T,b)
                    return(x_sol)   
            else:
                return "A y b son incompatibles"
        else:
            return "Problemas con b"
    else:
        return "Problemas con A"

Se vuelve a correr el ejercicio con las modificaciones al código

Se vuelve a hacer la prueba con el vector b=[1,1]


In [299]:
x1 = Solucion(A1,b1)
x1


[[ 1.]
 [ 0.]]
La Matriz A es singular

Se prueba con vectores b diferentes:


In [300]:
b2=np.array([[1],[2]])
x2 = Solucion(A1,b2)
x2


[[ 1.5]
 [ 0. ]]
La Matriz A es singular

In [307]:
b3=np.array([[1],[0]])
x3 = Solucion(A1,b3)
x3


[[ 0.5]
 [ 0. ]]
La Matriz A es singular

Las soluciones parecen ser de la forma [x,0]

La Imagen de A, son todas las posibles combinaciones lineales de las columnas de A.

Cambiar el valor de la posición [2,2] de la Matriz A1 por el valor 1e-32.


In [333]:
A2 =np.array([[1,1],[0,1e-32]])
A2


Out[333]:
array([[  1.00000000e+00,   1.00000000e+00],
       [  0.00000000e+00,   1.00000000e-32]])

In [334]:
x4 = Solucion(A2,b1)
x4


Out[334]:
array([[ -1.00000000e+32],
       [  1.00000000e+32]])

In [335]:
x5 = Solucion(A2,b2)
x5


Out[335]:
array([[ -2.00000000e+32],
       [  2.00000000e+32]])

In [336]:
x6 = Solucion(A2,b3)
x6


Out[336]:
array([[  1.00000000e+00],
       [ -1.11022302e-16]])

In [337]:
b4=np.array([[5],[7]])
x7 = Solucion(A1,b4)
x7


[[ 6.]
 [ 0.]]
La Matriz A es singular

In [338]:
x8 = Solucion(A2,b4)
x8


Out[338]:
array([[ -7.00000000e+32],
       [  7.00000000e+32]])

In [341]:
b5=np.array([[5],[1e-13]])
x9 = Solucion(A1,b5)
x9


[[ 2.5]
 [ 0. ]]
La Matriz A es singular

In [343]:
x10 = Solucion(A2,b5)
x10


Out[343]:
array([[ -1.00000000e+19],
       [  1.00000000e+19]])

Al cambiar la entrada [2,2] de la matriz A original de 0 a 1-e32 y correr para diferentes valores de b, se observa que el resultado esta determinado por la elección de b[2,1] sin importar el valor de b[1,1]. En este caso las soluciones parecen ser de la forma [-x,x] donde x=b[2,1].

En comparación con la matriz original, el valor de b[2,1] ahora no es cero, sino que toma un valor muy grande. El pequeño cambio en la Matriz A hace que las soluciones sean muy diferenetes.

Ejercicio 3

Utilizar la paquetería pandas para trabajar el ajuste por mínimos cuadrados a un conjunto de datos.


In [344]:
import numpy as np
from pandas import DataFrame
import pandas as pd
import os
  • Programar un script que lea el archivo study_vs_sat.csv y lo almacene como un dataframe en pandas.

In [371]:
tabla = pd.read_csv("/Users/valeriaalvarez/Documents/tarea2.csv")
tabla


Out[371]:
study_hours sat_score
0 4 390
1 9 580
2 10 650
3 14 730
4 4 410
5 7 530
6 12 600
7 22 790
8 1 350
9 3 400
10 8 590
11 11 640
12 5 450
13 6 520
14 10 690
15 11 690
16 16 770
17 13 700
18 13 730
19 10 640
  • Plantear un problema de optimización que intente hacer una aproximación de la forma: sat_score ~ alpha + beta*study_hoursi, con un valor para cada individuo

In [378]:
df=pd.DataFrame(tabla)
df


Out[378]:
study_hours sat_score
0 4 390
1 9 580
2 10 650
3 14 730
4 4 410
5 7 530
6 12 600
7 22 790
8 1 350
9 3 400
10 8 590
11 11 640
12 5 450
13 6 520
14 10 690
15 11 690
16 16 770
17 13 700
18 13 730
19 10 640
  • Plantear como un problema de optimización que intente hacer una aproximación de la forma sat_score ~ alpha + beta*study_hours minimizando la suma de los errores de predicción al cuadrado. Pueden consultar este link https://en.wikipedia.org/wiki/Simple_linear_regression ¿Cuál es el gradiente de la función que se quiere optimizar (hint: las variables que queremos optimizar son alpha y beta)?

Antes de hacer las funciones, se prueba con el resultado de las librerías de python.


In [417]:
import statsmodels.formula.api as sm

In [397]:
result = sm.ols(formula="sat_score ~ study_hours", data=df).fit()
print (result.params)


Intercept      353.164879
study_hours     25.326468
dtype: float64

In [430]:
study_hours=tabla["study_hours"]
sat_score=tabla["sat_score"]

In [431]:
Suma_study_hours = sum(study_hours)
Suma_sat_score = sum(sat_score)
Suma_hours_score = sum(study_hours*sat_score) 
Suma_hours_2 = sum(study_hours**2)
Suma_score_2 = sum(sat_score**2)
obs = len(study_hours)

In [432]:
Alpha = Suma_sat_score/obs - (Beta*Suma_study_hours)/obs
Alpha


Out[432]:
353.1648794988852

In [433]:
Beta = (obs*Suma_hours_score - Suma_study_hours*Suma_sat_score)/(obs*Suma_hours_2 - Suma_study_hours**2)
Beta


Out[433]:
25.326467777895743
  • Programar una función que reciba valores de alpha, beta y el vector sat_score y devuelva un vector array de numpy de predicciones alpha + beta*study_hours_i, con un valor por cada individuo

In [434]:
def mco(b0, b1, X):
    rows = len(X)
    sat_score_e1 = np.array([b0 + b1*X[i] for i in range(rows)])
    return(sat_score_e1)

In [435]:
sat_score_e1=mco(Alpha, Beta, study_hours)
sat_score_e1


Out[435]:
array([ 454.47075061,  581.1030895 ,  606.42955728,  707.73542839,
        454.47075061,  530.45015394,  657.08249283,  910.34717061,
        378.49134728,  429.14428283,  555.77662172,  631.75602506,
        479.79721839,  505.12368617,  606.42955728,  631.75602506,
        758.38836395,  682.40896061,  682.40896061,  606.42955728])
  • Definir un numpy array X de dos columnas, la primera con unos en todas sus entradas y la segunda con la variable study_hours. Observen que X[alpha,beta] devuelve alpha + betastudy_hours_i en cada entrada y que entonces el problema se vuelve sat_score ~ X*[alpha,beta].

In [405]:
len(df["study_hours"])
col_b0 = [1 for i in range(len(x))]
X_m = np.array([col_b0,df["study_hours"]])
X_m


Out[405]:
array([[ 1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,
         1,  1,  1],
       [ 4,  9, 10, 14,  4,  7, 12, 22,  1,  3,  8, 11,  5,  6, 10, 11, 16,
        13, 13, 10]])

In [436]:
parametros = np.array([Alpha,Beta])
sat_score_e2=np.dot(X_m.T,parametros)
sat_score_e2


Out[436]:
array([ 454.47075061,  581.1030895 ,  606.42955728,  707.73542839,
        454.47075061,  530.45015394,  657.08249283,  910.34717061,
        378.49134728,  429.14428283,  555.77662172,  631.75602506,
        479.79721839,  505.12368617,  606.42955728,  631.75602506,
        758.38836395,  682.40896061,  682.40896061,  606.42955728])
  • Calcular la pseudoinversa B=(X^TX)^(-1)X^T*sat_score para obtener alpha y beta soluciones

In [437]:
np.linalg.pinv((X_m).dot(X_m.T)).dot(X_m).dot(sat_score)


Out[437]:
array([ 353.1648795 ,   25.32646778])
  • Comparar la solución para obtener los parámetros anteriormente planteada y la que se obtienen al hacel la multiplicación de matrices directa.

La solución es la misma excepto por los decimales que se redondean en la de la pseudoinversa.

  • Avanzado: Usen la libreria matplotlib para visualizar las predicciones con alpha y beta solución contra los valores reales de sat_score.

In [438]:
import matplotlib.pyplot as plt

In [456]:
plt.subplot(223)
plt.scatter(study_hours, sat_score_e1, label="pronosticado")
plt.scatter(study_hours, sat_score, label="observado")
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.ylabel('sat score')
plt.xlabel('study hours')
plt.title('Ajuste Modelo')


Out[456]:
<matplotlib.text.Text at 0x11c3e5d30>

In [457]:
print(plt.show())


None
  • Muy avanzado: Programen el método de descenso gradiente para obtener alpha y beta por vía de un método numérico de optimización. Experimenten con distintos learning rates (tamaños de paso).

In [ ]:


In [ ]: