Medidor de energia con Arduino


El autor de esta idea perteneze a un pueblo de Odisha, India, donde los cortes de energía frecuentes son muy comunes. Obstaculiza la vida de todos. Durante su de infancia, continuar con los estudios después del anochecer fue un verdadero desafío. Debido a este problema diseño un sistema solar para su casa de forma experimental. Uso un panel solar de 10 Watt, 6V para iluminar algunos LED brillantes. Después de enfrentar muchas dificultades, el proyecto tuvo éxito. Luego decidio monitorear el voltaje, la corriente, la potencia y la energía involucradas en el sistema. Esto trajo la idea de diseñar un MEDIDOR DE ENERGÍA. Uso ARDUINO como el corazón de este proyecto porque es muy fácil escribir código en su IDE y hay una gran cantidad de bibliotecas de código abierto disponibles en Internet que se pueden usar de acuerdo con el requisito.

Características:
1-Monitoreo de energía por pantalla LCD
2. Envio de lecturas a través de Internet (carga Xively)
3. Registro de datos en una tarjeta SD

Piezas necesarias:

Piezas necesarias:

Potencia y energía

  1. ARDUINO UNO
  2. ARDUINO ETHERNET SHIELD
  3. LCD DE CARACTERES 16×2
  4. SENSOR DE CORRIENTE ACS 712
  5. RESISTENCIAS (10k, 330ohm)
  6. POTENCIOMETRO 10K
  7. CABLES DE PUENTE
  8. CABLE ETHERNET
  9. PANEL DE PAN

Potencia:
potencia es el producto del voltaje (voltios) y la corriente (Amp).P = VxI La unidad de potencia es Watt o KW

Energía:
Energía es producto de la potencia (Watt) y el tiempo (Hora) E = Pxt
Unidad de energía es Watt Hora o Kilovatios hora (kWh)
De la fórmula anterior, queda claro que para medir la energía necesitamos tres parámetros:
1. Voltaje
2. Corriente
3. Tiempo

Medición de voltaje

Medida de voltaje

El voltaje se mide con la ayuda de un circuito divisor de voltaje.Como el voltaje de entrada del pin analógico ARDUINO está restringido a 5V, diseño el divisor de voltaje de tal manera que el voltaje de salida debería ser inferior a 5V.Labatería utilizada para almacenar la potencia del panel solar es de 6 V, 5,5 Ah, por lo que tenemps que reducir este 6,5 V a un voltaje inferior a 5 V.
Usaremos R1 = 10k y R2 = 10K. El valor de R1 y R2 puede ser menor, pero el problema es que cuando la resistencia es baja, el flujo de corriente más alto lo atraviesa, como resultado, una gran cantidad de energía (P = I ^ 2R) se disipa en forma de calor. Por lo tanto, se pueden elegir diferentes valores de resistencia, pero se debe tener cuidado para minimizar la pérdida de potencia a través de la resistencia.

Vout = R2 / (R1 + R2) * Vbat
Vbat = 6.5 cuando está completamente cargado
R1 = 10k y R2 = 10k
Vout = 10 / (10 + 10) * 6.5 = 3.25v que es inferior a 5v y adecuado para el pin analógico ARDUINO

NOTA
Hemos mostrado que la batería de 9 voltios en el circuito de la placa desnuda es solo por ejemplo para conectar los cables. Se utiliza una batería de plomo-ácido de 6 voltios y 5,5 Ah.

Calibración de voltaje:
Cuando la batería está completamente cargada (6.5v) obtendremos un Vout = 3.25v y un valor más bajo para otro voltaje de batería más bajo.

ARDUINO ADC convierte la señal analógica a la correspondiente aproximación digital.
Cuando el voltaje de la batería es de 6.5v, obtuve 3.25v del divisor de voltaje y sample1 = 696 en el monitor en serie, donde sample1 es el valor ADC corresponde a 3.25v

Para una mejor comprensión, he adjuntado la simulación en tiempo real del circuito 123D para la calibración de medición de voltaje:

3.25v equivalente a 696
1 es equivalente a 3.25 / 696 = 4.669mv
Vout = (4.669 * sample1) / 1000 volt
Voltaje real de la batería = (2 * Vout) volt

CÓDIGO ARDUINO:

// tomando 150 muestras del divisor de voltaje con un intervalo de 2seg y luego promediar los datos de muestras

for (int i = 0; i <150; i ++)
{
sample1 = sample1 + analogRead (A2); // lee el voltaje del circuito divisor
delay (2);
}
muestra1 = muestra1 / 150;
voltaje = 4.669 * 2 * muestra1 / 1000;


Medición de corriente

Medida de corriente
Medida de corriente

Para la medición de corriente utilicé un sensor de corriente de efecto Hall ACS 712 (20 A). Hay diferentes sensores de rango de corriente ACS712 disponibles en el mercado, así que elija según sus necesidades. En el diagrama de la placa de pruebas, he mostrado el LED como una carga, pero la carga real es diferente.
PRINCIPIO DE FUNCIONAMIENTO:

El efecto Hall es la producción de una diferencia de voltaje (el voltaje de Hall) a través de un conductor eléctrico, transversal a una corriente eléctrica en el conductor y un campo magnético perpendicular a la corriente.
Para saber más sobre el sensor de efecto Hall, haga clic aquí.
La hoja de datos del sensor ACS 712 se encuentra aquí.

De la hoja de datos
1. El ACS 712 mide 20 amperios positivos y negativos, correspondientes a la salida analógica 100 mV / A
2. Ninguna corriente de prueba a través del voltaje de salida es VCC / 2 = 5v / 2 = 2.5V

