GDIDP_AMAT_16_Practica5

hace 385 días por etlopez16

1. REGRESIÓN

 

Regresión Lineal

Supongamos que tenemos el siguiente conjunto de datos:

puntos=[(1,1),(2,2.5),(3,2.75),(4,4.15),(5,4.95),(6,6.1),(7,8),(8,7.5),(9,9.1),(10,9.95)] points(puntos,color='purple',size=10,figsize=4) 
       

Vamos a adaptar un modelo de regresión lineal a los mismos. Para ello usaremos la función find_fit de Sage, que busca los parámetros adecuados para ajustar un modelo a un conjunto de datos, mediante el método de mínimos cuadrados.

var('a,b') model(x)=a+b*x f1=find_fit(puntos,model,solution_dict=True) print 'La recta de regresion es', f1[a]+f1[b]*x points(puntos,color='purple',size=20,figsize=4) + plot(model(a=f1[a],b=f1[b]),(x,0,10),color='red',thickness=2) 
       
La recta de regresion es 0.9890909090908853*x + 0.15999999999816839
La recta de regresion es 0.9890909090908853*x + 0.15999999999816839
f1[a]+f1[b]*1 
       
1.1490909090890535
1.1490909090890535

Análisis del error:  

residuos=vector(RDF,10); sumares=0; for i in [0..9]: residuos[i]=f1[a]+f1[b]*puntos[i][0]-puntos[i][1]; sumares=sumares+residuos[i]^2 print 'Residuos: ', residuos; print 'Suma de errores al cuadrado: ',sumares 
       
Residuos:  (0.149090909089, -0.36181818182, 0.377272727271,
-0.0336363636383, 0.155454545453, -0.00545454545652, -0.916363636366,
0.572727272725, -0.0381818181839, 0.100909090907)
Suma de errores al cuadrado:  1.50018181818
Residuos:  (0.149090909089, -0.36181818182, 0.377272727271, -0.0336363636383, 0.155454545453, -0.00545454545652, -0.916363636366, 0.572727272725, -0.0381818181839, 0.100909090907)
Suma de errores al cuadrado:  1.50018181818
 
       

Regresión cuadrática

Veamos cómo funciona la regresión cuadrática:

Ajustamos un modelo cuadrático:

var('a,b,c') model2(x)=a+b*x+c*x^2 f2=find_fit(puntos,model2,solution_dict=True) graf2=points(puntos,color='purple',size=20) +plot(model2(a=f2[a],b=f2[b],c=f2[c]),(x,0,10),color='red',thickness=2) show(graf2,figsize=4) print 'Regresión cuadrática', f2[a]+f2[b]*x+f2[c]*x^2 
       
Regresión cuadrática -0.0075757575779553665*x^2 + 1.0724242424244026*x -
0.006666666668869725
Regresión cuadrática -0.0075757575779553665*x^2 + 1.0724242424244026*x - 0.006666666668869725

Análisis del error:

residuos=vector(RDF,10); sumares=0; for i in [0..9]: residuos[i]=f2[a]+f2[b]*puntos[i][0]+f2[c]*(puntos[i][0])^2-puntos[i][1]; sumares=sumares+residuos[i]^2 print 'Residuos o errores: ', residuos; print 'Suma de errores al cuadrado: ',sumares 
       
Residuos o errores:  (0.0581818181776, -0.392121212132, 0.392424242403,
0.0118181817815, 0.216060606004, 0.0551515150712, -0.870909091018,
0.587878787737, -0.0684848486636, 0.00999999977962)
Suma de errores al cuadrado:  1.46987878788
Residuos o errores:  (0.0581818181776, -0.392121212132, 0.392424242403, 0.0118181817815, 0.216060606004, 0.0551515150712, -0.870909091018, 0.587878787737, -0.0684848486636, 0.00999999977962)
Suma de errores al cuadrado:  1.46987878788

OTROS MODELOS PARA AJUSTAR LOS PUNTOS

var('a,b,c,d') model3(x)=a+b*x+c*x^2+d*x^3 f3=find_fit(puntos,model3,solution_dict=True) graf3=points(puntos,color='purple',size=20) +plot(model3(a=f3[a],b=f3[b],c=f3[c],d=f3[d]),(x,0,10),color='red',thickness=2) show(graf3,figsize=6) print 'Regresión cúbica', f3[a]+f3[b]*x+f3[c]*x^2+f3[d]*x^3 
       
Regresión cúbica -0.0030885780907679017*x^3 + 0.04338578088372602*x^2 +
0.8373834498829563*x + 0.25833333333188
Regresión cúbica -0.0030885780907679017*x^3 + 0.04338578088372602*x^2 + 0.8373834498829563*x + 0.25833333333188
residuos=vector(RDF,10); sumares=0; for i in [0..9]: residuos[i]=f3[a]+f3[b]*puntos[i][0]+f3[c]*(puntos[i][0])^2+f3[d]*(puntos[i][0])^3-puntos[i][1]; sumares=sumares+residuos[i]^2 print 'Residuos o errores: ', residuos; print 'Suma de errores al cuadrado: ',sumares 
       
