Simulación – Procesos Estocásticos

En un sistema en que puede cambiar el estado, si definimos los estados posibles cómo 1, 2, 3, ... , n. Los estados cómo X1, X2, X3, ... , Xn.

Los cambios de estado sucederán en instantes discretos de tiempo t1, t2, t3, ... Estos procesos de cambio en el sistema se denominan “Cadena de Markov” y satifacen la siguiente propiedad -principal propiedad de una cadena de Markov- El siguiente estado en el proceso depende sólo del estado actual y de ningún otro estado previo.

Las transiciónes entre estados en la cadena de Markov son determinados por una matriz de transición P = || pij ||

Donde

pij son las probabilidades de que el proceso cambie su estado i al estado j -en cualquier instante- en que la transición sea posible.

En la mayoría de cadenas de Markov se puede determinar la distribución final

P{i} significa la probabilidad de que el proceso está en el estado i en un instante arbitrario de tiempo. Las probabilidades estacionarias πi se encuentran mediante el siguiente sistema lineal.

Procesamiento estadístico

Cada simulación del experimento debe dar una cadena de resultados que difiere de otras, entonces que información extraeriamos para analizar y que se podría utilizar para realizar comparaciones?

En la práctica se utiliza la distribución estacionaria, que puede ser obtenida solo después de una larga simulación de un proceso. Se empieza escogiendo un valor T de un periodo de tiempo y suponemos que después de la evolución durante ese periodo, el proceso será casi estacionario y también escogemos un número N que determinará el número de experimentos a realizar.

Algoritmo

Después de un periodo con duración T se mantiene un régimen estacionario de manera que no es necesario realizar la simulación desde el inicio hasta T, se la puede realizar hasta obtener una muestra -números de cambio de estado- igual a N.

  1. Simular el proceso durante un periodo con duración T y así hasta el siguiente cambio de estado, de manera que tenemos un estado i durante un periodo de tiempo t.
  2. Continuar con la simulación del proceso hasta el próximo cambio de estado i_new durante un periodo de tiempo t_new. Así está claro que el proceso estuvo en un estado i durante el periodo de tiempo (t_new - t).
  3. Reasignar i := i_new y t := t_new
  4. Si se han completado N iteraciones, se da por terminado el proceso.

Por cada periodo en que se realizó este proceso hay que dividir la duración para el total durante todas las iteraciones.

Descripción

Simular los cambios en el clima durante un periodo de tiempo con la siguiente matriz de transición.

CLEARCLOUDYOVERCAST
CLEAR-0.40.30.1
CLOUDY0.4-0.80.4
OVERCAST0.10.4-0.5

Resultado

Implementación

enum MyWeather
{
     Clear, Cloudy, Overcast 
}
int i = 0;
int countClear = 0;
int countCloudy = 0;
int countOvercast = 0;

int w = 4;
double sum = 0;

int t = int.Parse(textBoxTime.Text);
int n = int.Parse(textBoxIterations.Text);

// A -- Clear = 0
// B -- Cloudy = 1
// C -- Overcast = 2
double w1 = -0.4; //  Clear->Clear
double w2 = -0.8; //  Cloudy->Cloudy
double w3 = -0.5;  // Overcast->Overcast

double w4 = 0.4;  // Cloudy->Clear   -- Overcast->Cloudy -- Cloudy->Overcast
double w5 = 0.1;  // Overcast->Clear -- Clear->Overcast
double w6 = 0.3;  // Clear->Cloudy

double p1 = 0.33; //in any stage it can start Clear, Cloudy or Overcast
double p2 = 0;
double p3 = 0;
double p4 = 0;