Calibración:
La lectura analógica produce un valor de 0-1023, lo que equivale a 0v a 5v Por lo que la lectura analógica 1 = (5/1024) V = 4.89. Valor mv = (4.89 * Valor de lectura analógica) / 1000 V. Pero según las hojas de datos, la compensación es 2.5V (Cuando la corriente es cero, obtendrá 2.5V de la salida del sensor)
Valor real = (valor-2.5) V Corriente en amperios = real valor * 10

CÓDIGO ARDUINO:

// tomar 150 muestras de sensores con un intervalo de 2 segundos y luego promediar los datos de muestras recolectados
for (int i = 0; i <150; i ++)
{
sample2 + = analogRead (A3); // lee la corriente del sensor
delay (2);
}
muestra2 = muestra2 / 150;
val = (5.0 * muestra2) /1024.0;
actualval = val-2.5; // el voltaje de compensación es 2.5v
amps = actualval * 10;

Medición del tiempo

Para la medición del tiempo no se necesita ningún hardware externo, ya que ARDUINO tiene un temporizador incorporado.

La función millis () devuelve el número de milisegundos desde que la placa Arduino comenzó a ejecutar el programa actual.

CÓDIGO ARDUINO:

milisec largo = milis (); // calcula el tiempo en milisegundos
long time = milisec / 1000; // convierte milisegundos a segundos



Cómo ARDUINO calcula la potencia y la e


totamps = totamps + amperios; // calcular amperios totales
avgamps = totamps / time; // amperios promedio
amphr = (avgamps * tiempo) / 3600; // amperio-hora
vatio = voltaje * amperios; // potencia = voltaje *
energía actual = (vatios * tiempo) / 3600; Watt-sec se convierte nuevamente a Watt-Hr dividiendo 1 hora (3600 segundos)
// energía = (vatio * tiempo) / (1000 * 3600); para leer en kWhA

Salida visual

Todos los resultados se pueden visualizar en el monitor en serie o utilizando una pantalla LCD.
Utilicé una pantalla LCD de 16×2 caracteres para mostrar todos los resultados obtenidos en los pasos anteriores. Para los esquemas, consulte el circuito de la placa de pan que se muestra arriba.

Conecte la pantalla LCD con ARDUINO como se indica a continuación:

LCD -> Arduino
1. VSS -> Arduino GND
2. VDD -> Arduino + 5v
3. VO -> Arduino GND pin + Resistencia o potenciómetro
4. RS -> Arduino pin 8
5. RW -> Pin 7 de Arduino
6. E -> Pin 6 de Arduino
7. D0 -> Arduino – No conectado
8. D1 -> Arduino – No conectado
9. D2 -> Arduino – No conectado
10. D3 -> Arduino – No conectado
11 . D4 -> Pin 5 de Arduino
12. D5 -> Pin 4 de Arduino
13. D6 -> Pin 3 de Arduino
14. D7 -> Pin 2 de Arduino
15. A -> Pin 13 de Arduino + Resistencia (potencia de luz de fondo)
16. K -> Arduino GND (tierra de luz de fondo)

CÓDIGO ARDUINO:

Para monitor serie:

Serial.print ("VOLTAJE:");
Serial.print (voltaje);

Serial.println ("Volt");
Serial.print ("ACTUAL:");
Serial.print (amperios);
Serial.println ("Amperios");
Serial.print ("POWER:");
Serial.print (vatios);
Serial.println ("Watt");
Serial.print ("ENERGÍA CONSUMIDA:");
Serial.print (energía);
Serial.println ("Watt-Hora");
Serial.println (""); // imprime los siguientes conjuntos de parámetros después de un
retraso de línea en blanco (2000);

Para LCD:

Para la pantalla LCD, primero debe importar la biblioteca “LiquidCrystal” en el código aquí

El siguiente código es un formato para mostrar en LCD todo el cálculo de potencia y energía

#include <LiquidCrystal.h>
lcd (8, 7, 6, 5, 4, 3, 2);
int backLight = 9;configuración vacía ()
{
pinMode (luz de fondo, SALIDA); // establece el pin 9 como salida
analogWrite (backLight, 150); // controla la intensidad de la luz de fondo 0-254
lcd.begin (16,2); // columnas filas. tamaño de la pantalla
lcd.clear (); // limpia la pantalla
}
void loop ()
{
lcd.setCursor (16,1); // coloca el cursor fuera del recuento de la pantalla
lcd.print (""); // imprime un carácter vacío
delay (600);////////////////////////////////////////// Imprime potencia y energía en una pantalla LCD / /////////////////////////////////////////////////
lcd.setCursor (1,0); // coloca el cursor en la 1ª columna y la 1ª filalcd.print (vatios);
lcd.print ("W");
lcd.print (voltaje);
lcd.print ("V");
lcd.setCursor (1,1); // coloca el cursor en la 1ª columna y la 2ª fila
lcd.print (energía);
lcd.print ("WH");
lcd.print (amperios);
lcd.print ("A");
}

Carga de datos a Xively.com

Carga de datos a Xively.com
Carga de datos a Xively.com

Elija un nombre de usuario, contraseña, establezca su dirección Recibirá un correo electrónico de confirmación; luego haga clic en el enlace de activación para activar su cuenta. Después de abrir correctamente la cuenta, será redirigido a la página de dispositivos de desarrollo .Haga clic en el cuadro + Agregar dispositivo .Dé un nombre a su dispositivo y una descripción (por ejemplo,MONITOREO DE ENERGÍA). Elija datos privados o públicos (elijo privados)
Haga clic en Agregar dispositivo-Después de agregar el dispositivo, se le redirige a una nueva página donde hay mucha información importante.
ID de producto,secreto de producto, número de serie, código de activación·


