Metodos Numericos: abril 2017

Justificación y propósito

La interpolación es utilizada para hacer predicciones de un valor que se encuentra entre una serie de puntos obtenidos de observaciones experimentales o empíricas. En otras palabras, valores que no pudieron ser medidos u observados en alguna clase de experimento.

El algoritmo de interpolación encuentra una función que cruza todos los puntos dados, y de esta manera, es posible evaluar cualquier valor dentro del valor mínimo y el máximo de dichos puntos para estimar cualquier otro valor. Sin embargo, no significa que no tenga error, de hecho, la función creada no describe el comportamiento real de la información utilizada para la interpolación. Por ejemplo, si se mide la temperatura de cierta región en un lapso de tiempo y se crea una función que, para estos datos, dicha función podrá describir de manera aproximada los datos intermedios, pero los datos fuera del rango pueden no tener sentido o tener un incremento o decremento descontrolado.

Volviendo a algo más específico, la razón por la que se utiliza la interpolación por Newton o Lagrange, es porque estos pueden crear una función con n puntos y, a diferencia de una línea recta (interpolación linear) o una simple parábola (interpolación cuadrática), tiende a generar una mejor precisión en las predicciones de valores intermedios.

Antecedentes

Imaginemos que tenemos 2 puntos:

Podemos dibujar una línea recta entre ellos:

Y elegir un punto arbitrario entre estos puntos:

Así, podemos calcular varias distancias:

Entonces (tomando en cuenta los 2 triángulos que se forman), podemos hacer la siguiente igualdad:

Y resolviendo para f(x):

Explicación gráfica de la interpolación

La interpolación consiste en encontrar una función que cruce por una cierta cantidad de puntos.

En el caso de la interpolación por polinomios de Newton y Lagrange, se puede utilizar cualquier cantidad de puntos. Esto, debido a que estos métodos buscan un polinomio cuyo grado depende de la cantidad de valores conocidos, en otras palabras, utiliza una fórmula generalizada que se puede aplicar sin importar la cantidad de puntos y genera un polinomio que se ajusta a dichos puntos.

El funcionamiento gráfico es exactamente el mismo, una vez que se obtiene el polinomio, este ya se puede utilizar para evaluar el valor intermedio deseado.

Como se muestra en la siguiente imagen, con tres puntos se creó un polinomio de tercer grado que se ajusta a los puntos. Posteriormente, se evalúa en un punto intermedio y se da un valor cercano al real.

Codigo en C++

#include <iostream>
#include <cmath>
#include <vector>

using namespace std;
double EPSILON = 0.00001;       //Error


void interpolacionLineal( double y1, double y2, double x1, double x2, double x ) {

    double y = 0.0;

    y = y1 + ( y2 - y1 ) * ( x - x1 ) / ( x2 - x1 );

    cout << "Valor interpolado es (y): " << y << endl;
}

void interpolacionCuadrtica (  double y0, double y1, double y2, double x0, double x1, double x2 ) {

    double b0 = y0;

    // Diferencia dividida discreta
    double b1 = ( y1 - y0 ) / ( x1 - x0 );

    double b2 = ( ( ( y2 - y1 ) / ( x2 - x1 ) ) - b1 ) / ( x2 - x0);

    double a = b2;
    double b = ( b1 - b2 * x1 - b2 * x0 );
    double c = ( b0 - b1 * x0 + b2 * x0 * x1 );

    cout << "\n++++++++++++++++++++++++++++++" << endl;
    cout << "b0: " << b0 << "\tb1: " << b1 << "\tb2: " << b2 << endl;
    cout << "Funcion cuadratica: " << a << "x^2 + " << b << "x + " << c << endl << endl;
}