Residuos o errores:  (0.13601398600779446, -0.41806526809344646,
0.3275641024835494, -0.04562937080582419, 0.19382284349382317,
0.07738927683788788, -0.8134615393182418, 0.6527389264808292,
-0.04254079430950597, -0.06783217023385646)
Suma de errores al cuadrado:  1.44041375291
Residuos o errores:  (0.13601398600779446, -0.41806526809344646, 0.3275641024835494, -0.04562937080582419, 0.19382284349382317, 0.07738927683788788, -0.8134615393182418, 0.6527389264808292, -0.04254079430950597, -0.06783217023385646)
Suma de errores al cuadrado:  1.44041375291
datos=[(91,39.433942),(81,38.341728),(70,35.025399),(60,31.730170),(50,29.552647),(40,27.835548),(30,26.107849),(20,24.203518),(10,23.202388),(0,22.094660)]; points(datos,size=30,figsize=4) 
       

                                
                            

                                

Si ajustamos con un spline cúbico:

s=spline(datos) ss=plot(s,(x,0,91),color="blue") cp=points(datos,size=40,color="blue") show(ss+cp,figsize=4) 
       

                                
                            

                                

Si ajustamos con una función del tipo $a/(b+c e^{-ax+91})$

var('a b c x') assume(a>0, b>0.0001) 
       
model2(x,a,b,c) = a/(b+c*exp(-a*x/91)) 
       
fit=find_fit(datos, model2, parameters = [a, b, c], variables = [x], solution_dict = True) print 'a= ',fit[a]; print 'b= ',fit[b]; print 'c= ',fit[c]; 
       
a=  0.397278660493
b=  -0.00742890820925
c=  0.0257466654999
a=  0.397278660493
b=  -0.00742890820925
c=  0.0257466654999
graf=cp+plot(model2(a=fit[a],b=fit[b],c=fit[c]),(x,0,91),color='red',thickness=2) show(graf, figsize=4) 
       

                                
                            

                                
show(ss+cp+graf,figsize=4) 
       

                                
                            

                                

2. Integración Numérica

Vamos a calcular numéricamente la integral $\int_0^1 1+e^{-x}\sin(4x)\,dx$.

 Método del trapecio

El método del trapecio consiste en calcular la integral como el área del trapecio formado por el eje x y los puntos $(a,f(a))$, $(b,f(b))$. Aplicamos la fórmula:

f(x)=1+exp(-x)*sin(4*x) I1=(f(0.0)+f(1.0))/2.0;I1 
       

Vamos a dibujar la función y el trapecio anterior.

a=plot(f(x),(x,0,1),color='blue',thickness=4,ymin=0,ymax=2) v=([(0,f(0)),(1,f(1))]) b=point2d(v,size=120,color='black') c=line([(1,0),(1,f(1))],color='black',thickness=2) show(a+b+c+plot(f(0)+(x-0)*(f(1)-f(0)),(x,0,1),fill='axis', color='green',fillcolor='green',thickness=4)) #show(a+b+line([(1,0),(1,f(1))],color='black',thickness=4)+line([(0,f(0)),(1,f(1))],color='red',thickness=4)) 
       

 

Método de Simpson

Por el método de Simpson, la integral se aproxima por

$$I=\frac{h}{3}\left(f(a)+4 f(x_1)+f(b) \right)$$

Es decir:

f(x)=1+exp(-x)*sin(4*x) I2=(f(0.0)+4*f(0.5)+f(1.0))/6;I2 
       

Recordemos que consiste en interpolar la función en los extremos y el punto medio y calcular la integral de dicho polinomio. El polinomio de interpolación es:

R = PolynomialRing(RDF, 'x') px=vector([0,0.5,1]);py=vector([f(0),f(0.5),f(1)]); puntos=zip(px,py); Q(x) = R.lagrange_polynomial(puntos);expand(Q) 
       

Dibujamos ahora la función, el polinomio a integrar y el área que define la integral.

a=plot(f(x),(x,0,1),color='blue',thickness=4,ymin=0,ymax=2) v=([(0,f(0)),(0.5,f(0.5)),(1,f(1))]) b=point(v,size=120,color='black') c=line([(1,0),(1,f(1))],color='black',thickness=2) d=plot(Q(x),(x,0,1),color='green',fill='axis',fillcolor='green',thickness=4) show(a+b+c+d) 
       

Veamos que la integral obtenida por el método de Simpon según la fórmula $A=\frac{h}{3}(f(a)+4f(a+h)+f(b))$ coincide con la integral del polinomio de interpolación:

Q.integrate(x,0,1).n() 
       

 Trapecio compuesto

Vamos a calcular la integral mediante el método del trapecio compuesto dividiendo el intervalo en 5 subintervalos.

f(x)=1+exp(-x)*sin(4*x) I3=(f(0.0)+2*f(0.2)+2*f(0.4)+2*f(0.6)+2*f(0.8)+f(1.0))/10;I3 
       

Y representamos la función y el área calculada

a=plot(f(x),(x,0,1),color='blue',thickness=4,ymin=0,ymax=2) v=[(k,f(k)) for k in [0,0.2..1]] b=point2d(v,size=60,color='black') c=line([(1,0),(1,f(1))],color='black',thickness=2) l=sum([plot(f(k)+(x-k)*(f(k+0.2)-f(k))*5,(x,k,k+0.2),fill='axis', color='green',fillcolor='green',thickness=4) for k in [0,0.2.. 0.8]]) a+b+c+l 
       

 Integración con Sage

Podemos calcular la integral exacta utilizando el comando "integrate". Si hacemos f.integrate(x,a,b), nos devuelve el valor (exacto) de la integral de la función $f$ entre $a$ y $b$.

f.integrate(x,0,1).n() 
       

En sage tenemos el comando numerical_integral que nos permite calcular la integral numérica de una función (devuelve la integral numerica y una estimación del error cometido):

numerical_integral(f,0,1)