Feed ID, Feed URL, API End Point (Feed ID se usa en el código ARDUINO)


Agregar canales (elijo ENERGÍA y POTENCIA, pero puede elegir según su elección) Proporcione la unidad y el símbolo para el parámetro·Agrega tu ubicación·Claves API (utilizadas en el código ARDUINO, evite compartir este número)·Disparadores (hacer ping a una página web cuando ocurre un evento, como cuando el consumo de energía supera un cierto límite)

Código Xively y ARDUINO

Archivos adjuntos

Registro de datos en una tarjeta SD

Para almacenar datos en una tarjeta SD, debe importar la biblioteca SD.


El código para almacenar datos en una tarjeta SD se escribe por separado, ya que no tengo suficiente memoria en mi ARDUINO UNO después de escribir el código para la pantalla LCD y cargar los datos xively.com. Pero estoy tratando de mejorar el código de la versión beta para que un solo código pueda contener todas las funciones (pantalla LCD, carga de datos Xively y almacenamiento de datos en una tarjeta SD).
El código para el registro de datos se adjunta a continuación.


/*
  SD card datalogger
 
 This example shows how to log data from three analog sensors 
 to an SD card using the SD library.
 	
 The circuit:
 * SD card attached to SPI bus as follows:
 ** UNO:  MOSI - pin 11, MISO - pin 12, CLK - pin 13, CS - pin 4 (CS pin can be changed)
  and pin #10 (SS) must be an output
 ** Mega:  MOSI - pin 51, MISO - pin 50, CLK - pin 52, CS - pin 4 (CS pin can be changed)
  and pin #52 (SS) must be an output
 ** Leonardo: Connect to hardware SPI via the ICSP header
 		Pin 4 used here for consistency with other Arduino examples
 
original code was created on 24 Nov 2010 and  modified 9 Apr 2012 by Tom Igoe
I (Debasish Dutta) was again modified for my energy monitoring requirement on 14/01/14
podermos 
 
 This example code is in the public domain.
 	 
 */

#include <SD.h>

// On the Ethernet Shield, CS is pin 4. Note that even if it's not
// used as the CS pin, the hardware CS pin (10 on most Arduino boards,
// 53 on the Mega) must be left as an output or the SD library
// functions will not work.
const int chipSelect = 4;

File dataFile;


float sample1=0; // for voltage
float sample2=0; // for current
float voltage=0.0;
float val; // current callibration
float actualval; // read the actual current from ACS 712
float amps=0.0;
float totamps=0.0; 
float avgamps=0.0;
float amphr=0.0;
float watt=0.0;
float energy=0.0; 

void setup()
{
 // Open serial communications and wait for port to open:
  Serial.begin(9600);
   while (!Serial) {
    ; // wait for serial port to connect. Needed for Leonardo only
  }


  Serial.print("Initializing SD card...");
  // make sure that the default chip select pin is set to
  // output, even if you don't use it:
  pinMode(SS, OUTPUT);
  
  // see if the card is present and can be initialized:
  if (!SD.begin(chipSelect)) {
    Serial.println("Card failed, or not present");
    // don't do anything more:
    while (1) ;
  }
  Serial.println("card initialized.");
  
  // Open up the file we're going to log to!
  dataFile = SD.open("energy.txt", FILE_WRITE);
  if (! dataFile) {
    Serial.println("error opening energy.txt");
    // Wait forever since we cant write data
    while (1) ;
  }
}

void loop()

{
   
 long milisec = millis(); // calculate time in milisec
 long time=milisec/1000; // convert time to sec

 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// 
                      /// taking 150 samples from sensors with a inerval of 2sec and then average the samples data collected
  for(int i=0;i<150;i++)
  {
    sample1+=analogRead(A2);  //read the voltage from the sensor
    sample2+=analogRead(A3); //read the current from sensor
    delay(2);
  }
   sample1=sample1/150; 
   sample2=sample2/150;
   
  ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// 
                                                 /////voltage calculation//////////////////////
  
   voltage=4.669*2*sample1/1000; //  callibration // 3.25 from voltage div is eqv 696 in A0 reading 
                                 // multiply 2 to get actual voltage//  I used two 1k resistor to read 6.36v battery volt
                                 
    //////////////////////////////////////////////// current calculation //////////////////////
 val =(5.0*sample2)/1024.0; 
 actualval =val-2.5; // offset voltage is 2.5v 
 amps =actualval*10;// 10 is multiplied as 100mv/A ( from data sheet )
 totamps=totamps+amps; // total amperes 
 avgamps=totamps/time; // average amperess
 amphr=(avgamps*time);  // ampere hour
 watt =voltage*amps;    // power=voltage*current
 energy=(watt*time)/3600;      // energy = power*timein Watt-sec ///again convert to Watt-Hr by dividing 3600sec
 
  ////////////////////////////////////////////////////////
  // make a string for assembling the data to log:
String dataString = "";
int parameter[4]={voltage,amps,watt,energy}; // here parameters are power,energy,watt-hour and current
  // read 4 parameters and append to the string:
  for (int i = 0; i < 4; i++) 
  {
    int sensor = parameter[i];
    dataString += String(sensor);
    if (i < 4) 
    {
      dataString += ","; 
    }
  }

  dataFile.println(dataString);

  // print to the serial port too:
  Serial.println(dataString);
  
  // The following line will 'save' the file to the SD card after every
  // line of data - this will use more power and slow down how much data
  // you can read but it's safer! 
  // If you want to speed up the system, remove the call to flush() and it
  // will save the file only every 512 bytes - every time a sector on the 
  // SD card is filled with data.
  dataFile.flush();
  
  // Take 1 measurement every 500 milliseconds
  delay(500);
}