/**
    Metodo utilizado en la interpolacion Polinomial de Newton
*/
double divisionDiscreta( double *y, double *x, int n, int i, int j, vector<double> *B ) {

    if ( n == 1 ) {

        if ( i == 0 ) {

            B->push_back( y[i] );
            cout << "f(x0): " << y[i] << endl;
        }

        return y[i];
    } else {
        //cout << "xi: " << i << endl;
        //cout << "xj: " << j << endl;
        double val = ( divisionDiscreta( y, x, n - 1, i, j + 1, B ) - divisionDiscreta( y, x, n - 1, i - 1, j, B ) ) / ( x[i] - x[j] );
        if ( j == 0 ) {
            B->push_back( val );
        }
        cout << "f[ x" << i << ", x" << j << " ]: " << val << endl;
        return val;
    }
}

void interpolacionPolinomialNewton() {
    /* Nota: La suma de los binomios son conmutativas */
    double y[] =
    //{2.31, 3.36, 4.59, 6};
    //{0, 1.386294, 1.791759, 1.609438};
    {0.1353, 7.3890, 20.0855};

    double x[] =
    //{0.1, 0.4, 0.7, 1};
    //{1, 4, 6, 5};
    {-2, 2, 3};

    double x_inc = -1;

    vector<double> B;

    int n = sizeof( y ) / sizeof( y[0] );

    divisionDiscreta( y, x, n, n - 1, 0, &B );

    //
    int j = B.size();
    double v = 0;

    cout << "Polinomio: ";                      //Print Function
    while ( j --> 0 ) {

        cout << B[j];                           //Print Function

        double temp = B[j];
        int k = j;
        while ( k --> 0 ) {

            cout << "(x - " << x[k] << ")";     //Print Function
            temp *= ( x_inc - x[k] );
        }

        if ( j != 0)                            //Print Function
            cout << " + ";                      //Print Function

        v += temp;
    }

    cout << "\nf(x) = " << v << endl;

    // El polinomio es: b0 + b1(x - x0) + b2(x - x0)(x - x1) + b3 (x - x0)(x - x1)(x - x2)
}

/**
    Formula utilizda en la iterpolacion de Lagrange
*/
double L( int i, double x_incognita, double *x, int n ) {

    int j = 0;
    double val = 1;

    while ( j < n ) {

        if ( i != j ) {

            val *= ( x_incognita - x[j] ) / ( x[i] - x[j] );
        }

        j++;
    }

    return val;
}

double f( double *y, double *x, double x_incognita, int n ) {

    int i = 0;

    double val = 0;

    while ( i < n ) {

        val += L( i, x_incognita, x, n ) * y[i];
        i++;
    }

    return val;
}

void InterpolacionLagrange() {

    double y[] =
    //{2.31, 3.36, 4.59, 6};
    {0, 1.386294, 1.791759, 1.609438};

    double x[] =
    //{0.1, 0.4, 0.7, 1};
    {1, 4, 6, 5};

    int n = sizeof( y ) / sizeof( y[0] );
    double x_incognita = 10;

    cout << "y calculada: " << f( y, x, x_incognita, n ) << endl;;
}

int main() {
    double y0, y1, y2, x0, x1, x2;

    // y1 = 86.4;
    // y2 = 92.5;

    // x1 = 1.0;
    // x2 = 1.5;

    // x = 1.25;

    // interpolacionLineal ( y1, y2, x1, x2, x );

    y0 = 0;
    y1 = 1.386294;
    y2 = 1.791759;

    x0 = 1;
    x1 = 4;
    x2 = 6;

    //interpolacionCuadrtica( y0, y1, y2, x0, x1, x2 );

    interpolacionPolinomialNewton();

    //InterpolacionLagrange();

    return 0;
}

Ejemplos

Input

Valor x deseado y f(x) verdadero

Output

f(x) = ln(x)

x	y
1	0
4	1.386294
5	1.791759
6	1.609438

X = 2

f(2) = 0.69315

X = 3

f(3) = 1.09861

Newton:

X = 2

f(x) = 0.628767

X = 3

f(x) = 1.07513

Lagrange:

X = 2

f(x) = 0.628767

X = 3

f(x) = 1.07513

f(x) = log x

