<jupyter_start><jupyter_text>Omitir por el momento<jupyter_code>def calc_error(res,e,names=[]):
    tmp_i = np.zeros(len(res.x))
    etol = e * max(1, abs(res.fun))
    for i in range(len(res.x)):
        tmp_i[i] = 1.0
        hess_inv_i = res.hess_inv(tmp_i)[i]
        uncertainty_i = np.sqrt(etol * hess_inv_i)
        tmp_i[i] = 0.0
        if len(names) > 0:
            print('{0} = {1:12.4e} ± {2:.1e}'.format(names[i], res.x[i], uncertainty_i))
        else:
            print('x^{0} = {1:12.4e} ± {2:.1e}'.format(i, res.x[i], uncertainty_i))<jupyter_output><empty_output><jupyter_text>Revisar bases de datos conocidas[SciKit Learn](https://scikit-learn.org/stable/datasets/index.html)<jupyter_code>A = pd.DataFrame(housing['data'], columns=housing['feature_names'])
b = pd.DataFrame(housing['target'], columns=["MEDV"])
A.insert(0,'c0',1)
 # -36.9419202
 # -3.69419202 -> numerico de precision
A<jupyter_output><empty_output><jupyter_text>$$\mathbb{A}\cdot\vec{x} + c_0 \cdot \mathbb{1} = \vec{b} $$Tenemos que utilizar minimimos cuadrados?$$MSE=\min_{\vec{x}} || \vec{b} - \mathbb{A}\cdot\vec{x} ||^2$$La funcion de costo -> Cost function$$\vec{x} = (\mathbb{A}^T\cdot\mathbb{A})^{-1}\cdot\mathbb{A}^T\vec{b} $$y = rho()*sigma_y/sigma_x$$X = x_1+x_2+x_3+\dots$$Mi intención,$$CRIM\times x_1+ZN\times x_2+INDUS\times x_3+\dots = b$$<jupyter_code>Anumpy = np.matrix(A.to_numpy(dtype=np.float64)) # ""
bnumpy = np.matrix(b.to_numpy(dtype=np.float64))
# %%timeit
X_opt = np.linalg.inv(Anumpy.T*Anumpy)*Anumpy.T*bnumpy
print(X_opt.T)
np.array(Anumpy*X_opt).flatten()
sns.scatterplot(np.array(Anumpy*X_opt).flatten(),np.array(bnumpy).flatten())
plt.plot([0,10],[0,10])
predicciones = np.array(Anumpy*X_opt).flatten()
observaciones = np.array(bnumpy).flatten()
bbarrita = observaciones.mean()
r2 = 1-(np.linalg.norm(predicciones - observaciones)**2)/(np.linalg.norm(observaciones - bbarrita)**2)
r2<jupyter_output><empty_output><jupyter_text>Implementemos nuestro algoritmo de Optimización por <NAME>La idea es encontrar $X_{opt}$ de forma iterativa sin invertir la matriz<jupyter_code># Definamos la func de costo (debe devolver un escalar)
def L(x,A,b):
    # (b_pred-b_obs)^2
    # m es el numero de datos
    # n es el numero de parametros == 13
    m,n = A.shape
    X = np.matrix(x).T
    DeltaB=(A*X-b) # b gorro - b
    return (DeltaB.T*DeltaB)[0,0]/m # matriz 1x1
def dLdx(x,A,b): # gradiente de L
    # (b_pred-b_obs)^2
    # m es el numero de datos
    # n es el numero de parametros == 13
    m,n = A.shape
    X = np.matrix(x).T
    DeltaB=(A*X-b) 
    return (2/m)*np.array(A.T*DeltaB).flatten() # como un vector [1,2,3] y no como un "vector matriz" [[1],[2],[3]]
# encontrar una forma iterativa de actualizar X para ir minimizando la funcion de costo L
x = np.zeros(Anumpy.shape[1])
epsilon = 2e-8
cost = []
N = 100
for it in range(N):
    x = x - epsilon*dLdx(x,Anumpy,bnumpy)
    cost.append(L(x,Anumpy,bnumpy))
sns.scatterplot(np.arange(N)+1,cost)
# plt.xscale('log')
# plt.yscale('log')
L(np.array(X_opt).flatten(),Anumpy,bnumpy)
np.array(X_opt).flatten()
# %%timeit
res1 = sp.optimize.minimize(fun=L,x0=np.zeros(Anumpy.shape[1]), args = (Anumpy,bnumpy), tol=1e-10)
res1['x']
e=1e-10
# %%timeit
res2 = sp.optimize.minimize(fun=L,jac=dLdx, x0=np.zeros(Anumpy.shape[1]), args = (Anumpy,bnumpy), method='Newton-CG', tol=e)
res2['x']
# %%timeit
res3 = sp.optimize.minimize(fun=L,jac=dLdx, x0=np.zeros(Anumpy.shape[1]), args = (Anumpy,bnumpy), method='L-BFGS-B', tol=e)
res3['x']
L(res1['x'],Anumpy,bnumpy)
L(res2.x,Anumpy,bnumpy)
L(res3.x,Anumpy,bnumpy)
L(np.array(X_opt).flatten(),Anumpy,bnumpy)
calc_error(res3,e,names=[])<jupyter_output>x^0 =  -3.6932e+00 ± 5.1e-02
x^1 =   4.3667e-01 ± 7.6e-04
x^2 =   9.4349e-03 ± 9.0e-05
x^3 =  -1.0752e-01 ± 4.0e-04
x^4 =   6.4586e-01 ± 3.9e-03
x^5 =  -4.0173e-06 ± 4.4e-07
x^6 =  -3.7768e-03 ± 1.8e-05
x^7 =  -4.2132e-01 ± 5.3e-03
x^8 =  -4.3444e-01 ± 5.9e-03<jupyter_text>Fit a funciones no lineales<jupyter_code>print("Nuestra caja experimental!")
a = 3/2
b = 4
c = -3
N=100
x = np.linspace(0.2,10,N) # 0.2 y 10 en N elementos igualmente espaciados
y = a/(1+np.exp(c*(x-b))) # modelo teorico -> principio fisico, matematico, biologico...
x1 = x + np.random.exponential(0.01,size=N)
y1 = y + np.random.normal(0,0.05,size=N) # ruido gaussiano
x2 = x + np.random.normal(0,0.03,size=N)
y2 = y + np.random.exponential(0.05,size=N) # ruido exponencial
sns.scatterplot(x,y)
sns.scatterplot(x,y2)
sns.scatterplot(x1,y1)
sns.scatterplot(x2,y2)<jupyter_output>/opt/conda3/envs/fisi2028/lib/python3.9/site-packages/seaborn/_decorators.py:36: FutureWarning: Pass the following variables as keyword args: x, y. From version 0.12, the only valid positional argument will be `data`, and passing other arguments without an explicit keyword will result in an error or misinterpretation.
  warnings.warn(<jupyter_text>¿Cómo hacer el fit a la función?$$f(x) = a\frac{1}{1+e^{b(x-c)}}$$<jupyter_code>#como encuentro yo a, b y c?
# cual seria nuestra funcion de costo
def f(parametros,x): # parametros es un vector de 3 componentes
    return parametros[0]/(1+np.exp(parametros[1]*(x-parametros[2])))
def Lfit(parametros,x,y): # funcion de costo MSE (No es la mejor!)
    # L = promedio sobre todos los puntos (f(a,b,c;x)-y)^2
    # parametros np.array([a,b,c])
    deltaY=f(parametros,x) - y
    return np.dot(deltaY,deltaY)/len(y)
print("Ajuste para el primer set: x1,y1")
e=1e-8
# ansatz: a=1,b=0,c=0
res1 = sp.optimize.minimize(fun=Lfit, x0=np.array([1,0,0]), args = (x2,y2), method='L-BFGS-B', tol=e)
res1
y1_pred = f(res1.x,x1)
sns.scatterplot(x1,y1)
plt.plot(x1,y1_pred,color='r')
r2 = 1-np.sum((y1_pred-y1)**2)/np.sum((y1-y1.mean())**2)
r2
calc_error(res1,e,names=['a','b','c'])<jupyter_output>a =   1.5596e+00 ± 1.1e-04
b =  -2.7094e+00 ± 7.4e-04
c =   3.9522e+00 ± 2.8e-04<jupyter_text>¿Cómo hacer el fit a la función?$$f(x) = a + b\tanh(c(x-d))$$<jupyter_code>#como encuentro yo a, b y c?
# cual seria nuestra funcion de costo
def ftilde(parametros,x):
    return parametros[0]+parametros[1]*np.tanh(parametros[2]*(x-parametros[3]))
def Lfit(parametros,x,y):
    # L = promedio sobre todos los puntos (f(a,b,c,d;x)-y)^2
    # parametros np.array([a,b,c,d])
    deltaY=ftilde(parametros,x) - y
    return np.dot(deltaY,deltaY)/len(y)
print("Ajuste para el primer set: x1,y1")
e=1e-8
# ansatz: a=0,b=1,c=0,d=0
res1 = sp.optimize.minimize(fun=Lfit, x0=np.array([0,1,0,0]), args = (x1,y1), method='L-BFGS-B', tol=e)
res1
y1_pred = ftilde(res1.x,x1)
sns.scatterplot(x1,y1)
plt.plot(x1,y1_pred,color='r')
r2 = 1-np.sum((y1_pred-y1)**2)/np.sum((y1-y1.mean())**2)
r2
calc_error(res1,e,names=['a','b','c','d'])<jupyter_output>a =   7.4356e-01 ± 8.7e-05
b =   7.4804e-01 ± 9.5e-05
c =   1.5705e+00 ± 6.4e-04
d =   4.0267e+00 ± 3.1e-04
{'fraction_non_alphanumeric': 0.19716646989374262, 'fraction_numerical': 0.06301652892561983, 'mean_word_length': 4.624066390041493, 'pattern_counts': {'":': 0, '<': 26, '<?xml version=': 0, '>': 30, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '373518', 'license': 'MIT', 'n_tokens_mistral': 3310, 'n_tokens_neox': 3142, 'n_words': 552, 'path': 'ejercicios/semana7-8/sistemas-lineales.ipynb', 'repo_name': 'SamuelCanas/FISI2028-202120', 'size': '6784'}