Fuente https://www.instructables.com/ARDUINO-ENERGY-METER/







Estimación de máxima verosimilitud en Python


¿Qué pasa si una relación lineal no es un supuesto apropiado para nuestro modelo?

Una alternativa ampliamente propuesta por Thomas J. Sargent y John Stachurski es la estimación de máxima verosimilitud, que implica especificar una clase de distribuciones, indexadas por parámetros desconocidos, y luego usar los datos para precisar los valores de estos parámetros.

El beneficio relativo a la regresión lineal es que permite una mayor flexibilidad en las relaciones probabilísticas entre variables.

Aquí ilustramos la máxima probabilidad replicando el artículo de Daniel Treisman (2016) sobre los multimillonarios de Rusia , que conecta la cantidad de multimillonarios en un país con sus características económicas concluyendo que Rusia tiene una mayor cantidad de multimillonarios de lo que predicen factores económicos como el tamaño del mercado y la tasa impositiva.

Requeriremos las siguientes importaciones:

%matplotlib inline
import matplotlib.pyplot as plt
plt.rcParams["figure.figsize"] = (11, 5)  #set default figure size
import numpy as np
from numpy import exp
from scipy.special import factorial
import pandas as pd
from mpl_toolkits.mplot3d import Axes3D
import statsmodels.api as sm
from statsmodels.api import Poisson
from scipy import stats
from scipy.stats import norm
from statsmodels.iolib.summary2 import summary_col

 Configuración y suposiciones 

Consideremos los pasos que debemos seguir en la estimación de máxima verosimilitud y cómo pertenecen a este estudio.

Flujo de ideas 

El primer paso con la estimación de máxima verosimilitud es elegir la distribución de probabilidad que se cree que genera los datos.

Más precisamente, necesitamos suponer qué clase de distribución paramétrica está generando los datos.

  • por ejemplo, la clase de todas las distribuciones normales o la clase de todas las distribuciones gamma.

Cada una de estas clases es una familia de distribuciones indexadas por un número finito de parámetros.

  • Por ejemplo, la clase de distribuciones normales es una familia de distribuciones indexadas por su media. μ∈(−∞,∞) y desviación estándar σ∈(0,∞).

Dejaremos que los datos seleccionen un elemento particular de la clase definiendo los parámetros.

Las estimaciones de parámetros así producidas se denominarán estimaciones de máxima verosimilitud .

 Contando multimillonarios 

Treisman está interesado en estimar el número de multimillonarios en diferentes países.El número de multimillonarios se valora en números enteros.Por lo tanto, consideramos distribuciones que toman valores solo en los enteros no negativos.

(Esta es una de las razones por las que la regresión por mínimos cuadrados no es la mejor herramienta para el problema actual, ya que la variable dependiente en la regresión lineal no está restringida a valores enteros)

Una distribución entera es la distribución de Poisson , cuya función de masa de probabilidad (pmf) esf(y)=μyy!e−μ,y=0,1,2,…,∞

Podemos trazar la distribución de Poisson sobre y para diferentes valores de μ como sigue

poisson_pmf = lambda y, μ: μ**y / factorial(y) * exp(-μ)
y_values = range(0, 25)

fig, ax = plt.subplots(figsize=(12, 8))

for μ in [1, 5, 10]:
    distribution = []
    for y_i in y_values:
        distribution.append(poisson_pmf(y_i, μ))
    ax.plot(y_values,
            distribution,
            label=f'$\mu$={μ}',
            alpha=0.5,
            marker='o',
            markersize=8)

ax.grid()
ax.set_xlabel('$y$', fontsize=14)
ax.set_ylabel('$f(y \mid \mu)$', fontsize=14)
ax.axis(xmin=0, ymin=0)
ax.legend(fontsize=14)

plt.show()

_images / mle_3_0.png

Observe que la distribución de Poisson comienza a parecerse a una distribución normal como la media de y aumenta.

Echemos un vistazo a la distribución de los datos con los que trabajaremos en esta conferencia.

La principal fuente de datos de Treisman son las clasificaciones anuales de multimillonarios de Forbes y su patrimonio neto estimado.

El conjunto de datos mle/fp.dtase puede descargar desde aquí o desde su página AER .

pd.options.display.max_columns = 10

# Load in data and view
df = pd.read_stata('https://github.com/QuantEcon/lecture-python/blob/master/source/_static/lecture_specific/mle/fp.dta?raw=true')
df.head()

paísccodeañoañoentumecidotopint08rintrnoyrsroflawnrrents
0Estados Unidos2.01990.021990.0Yaya39.7999994.98840520,01,61Yaya
1Estados Unidos2.01991.021991.0Yaya39.7999994.98840520,01,61Yaya
2Estados Unidos2.01992,021992.0Yaya39.7999994.98840520,01,61Yaya
3Estados Unidos2.01993.021993.0Yaya39.7999994.98840520,01,61Yaya
4Estados Unidos2.01994.021994.0Yaya39.7999994.98840520,01,61Yaya

5 filas × 36 columnas

Usando un histograma, podemos ver la distribución del número de multimillonarios por país numbil0, en 2008 (los Estados Unidos se descartan para fines de trazado)

numbil0_2008 = df[(df['year'] == 2008) & (
    df['country'] != 'United States')].loc[:, 'numbil0']

plt.subplots(figsize=(12, 8))
plt.hist(numbil0_2008, bins=30)
plt.xlim(left=0)
plt.grid()
plt.xlabel('Number of billionaires in 2008')
plt.ylabel('Count')
plt.show()

_images / mle_7_0.png