x	y
8	0.90309
9	0.9542425
11	1.0413927

X = 10

f(x) = 1

Newton:

X = 10

f(x) = 1.00034

Lagrange:

X = 10

f(x) = 1.00034

f(x) = e^x

x	y
-2	0.1353
2	7.3890
3	20.0855

X = -1

f(x) = 0.3678

X = 0

f(x) = 1

X = 1

f(x) = 2.7182

Newton:

X = -1

f(x) = -4.58112

X = 0

f(x) = -4.94431

X = 1

f(x) = -0.95427

Lagrange:

X = -1

f(x) = -4.58112

X = 0

f(x) = -4.94431

X = 1

f(x) = -0.95427

Ajuste de Curvas por Mínimos cuadrados

Regresion Lineal

Antecedentes

El ajuste de Curvas consiste en encontrar una función, tal que, esta se aproxime a los valores discretos tabulados por mediciones experimentales o bien desde el cálculo con valores específicos.

Es un método que facilita el cálculo de valores utilizados en estadística como la Media aritmética, desviación estándar y varianza, coeficiente de variación, entre otros; a través de una lectura de datos donde se automatiza el cálculo de dichos valores y aproxima un comportamiento de una muestra.

Justificación y propósito

El propósito de el ajuste de curvas es encontrar, como dicho anteriormente, una función que se ajuste a ciertos valores ya obtenidos. Esto, con el objeto de dar una aproximación a valores intermedios entre los valores discretizados así como también dar aproximaciones o intentar predecir valores fuera del rango de una nube de datos.

Esto nos permite modelar mediciones experimentales o casos donde todavía no existen tablas o datos. Además, con la información obtenida de estos métodos, se puede caracterizar a una población con base a una muestra que a su vez, permitiría también hacer análisis estadísticos como el análisis de tendencias o pruebas de hipótesis.

En el caso de la regresión lineal, como su mismo nombre lo dice, busca ajustar y describir una línea recta como función para describir el comportamiento de un conjunto de valores u observaciones experimentales.

Diagrama de Flujo

Codigo en c++

#include <iostream>
#include <fstream>
#include <cmath>

using namespace std;

int main(){
    ifstream p_fileX;
    ifstream p_fileY;
    p_fileX.open("Data/Xdata.dat");
    p_fileY.open("Data/Ydata.dat");

    double XY = 0.0, X = 0.0, Y = 0.0, X2 = 0.0, n = 0.0;
    double a0, a1;      //Coeficiente para la fucion de la recta
    double yprom;       //Promedio de las Ys
    double x[31];       //Valores de x
    double y[31];       //Valores de y
    
    int i = 0;

    if ( p_fileX.is_open() ) {
        while ( !p_fileX.eof() && !p_fileY.eof() ) {
            
            p_fileX >> x[i];
                
            p_fileY >> y[i];
            
            cout << x[i] << ", " << y[i] << endl;
            
            X += x[i];
            Y += y[i];
            XY += x[i] * y[i];
            X2 += x[i] * x[i];
            n += 1.0;
            i++;
        }
        
        a1 = ( n * XY - X * Y ) / ( n * X2 - X * X );
        a0 = Y / n - a1 * X / n;
        cout << "y = (" << a0 << ") + (" << a1 << ")x" << endl;
    }
    p_fileX.close();
    p_fileY.close();
    
    yprom = Y / n;
    double st = 0.0;    //suma total de cuadrados alrededor de la media
    double sr = 0.0;    //suma de cuadrados de residuos, alrededor dela linea de regresion
    double Sxy = 0.0;   //Error de estandard de estimacion
    double r = 0.0;     //Coeficiente de correlacion
    
    //Calcular Sr y St
    for ( int j = 0; j < n; j++ ) {
        sr += pow (  ( y[j] - a0 - a1 * x[j] ), 2 );
        st += pow ( ( y[j] - yprom ), 2 );
    }
    
    //VAraicion estndar St / n-1
    //Coeficiente de Variancion c.v = desviacion estantar / Ypromedio
    
    //Calcular Sxy
    Sxy = sqrt( sr / (n - 2) );
    
    //Calcular coeficiente de correlacion r
    r = sqrt ( ( st - sr ) / st );
    
    //Print
    //determinacion es r al cuadrado
    cout << "Sr: " << sr << "\tSt: " << st << endl;
    cout << "Error estandar: " << Sxy << "\tCorrelacion(r) :" << r << endl;
    
    return 0;
}