int[] weather = new int[w];
for (int it = 1; it <= n; it++)
{
    labelCurrentIteration.Invoke((MethodInvoker)delegate
    {
        labelCurrentIteration.Text = "Current Iteration: " + it.ToString();
        labelCurrentIteration.Update();
    });

    weather[0] = rand.Next(0, 3);
    weather[1] = rand.Next(0, 3);
    weather[2] = rand.Next(0, 3);
    weather[3] = rand.Next(0, 3);
    
    labelCurrentTime.Invoke((MethodInvoker)delegate
    {
        labelWeather.Text = "Weather for the next hour: " +
                           (MyWeather)weather[0] + " -> " +
                           (MyWeather)weather[1] + " -> " +
                           (MyWeather)weather[2] + " -> " +
                           (MyWeather)weather[3];
        labelWeather.Update();
    });

    for (int wi = 0; wi < weather.Length - 1; wi++)
    {
        if (weather[wi] == 0) // Clear
        {
            countClear++;

            if (weather[wi + 1] == 0) // Clear
            {
                if (wi+1 == 1) p2 = w1;
                else if (wi+1 == 2) p3 = w1;
                else if (wi + 1 == 3) p4 = w1;
            }
            else if (weather[wi + 1] == 1) // Cloudy
            {
                if (wi + 1 == 1) p2 = w6;
                else if (wi + 1 == 2) p3 = w6;
                else if (wi + 1 == 3) p4 = w6;
            }
            else if (weather[wi + 1] == 2) // Overcast
            {
                if (wi + 1 == 1) p2 = w5;
                else if (wi + 1 == 2) p3 = w5;
                else if (wi + 1 == 3) p4 = w5;
            }
        }
        if (weather[wi] == 1) // Cloudy
        {
            countCloudy++;

            if (weather[wi + 1] == 0 || weather[wi + 1] == 2) // Clear OR Overcast
            {
                if (wi + 1 == 1) p2 = w4;
                else if (wi + 1 == 2) p3 = w4;
                else if (wi + 1 == 3) p4 = w4;
            }
            else if (weather[wi + 1] == 1) // Cloudy
            {
                if (wi + 1 == 1) p2 = w2;
                else if (wi + 1 == 2) p3 = w2;
                else if (wi + 1 == 3) p4 = w2;
            }
        }
        if (weather[wi] == 2) // Overcast
        {
            countOvercast++;

            if (weather[wi + 1] == 0) // Clear
            {
                if (wi + 1 == 1) p2 = w5;
                else if (wi + 1 == 2) p3 = w5;
                else if (wi + 1 == 3) p4 = w5;
            }
            else if (weather[wi + 1] == 1) // Cloudy
            {
                if (wi + 1 == 1) p2 = w4;
                else if (wi + 1 == 2) p3 = w4;
                else if (wi + 1 == 3) p4 = w4;
            }
            else if (weather[wi + 1] == 2) // Overcast
            {
                if (wi + 1 == 1) p2 = w3;
                else if (wi + 1 == 2) p3 = w3;
                else if (wi + 1 == 3) p4 = w3;
            }
        }
    }

    sum = p1 + p2 + p3 + p4;
    double[] e_probabilities = new double[9];
    double[] probabilities = new double[9];
    double[] experiments = new double[9];
    double[] frequency = new double[9];

    //Expected Value - µ
    double xp1 = 1 * p1;
    double xp2 = 2 * p2;
    double xp3 = 3 * p3;
    double xp4 = 4 * p4;

    double miu = xp1 + xp2 + xp3 + xp4;

    //Variance - µ
    double x2p1 = (1 * 1) * p1;
    double x2p2 = (2 * 2) * p2;
    double x2p3 = (3 * 3) * p3;
    double x2p4 = (4 * 4) * p4;

    double variance = x2p1 + x2p2 + x2p3 + x2p4;
    variance -= (miu * miu);

    experiments[0] = 0;
    experiments[1] = 0;
    experiments[2] = 0;
    experiments[3] = 0;
    
    e_probabilities[0] = xp1;
    e_probabilities[1] = xp2;
    e_probabilities[2] = xp3;
    e_probabilities[3] = xp4;
   
    probabilities[0] = p1;
    probabilities[1] = p2;
    probabilities[2] = p3;
    probabilities[3] = p4;

    while (i < n)
    {
        int r = rand.Next(0, 3); // Clear - Cloudy - Overcast

        if (r == 1) experiments[0]++;
        else if (r == 2) experiments[1]++;
        else if (r == 3) experiments[2]++;
        else if (r >= 4) experiments[3]++;

        i++;
    }

    i = 0;
    foreach (var item in experiments)
    {
        frequency[i++] = item / n;
    }
    
    double e_miu = 1 * frequency[0] +
                   2 * frequency[1] +
                   3 * frequency[2] +
                   4 * frequency[3];
 
    double e_variance = (1 * 1) * frequency[0] +
                        (2 * 2) * frequency[1] +
                        (3 * 3) * frequency[2] +
                        (4 * 4) * frequency[3];

    e_variance -= (e_miu * e_miu);

    //Errors
    double miuError = Math.Abs(e_miu - miu) / Math.Abs(miu);
    double varianceError = Math.Abs(e_variance - variance) / Math.Abs(variance);

    //Chi
    double chi = 0;
    double alfa = 0.05;

    for (i = 0; i < 4; i++)
    {
        chi += Math.Pow((frequency[i] - probabilities[i]), 2) / probabilities[i];
    }
}

Descargar Proyecto