A partir del histograma, parece que la suposición de Poisson no es irrazonable (aunque con un valor muy bajo μ y algunos valores atípicos).

Distribuciones condicionales 

En el artículo de Treisman, la variable dependiente: el número de multimillonarios yi en el pais i – se modela en función del PIB per cápita, el tamaño de la población y los años de pertenencia al GATT y la OMC.

Por tanto, la distribución de yi necesita estar condicionado al vector de variables explicativas xi.

La formulación estándar, el llamado modelo de regresión de Poisson , es la siguiente:(58.1)f(yi∣xi)=μiyiyi!e−μi;yi=0,1,2,…,∞.where μi=exp⁡(xi′β)=exp⁡(β0+β1xi1+…+βkxik)

Para ilustrar la idea de que la distribución de yi depende de xi ejecutemos una simulación simple.

Usamos nuestra poisson_pmffunción de arriba y valores arbitrarios para β y xi

y_values = range(0, 20)

# Define a parameter vector with estimates
β = np.array([0.26, 0.18, 0.25, -0.1, -0.22])

# Create some observations X
datasets = [np.array([0, 1, 1, 1, 2]),
            np.array([2, 3, 2, 4, 0]),
            np.array([3, 4, 5, 3, 2]),
            np.array([6, 5, 4, 4, 7])]


fig, ax = plt.subplots(figsize=(12, 8))

for X in datasets:
    μ = exp(X @ β)
    distribution = []
    for y_i in y_values:
        distribution.append(poisson_pmf(y_i, μ))
    ax.plot(y_values,
            distribution,
            label=f'$\mu_i$={μ:.1}',
            marker='o',
            markersize=8,
            alpha=0.5)

ax.grid()
ax.legend()
ax.set_xlabel('$y \mid x_i$')
ax.set_ylabel(r'$f(y \mid x_i; \beta )$')
ax.axis(xmin=0, ymin=0)
plt.show()

_images / mle_9_0.png

Podemos ver que la distribución de yi está condicionado a xi (μi ya no es constante).

Estimación de máxima verosimilitud 

En nuestro modelo para el número de multimillonarios, la distribución condicional contiene 4 (k=4) parámetros que necesitamos estimar.

Etiquetaremos todo nuestro vector de parámetros como β dóndeβ=[β0β1β2β3]

Para estimar el modelo usando MLE, queremos maximizar la probabilidad de que nuestra estimación β^ es el verdadero parámetro β.

Intuitivamente, queremos encontrar el β^ que mejor se adapte a nuestros datos.

Primero, necesitamos construir la función de verosimilitud L(β), que es similar a una función de densidad de probabilidad conjunta.

Supongamos que tenemos algunos datos yi={y1,y2} y yi∼f(yi).

Si y1 y y2 son independientes, la PMF conjunta de estos datos es f(y1,y2)=f(y1)⋅f(y2).

Si yi sigue una distribución de Poisson con λ=7, podemos visualizar el PMF conjunto así

def plot_joint_poisson(μ=7, y_n=20):
    yi_values = np.arange(0, y_n, 1)

    # Create coordinate points of X and Y
    X, Y = np.meshgrid(yi_values, yi_values)

    # Multiply distributions together
    Z = poisson_pmf(X, μ) * poisson_pmf(Y, μ)

    fig = plt.figure(figsize=(12, 8))
    ax = fig.add_subplot(111, projection='3d')
    ax.plot_surface(X, Y, Z.T, cmap='terrain', alpha=0.6)
    ax.scatter(X, Y, Z.T, color='black', alpha=0.5, linewidths=1)
    ax.set(xlabel='$y_1$', ylabel='$y_2$')
    ax.set_zlabel('$f(y_1, y_2)$', labelpad=10)
    plt.show()

plot_joint_poisson(μ=7, y_n=20)

_images / mle_11_0.png

De manera similar, la pmf conjunta de nuestros datos (que se distribuye como una distribución condicional de Poisson) se puede escribir comof(y1,y2,…,yn∣x1,x2,…,xn;β)=∏i=1nμiyiyi!e−μi

yi está condicionado a ambos valores de xi y los parámetros β.

La función de verosimilitud es la misma que la pmf conjunta, pero trata el parámetro β como una variable aleatoria y toma las observaciones (yi,xi) como se indicaL(β∣y1,y2,…,yn ; x1,x2,…,xn)=∏i=1nμiyiyi!e−μi=f(y1,y2,…,yn∣ x1,x2,…,xn;β)

Ahora que tenemos nuestra función de verosimilitud, queremos encontrar el β^ que produce el valor de máxima verosimilitudmaxβL(β)

Al hacerlo, generalmente es más fácil maximizar la probabilidad logarítmica (considere diferenciar f(x)=xexp⁡(x) vs. f(x)=log⁡(x)+x).

Dado que tomar un logaritmo es una transformación creciente monótona, un maximizador de la función de verosimilitud también será un maximizador de la función logarítmica de verosimilitud.

En nuestro caso, la probabilidad logarítmica eslog⁡L(β)= log⁡(f(y1;β)⋅f(y2;β)⋅…⋅f(yn;β))=∑i=1nlog⁡f(yi;β)=∑i=1nlog⁡(μiyiyi!e−μi)=∑i=1nyilog⁡μi−∑i=1nμi−∑i=1nlog⁡y!

El MLE del Poisson al Poisson para β^ se puede obtener resolviendomaxβ(∑i=1nyilog⁡μi−∑i=1nμi−∑i=1nlog⁡y!)

Sin embargo, no existe una solución analítica para el problema anterior: para encontrar el MLE necesitamos usar métodos numéricos.

 MLE con métodos numéricos 