Modelo Exponencial

Codigo en c++

#include <iostream>
#include <fstream>
#include <cmath>

using namespace std;

int main(){
    ifstream p_fileX;
    ifstream p_fileY;
    p_fileX.open("Data/Xdata.dat");
    p_fileY.open("Data/Ydata.dat");
    
        double XY = 0.0, X = 0.0, Y = 0.0, X2 = 0.0, n = 0.0;
    double a0, a1, ALPHA, BETA;      //Coeficiente para la fucion de la recta
    double yprom;       //Promedio de las Ys
    double x[31];       //Valores de x
    double y[31];       //Valores de y
    
    int i = 0;

    if ( p_fileX.is_open() ) {
        while ( !p_fileX.eof() && !p_fileY.eof() ) {
            
            p_fileX >> x[i];
                
            p_fileY >> y[i];
            
            cout << x[i] << ", " << y[i] << endl;
            
            y[i] = log ( y[i] );
            
            X += x[i];
            Y += y[i];
            XY += x[i] * y[i];
            X2 += x[i] * x[i];
            n += 1.0;
            i++;
        }
        
        a1 = ( n * XY - X * Y ) / ( n * X2 - X * X );
        a0 = Y / n - a1 * X / n;
        
        ALPHA = exp( a0 );
        BETA = a1;
        cout << "y = " << ALPHA << " * exp^" << BETA << "x" << endl;
    }
    p_fileX.close();
    p_fileY.close();
    
    yprom = Y / n;
    double st = 0.0;    //suma total de cuadrados alrededor de la media
    double sr = 0.0;    //suma de cuadrados de residuos, alrededor dela linea de regresion
    double Sxy = 0.0;   //Error de estandard de estimacion
    double r = 0.0;     //Coeficiente de correlacion
    
    //Calcular Sr y St
    for ( int j = 0; j < n; j++ ) {
        sr += pow ( exp ( y[j] ) - ALPHA * exp ( BETA * x[j] ) , 2 );
        st += pow ( exp ( y[j] ) - yprom, 2 );
    }
    
    //VAraicion estndar St / n-1
    //Coeficiente de Variancion c.v = desviacion estantar / Ypromedio
    
    //Calcular Sxy
    Sxy = sqrt( sr / (n - 2) );
    
    //Calcular coeficiente de correlacion r
    r = sqrt ( ( st - sr ) / st );
    
    
    //Print
    cout << "Sr: " << sr << "\tSt: " << st << endl;
    cout << "Error estandar: " << Sxy << "\tCorrelacion(r) :" << r << endl;
    
    return 0;
}

Modelo de Potencia

Codigo en c++

#include <iostream>
#include <fstream>
#include <cmath>

using namespace std;