Muchas distribuciones no tienen buenas soluciones analíticas y, por lo tanto, requieren métodos numéricos para resolver las estimaciones de los parámetros.

Uno de esos métodos numéricos es el algoritmo de Newton-Raphson.

Nuestro objetivo es encontrar la estimación de máxima verosimilitud β^.

A β^, la primera derivada de la función logarítmica de verosimilitud será igual a 0.

Ilustremos esto suponiendolog⁡L(β)=−(β−10)2−10

β = np.linspace(1, 20)
logL = -(β - 10) ** 2 - 10
dlogL = -2 * β + 20

fig, (ax1, ax2) = plt.subplots(2, sharex=True, figsize=(12, 8))

ax1.plot(β, logL, lw=2)
ax2.plot(β, dlogL, lw=2)

ax1.set_ylabel(r'$log \mathcal{L(\beta)}$',
               rotation=0,
               labelpad=35,
               fontsize=15)
ax2.set_ylabel(r'$\frac{dlog \mathcal{L(\beta)}}{d \beta}$ ',
               rotation=0,
               labelpad=35,
               fontsize=19)
ax2.set_xlabel(r'$\beta$', fontsize=15)
ax1.grid(), ax2.grid()
plt.axhline(c='black')
plt.show()

_images / mle_13_0.png

La gráfica muestra que el valor de máxima verosimilitud (la gráfica superior) ocurre cuando dlog⁡L(β)dβ=0 (la trama inferior).

Por lo tanto, la probabilidad se maximiza cuando β=10.

También podemos asegurarnos de que este valor sea un máximo (en lugar de un mínimo) comprobando que la segunda derivada (pendiente del gráfico inferior) sea negativa.

El algoritmo de Newton-Raphson encuentra un punto donde la primera derivada es 0.

Para usar el algoritmo, hacemos una estimación inicial del valor máximo, β0 (las estimaciones de los parámetros de MCO pueden ser una suposición razonable), entonces

  1. Utilice la regla de actualización para iterar el algoritmoβ(k+1)=β(k)−H−1(β(k))G(β(k))dónde:G(β(k))=dlog⁡L(β(k))dβ(k)H(β(k))=d2log⁡L(β(k))dβ(k)dβ(k)′
  2. Compruebe si β(k+1)−β(k)<tol
    • Si es verdadero, deje de iterar y configure β^=β(k+1)
    • Si es falso, actualice β(k+1)

Como puede verse en la ecuación de actualización, β(k+1)=β(k) sólo cuando G(β(k))=0es decir. donde la primera derivada es igual a 0.

(En la práctica, dejamos de iterar cuando la diferencia está por debajo de un pequeño umbral de tolerancia)

Intentemos implementar el algoritmo de Newton-Raphson.

Primero, crearemos una clase llamada PoissonRegressionpara que podamos volver a calcular fácilmente los valores de probabilidad de registro, gradiente y hessiano para cada iteración.

class PoissonRegression:

    def __init__(self, y, X, β):
        self.X = X
        self.n, self.k = X.shape
        # Reshape y as a n_by_1 column vector
        self.y = y.reshape(self.n,1)
        # Reshape β as a k_by_1 column vector
        self.β = β.reshape(self.k,1)

    def μ(self):
        return np.exp(self.X @ self.β)

    def logL(self):
        y = self.y
        μ = self.μ()
        return np.sum(y * np.log(μ) - μ - np.log(factorial(y)))

    def G(self):
        y = self.y
        μ = self.μ()
        return X.T @ (y - μ)

    def H(self):
        X = self.X
        μ = self.μ()
        return -(X.T @ (μ * X))

Nuestra función newton_raphsontomará un PoissonRegressionobjeto que tiene una suposición inicial del vector de parámetrosβ0.

El algoritmo actualizará el vector de parámetros de acuerdo con la regla de actualización y volverá a calcular el gradiente y las matrices hessianas en las estimaciones de los nuevos parámetros.

La iteración terminará cuando:

  • La diferencia entre el parámetro y el parámetro actualizado está por debajo de un nivel de tolerancia.
  • Se ha alcanzado el número máximo de iteraciones (lo que significa que no se logra la convergencia).

Para que podamos tener una idea de lo que sucede mientras se ejecuta el algoritmo, display=Truese agrega una opción para imprimir valores en cada iteración.

def newton_raphson(model, tol=1e-3, max_iter=1000, display=True):

    i = 0
    error = 100  # Initial error value

    # Print header of output
    if display:
        header = f'{"Iteration_k":<13}{"Log-likelihood":<16}{"θ":<60}'
        print(header)
        print("-" * len(header))

    # While loop runs while any value in error is greater
    # than the tolerance until max iterations are reached
    while np.any(error > tol) and i < max_iter:
        H, G = model.H(), model.G()
        β_new = model.β - (np.linalg.inv(H) @ G)
        error = β_new - model.β
        model.β = β_new

        # Print iterations
        if display:
            β_list = [f'{t:.3}' for t in list(model.β.flatten())]
            update = f'{i:<13}{model.logL():<16.8}{β_list}'
            print(update)

        i += 1

    print(f'Number of iterations: {i}')
    print(f'β_hat = {model.β.flatten()}')

    # Return a flat array for β (instead of a k_by_1 column vector)
    return model.β.flatten()

Probemos nuestro algoritmo con un pequeño conjunto de datos de 5 observaciones y 3 variables en X.

X = np.array([[1, 2, 5],
              [1, 1, 3],
              [1, 4, 2],
              [1, 5, 2],
              [1, 3, 1]])

y = np.array([1, 0, 1, 1, 0])

# Take a guess at initial βs
init_β = np.array([0.1, 0.1, 0.1])