int main(){
    ifstream p_fileX;
    ifstream p_fileY;
    p_fileX.open("Data/Xdata.dat");
    p_fileY.open("Data/Ydata.dat");

    double XY = 0.0, X = 0.0, Y = 0.0, X2 = 0.0, n = 0.0;
    double a0, a1, ALPHA, BETA;      //Coeficiente para la fucion de la recta
    double yprom;       //Promedio de las Ys
    double x[31];       //Valores de x
    double y[31];       //Valores de y
    
    int i = 0;

    if ( p_fileX.is_open() ) {
        while ( !p_fileX.eof() && !p_fileY.eof() ) {
            
            p_fileX >> x[i];
                
            p_fileY >> y[i];
            
            cout << x[i] << ", " << y[i] << endl;
            
            y[i] = log10 ( y[i] );
            x[i] = log10 ( x[i] );
            
            X += x[i];
            Y += y[i];
            XY += x[i] * y[i];
            X2 += x[i] * x[i];
            n += 1.0;
            i++;
        }
        
        a1 = ( n * XY - X * Y ) / ( n * X2 - X * X );
        a0 = Y / n - a1 * X / n;
        
        ALPHA = pow( 10,  a0 );
        BETA = a1;
        cout << "y = " << ALPHA << " * x^" << BETA << endl;
    }
    p_fileX.close();
    p_fileY.close();
    
    yprom = Y / n;
    double st = 0.0;    //suma total de cuadrados alrededor de la media
    double sr = 0.0;    //suma de cuadrados de residuos, alrededor dela linea de regresion
    double Sxy = 0.0;   //Error de estandard de estimacion
    double r = 0.0;     //Coeficiente de correlacion
    
    //Calcular Sr y St
    for ( int j = 0; j < n; j++ ) {
        sr += pow ( pow ( 10, y[j] ) - ALPHA * pow ( pow ( 10, x[j]) , BETA ) , 2 );
        st += pow ( pow ( 10, y[j] ) - yprom, 2 );
    }
    
    //VAraicion estndar St / n-1
    //Coeficiente de Variancion c.v = desviacion estantar / Ypromedio
    
    //Calcular Sxy
    Sxy = sqrt( sr / (n - 2) );
    
    //Calcular coeficiente de correlacion r
    r = sqrt ( ( st - sr ) / st );
    
    
    //Print
    cout << "Sr: " << sr << "\tSt: " << st << endl;
    cout << "Error estandar: " << Sxy << "\tCorrelacion(r) :" << r << endl;
    
    return 0;
}

Ejemplos

x	y
1	400
2	100
3	25
4	625
5	1.5625

Regresion Lineal:

Y = 311.875 - 27.1875x

Error estandar: 310.184

Coeficiente de correlación: 0.158015

Modelo Exponencial:

y = 1009.53 * exp^-0.925777x

Error estandar: 348.836

Coeficiente de correlación: 0.581741

Modelo de Potencia:

Y = 484.857 * x ^ -2.13467

Error estandar: 350.141

Coeficiente de correlación: 0.582954

Año	Valor en la cuenta
1996	5000
1997	5400
1998	5800
1999	6300
2000	6800
2001	7300
2002	7900
2003	8600
2004	9300
2005	10000
2006	11000

Regresion Lineal:

y = (-1.16755e+06) + (587.273)x

Error estandar: 262.159

Correlación(r) :0.991946

Modelo Exponencial:

y = 9.06531e-65 * exp^0.0781442x

Error estándar: 54.522

Correlación(r) :0.99998

Modelo de Potencia:

y = 0 * x^156.365

Error estandar: -nan

Correlación(r) :-nan

En este caso, el modelo de potencia falla debido a que a0 (coeficiente) es muy cercano a cero, y la computadora lo reconoce como 0.

x	y
2.5	13
3.5	11
5	8.5
6	8.2
7.5	7
10	6.2
12.5	5.2
15	4.8
17.5	4.6
20	4.3

Regresion Lineal:

y = (11.5696) + (-0.43112)x

Error estandar: 1.35745

Correlación(r) :0.898747

Modelo Exponencial:

y = 12.313 * exp^-0.0595688x

Error estandar: 1.07669

Correlación(r) :0.987185

Modelo de Potencia:

y = 21.1458 * x^-0.54029

Error estandar: 0.204067

Correlación(r) :0.999662

Metodos Numericos

domingo, 23 de abril de 2017

Interpolación por Polinomios de Newton y Lagrange

lunes, 3 de abril de 2017

Regresión Lineal

Antecedentes

Justificación y propósito

Diagrama de Flujo

Codigo en c++

Codigo en c++

Codigo en c++