# Create an object with Poisson model values
poi = PoissonRegression(y, X, β=init_β)

# Use newton_raphson to find the MLE
β_hat = newton_raphson(poi, display=True)

Iteration_k  Log-likelihood  θ                                                           
-----------------------------------------------------------------------------------------
0            -4.3447622      ['-1.49', '0.265', '0.244']
1            -3.5742413      ['-3.38', '0.528', '0.474']
2            -3.3999526      ['-5.06', '0.782', '0.702']
3            -3.3788646      ['-5.92', '0.909', '0.82']
4            -3.3783559      ['-6.07', '0.933', '0.843']
5            -3.3783555      ['-6.08', '0.933', '0.843']
Number of iterations: 6
β_hat = [-6.07848205  0.93340226  0.84329625]

Como se trataba de un modelo simple con pocas observaciones, el algoritmo logró la convergencia en solo 6 iteraciones.

Puede ver que con cada iteración, el valor de probabilidad logarítmica aumentó.

Recuerde, nuestro objetivo era maximizar la función de probabilidad logarítmica, que el algoritmo ha trabajado para lograr.

Además, tenga en cuenta que el aumento en log⁡L(β(k)) se vuelve más pequeño con cada iteración.

Esto se debe a que el gradiente se acerca a 0 a medida que alcanzamos el máximo y, por lo tanto, el numerador de nuestra ecuación de actualización se hace más pequeño.

El vector de gradiente debe estar cerca de 0 en β^

poi.G()

array([[-3.95169228e-07],
       [-1.00114805e-06],
       [-7.73114562e-07]])

El proceso iterativo se puede visualizar en el siguiente diagrama, donde el máximo se encuentra en β=10

logL = lambda x: -(x - 10) ** 2 - 10

def find_tangent(β, a=0.01):
    y1 = logL(β)
    y2 = logL(β+a)
    x = np.array([[β, 1], [β+a, 1]])
    m, c = np.linalg.lstsq(x, np.array([y1, y2]), rcond=None)[0]
    return m, c

β = np.linspace(2, 18)
fig, ax = plt.subplots(figsize=(12, 8))
ax.plot(β, logL(β), lw=2, c='black')

for β in [7, 8.5, 9.5, 10]:
    β_line = np.linspace(β-2, β+2)
    m, c = find_tangent(β)
    y = m * β_line + c
    ax.plot(β_line, y, '-', c='purple', alpha=0.8)
    ax.text(β+2.05, y[-1], f'$G({β}) = {abs(m):.0f}$', fontsize=12)
    ax.vlines(β, -24, logL(β), linestyles='--', alpha=0.5)
    ax.hlines(logL(β), 6, β, linestyles='--', alpha=0.5)

ax.set(ylim=(-24, -4), xlim=(6, 13))
ax.set_xlabel(r'$\beta$', fontsize=15)
ax.set_ylabel(r'$log \mathcal{L(\beta)}$',
               rotation=0,
               labelpad=25,
               fontsize=15)
ax.grid(alpha=0.3)
plt.show()

_images / mle_23_0.png

Tenga en cuenta que nuestra implementación del algoritmo de Newton-Raphson es bastante básica; para implementaciones más sólidas, consulte, por ejemplo, scipy.optimize .

Estimación de máxima verosimilitud con statsmodels

Ahora que sabemos lo que está pasando bajo el capó, podemos aplicar MLE a una aplicación interesante.

Usaremos el modelo de regresión de Poisson statsmodelspara obtener una salida más rica con errores estándar, valores de prueba y más.

statsmodels utiliza el mismo algoritmo anterior para encontrar las estimaciones de máxima verosimilitud.

Antes de comenzar, volvamos a estimar nuestro modelo simple con statsmodels para confirmar que obtenemos los mismos coeficientes y valor logarítmico de verosimilitud.

X = np.array([[1, 2, 5],
              [1, 1, 3],
              [1, 4, 2],
              [1, 5, 2],
              [1, 3, 1]])

y = np.array([1, 0, 1, 1, 0])

stats_poisson = Poisson(y, X).fit()
print(stats_poisson.summary())

Optimization terminated successfully.
         Current function value: 0.675671
         Iterations 7
                          Poisson Regression Results                          
==============================================================================
Dep. Variable:                      y   No. Observations:                    5
Model:                        Poisson   Df Residuals:                        2
Method:                           MLE   Df Model:                            2
Date:                Thu, 15 Jul 2021   Pseudo R-squ.:                  0.2546
Time:                        00:14:01   Log-Likelihood:                -3.3784
converged:                       True   LL-Null:                       -4.5325
Covariance Type:            nonrobust   LLR p-value:                    0.3153
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const         -6.0785      5.279     -1.151      0.250     -16.425       4.268
x1             0.9334      0.829      1.126      0.260      -0.691       2.558
x2             0.8433      0.798      1.057      0.291      -0.720       2.407
==============================================================================

Ahora repliquemos los resultados del artículo de Daniel Treisman, Los multimillonarios de Rusia , mencionado anteriormente en la conferencia.

Treisman comienza estimando la ecuación (58.1) , donde:

  • yi es number of billionairesi
  • xi1 es log⁡GDP per capitai
  • xi2 es log⁡populationi
  • xi3 es years in GATTi – años de membresía en el GATT y la OMC (para poder acceder a los mercados internacionales)

El documento solo considera el año 2008 para la estimación.

Configuraremos nuestras variables para la estimación de esa manera (debe tener los datos asignados dfdesde antes en la lección)

# Keep only year 2008
df = df[df['year'] == 2008]

# Add a constant
df['const'] = 1

# Variable sets
reg1 = ['const', 'lngdppc', 'lnpop', 'gattwto08']
reg2 = ['const', 'lngdppc', 'lnpop',
        'gattwto08', 'lnmcap08', 'rintr', 'topint08']
reg3 = ['const', 'lngdppc', 'lnpop', 'gattwto08', 'lnmcap08',
        'rintr', 'topint08', 'nrrents', 'roflaw']

Entonces podemos usar la Poissonfunción de statsmodelspara ajustar el modelo.

Usaremos errores estándar robustos como en el artículo del autor.

# Specify model
poisson_reg = sm.Poisson(df[['numbil0']], df[reg1],
                         missing='drop').fit(cov_type='HC0')
print(poisson_reg.summary())

Optimization terminated successfully.
         Current function value: 2.226090
         Iterations 9
                          Poisson Regression Results                          
==============================================================================
Dep. Variable:                numbil0   No. Observations:                  197
Model:                        Poisson   Df Residuals:                      193
Method:                           MLE   Df Model:                            3
Date:                Thu, 15 Jul 2021   Pseudo R-squ.:                  0.8574
Time:                        00:14:01   Log-Likelihood:                -438.54
converged:                       True   LL-Null:                       -3074.7
Covariance Type:                  HC0   LLR p-value:                     0.000
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const        -29.0495      2.578    -11.268      0.000     -34.103     -23.997
lngdppc        1.0839      0.138      7.834      0.000       0.813       1.355
lnpop          1.1714      0.097     12.024      0.000       0.980       1.362
gattwto08      0.0060      0.007      0.868      0.386      -0.008       0.019
==============================================================================

¡Éxito! El algoritmo pudo lograr la convergencia en 9 iteraciones.

Nuestro resultado indica que el PIB per cápita, la población y los años de membresía en el Acuerdo General sobre Aranceles Aduaneros y Comercio (GATT) están relacionados positivamente con la cantidad de multimillonarios que tiene un país, como se esperaba.

Estimemos también los modelos más completos del autor y los mostraremos en una sola tabla.

regs = [reg1, reg2, reg3]
reg_names = ['Model 1', 'Model 2', 'Model 3']
info_dict = {'Pseudo R-squared': lambda x: f"{x.prsquared:.2f}",
             'No. observations': lambda x: f"{int(x.nobs):d}"}
regressor_order = ['const',
                   'lngdppc',
                   'lnpop',
                   'gattwto08',
                   'lnmcap08',
                   'rintr',
                   'topint08',
                   'nrrents',
                   'roflaw']
results = []

for reg in regs:
    result = sm.Poisson(df[['numbil0']], df[reg],
                        missing='drop').fit(cov_type='HC0',
                                            maxiter=100, disp=0)
    results.append(result)

results_table = summary_col(results=results,
                            float_format='%0.3f',
                            stars=True,
                            model_names=reg_names,
                            info_dict=info_dict,
                            regressor_order=regressor_order)
results_table.add_title('Table 1 - Explaining the Number of Billionaires \
                        in 2008')
print(results_table)

Table 1 - Explaining the Number of Billionaires                         in 2008
=================================================
                  Model 1    Model 2    Model 3  
-------------------------------------------------
const            -29.050*** -19.444*** -20.858***
                 (2.578)    (4.820)    (4.255)   
lngdppc          1.084***   0.717***   0.737***  
                 (0.138)    (0.244)    (0.233)   
lnpop            1.171***   0.806***   0.929***  
                 (0.097)    (0.213)    (0.195)   
gattwto08        0.006      0.007      0.004     
                 (0.007)    (0.006)    (0.006)   
lnmcap08                    0.399**    0.286*    
                            (0.172)    (0.167)   
rintr                       -0.010     -0.009    
                            (0.010)    (0.010)   
topint08                    -0.051***  -0.058*** 
                            (0.011)    (0.012)   
nrrents                                -0.005    
                                       (0.010)   
roflaw                                 0.203     
                                       (0.372)   
Pseudo R-squared 0.86       0.90       0.90      
No. observations 197        131        131       
=================================================
Standard errors in parentheses.
* p<.1, ** p<.05, ***p<.01

El resultado sugiere que la frecuencia de multimillonarios se correlaciona positivamente con el PIB per cápita, el tamaño de la población, la capitalización del mercado de valores y se correlaciona negativamente con la tasa de impuesto sobre la renta marginal máxima.

Para analizar nuestros resultados por país, podemos trazar la diferencia entre los valores predichos y reales, luego ordenar de mayor a menor y trazar los primeros 15

data = ['const', 'lngdppc', 'lnpop', 'gattwto08', 'lnmcap08', 'rintr',
        'topint08', 'nrrents', 'roflaw', 'numbil0', 'country']
results_df = df[data].dropna()

# Use last model (model 3)
results_df['prediction'] = results[-1].predict()

# Calculate difference
results_df['difference'] = results_df['numbil0'] - results_df['prediction']

# Sort in descending order
results_df.sort_values('difference', ascending=False, inplace=True)

# Plot the first 15 data points
results_df[:15].plot('country', 'difference', kind='bar',
                    figsize=(12,8), legend=False)
plt.ylabel('Number of billionaires above predicted level')
plt.xlabel('Country')
plt.show()

_images / mle_33_0.png

Como podemos ver, Rusia tiene, con mucho, el mayor número de multimillonarios por encima de lo que predice el modelo (alrededor de 50 más de lo esperado).

Treisman usa este resultado empírico para discutir las posibles razones del exceso de multimillonarios de Rusia, incluido el origen de la riqueza en Rusia, el clima político y la historia de la privatización en los años posteriores a la URSS.

Fuente https://python.quantecon.org/mle.html