Monitorización de consumo energético en CC con el sensor INA226


En el pasado, muchos proyectos de amperímetros se basaban ​​en sensores de corriente de efecto Hall como el ACS712,  en amplificadores de detección de corriente de lado alto como el MAX4080SASA o fabricados con amplificadores operacionales , lo cual tiene un pequeño inconveniente:  estos sistemas tienen una salida analógica que luego debe digitalizarse para procesarse en un sistema de adquisición de datos como Arduino.

Para intentar superar esos inconvenientes existe un CI muy interesante como el INA226, chip fabricado por Texas Instruments, que es un sensor de 36V, 16-bit, ultra-precise i2c output current/voltage/power monitor w/alert.

 El sensor INA226 tiene una salida digital e incorpora un ADC de 16 bits por lo que se obtiene una gran exactitud y precisión. Mide la corriente y el voltaje y calcula la potencia mientras Arduino se comunica con el chip, pudiendo presentar las medidas en una pantalla LCD y almacenarlas en una tarjeta micro SD. 

Este chip opera con un voltaje máximo de 36 voltios, mientras que la corriente está limitada solo por la derivación utilizada. Hay algunas bibliotecas para el chip INA226, peros se puede utilizar la biblioteca Korneliusz Jarzebski que me parece bastante completa incluso si tuviera que hacer algunos cambios en dos funciones.

Existen numerosas aplicaciones posibles para esta herramienta de monitorización: dispositivos alimentados por batería como scooters o bicicletas asistidas por pedaleo, paneles fotovoltaicos, etc.

En las mediciones de corriente con el shunt hay dos formas de insertarlo:

  • Tierra (lado bajo): la derivación está conectada entre la carga y la tierra.
  • Hacia la fuente de alimentación (lado alto): el shunt está conectado entre la fuente de alimentación y la carga.

El circuito integrado INA226, de Texas Instruments, es un dispositivo digital que mide la corriente con una derivación de lado alto o lado bajo y también mide el voltaje, calcula la potencia y proporciona una alarma multifuncional. El diagrama de bloques de dicho integrado es el siguiente:

INA226 - SBOS547A

La resolución de la tensión de derivación es de 2,5 m V con una escala completa de 32768×2,5 = 81,92 mV. Para la tensión VBUS, la resolución es de 1,25 mV con un fondo de escala teórico de 40,96 V, incluso si no se deben superar los 36 V. La resolución de la potencia es 25 veces mayor que la de la corriente, con un fondo de escala que depende del shunt utilizado. De modo que el sistema tiene una precisión de medición notable.

El ADC interno se basa en un convertidor delta-sigma de 16 bits (ΔΣ) con una frecuencia de muestreo típica de 500 kHz (± 30%), por lo que también es adecuado para corrientes que cambian rápidamente con el tiempo.El chip tiene 10 pines y dimensiones muy pequeñas, con una caja DGS (VSSOP).

El INA226 puede proporcionar una alerta de hardware o software si una variable, seleccionada por el usuario, ha superado un límite. El usuario puede seleccionar una de las cinco funciones disponibles para monitorear y / o establecer el bit de conversión lista.  Las cinco funciones de alerta que se pueden monitorear son:

  • Sobre-límite de voltaje de derivación (SOL): excede el umbral de corriente máximo;
  • Bajo límite de voltaje de derivación (SUL): excede el umbral de corriente mínimo;
  •  Límite de voltaje de bus (BOL): excede el umbral de voltaje máximo;
  • Límite inferior de voltaje de bus (BUL): superando el umbral de voltaje mínimo;
  •  Power Over-Limit (POL): superando el umbral de potencia máxima;

La salida de alerta es de colector abierto y se puede conectar fácilmente a un dispositivo de bloqueo.  Para simplificar, podemos leer la alerta de tipo de sobre-límite de voltaje de derivación (SOL) a través del software y presentarla en la pantalla. Gracias a este sensor pues podemos calcular la potencia con un Arduino presentando estas en una pantalla LCD almacenándolas en una tarjeta micro SD. 

Este chip funciona con un voltaje máximo de 36 voltios, mientras que la corriente está limitada solo por la derivación utilizada. La resolución del voltaje de derivación es de 2.5 m V con una escala total de 32768×2.5 = 81.92mV. 

El dispositivo tiene dos pines de dirección, A0 y A1. La siguiente tabla enumera las conexiones de clavijas para cada una de las 16 direcciones posibles.

A1A0Dirección del esclavoA1A0Dirección del esclavo
GNDGND1000000SDAGND1001000
GNDVS1000001SDAVS1001001
GNDSDA1000010SDASDA1001010
GNDSCL1000011SDASCL1001011
VSGND1000100SCLGND1001100
VSVS1000101SCLVS1001101
VSSDA1000110SCLSDA1001110
VSSCL1000111SCLSCL1001111

El módulo monta dos resistencias pull-down (R2 y R3), por lo que la dirección es 0x40 si no conectamos los jumpers colocados en el lado opuesto al de los componentes (ver figura 2 a la derecha).Con la resistencia de derivación de 0.1 W (R100), montada en el módulo, hay una resolución de corriente de 2.5 m V / 0.1 = 0.025 mA, una escala completa de 81.92mV / 0.1 W = 819.2 mA y una resolución de potencia de 0.025 mA * 25 = 0,625 mW. 

Una opción muy aconsejable es montar un shunt externo para corrientes altas, desoldando el interno e insertar un filtro RC,

Pantalla LCD

LCDfunciónArduinoLCDfunciónArduino 
14Línea de bus de datos D7A35 5R / W – Lectura / EscrituraGnd
13Línea de bus de datos D6A24 4RS – Registrarse SeleccionarD7
12Línea de bus de datos D5D43VEE
11Línea de bus de datos D4D52VCC ( + 5V)VCC
6 6E – Habilitar señalD61VSS (GND)Gnd

El Ci se puede complementar con una pantalla LCD común de dos líneas de 16 caracteres con controlador compatible con Hitachi HD44780 y una retroiluminada equipada con diodos LED de alta eficiencia que consume alrededor de 20 mA, que puede reducierse a 10 mA con una resistencia externa en serie con la iluminaciona. Las conexiones con Arduino Nano son las siguientes:

El módulo micro SD

La siguiente tabla muestra las conexiones y pines dedicados a la microSD:

J2Módulo SDArduino
6 6CSD10
5 5SCKD13
4 4MOSID11
3MISOD12
2VCC+5 V
1GNDGND

Cálculos para la derivación utilizada

Podemos usar cualquier SHUNT de potencia. Como ejemplo el ICE 680R. Este shunt tiene una salida de 100 mV con una corriente de 25 A, por lo que Rs = 0.1 / 25 = 0.004 Ω.
Desoldand la derivación de 0.1 Ω, montada en el módulo INA226, y soldarmos un conensador de 1 MF y dos resistencias de filtro. Este recurso permite el filtrado RC, como sugiere el fabricante del chip.

Con el voltaje de derivación a escala completa en mente, la corriente máxima medible ahora se convierte en:Ifs = 81,92 mV / 4 mΩ = 20,48 A

La resolución actual será: 20,48 / 32768 = 0,625 mA.

La función de calibración tiene como parámetros de entrada la resistencia en derivación y la corriente máxima, en nuestro caso: rShunt = 0.004 Ω e iMaxExpected = 20.48 A.

Calculo de varias variables:

  • currentLSB = iMaxExpected / 32768 = 20,48 / 32768 = 0,625 [mA]
  • alibrationValue = 0,00512 / currentLSB / rShunt = 0,00512 / 0,000625 / 0,004 = 2,048
  •    powerLSB = currentLSB * 25 = 0.000625 * 25 = 15.625 [mW]

Otro dato para pasar al programa es el umbral actual de la función setShuntVoltageLimit .Si establecemos el límite superior del instrumento en 20 A, podemos calcular:setShuntVoltageLimit = 0,004 Ω * 20 A = 0,08 V

El montaje del circuito final

Podemos usar un Arduino Nano, pero por supuesto se puedne usar otras placas Arduino, como Arduino Uno o Arduino Pro, (pero esta placa es muy compacta, completa con un adaptador USB y se puede montar en un prototipo de pcb pretaladrado).

Arduino Nano tiene su propio regulador de 5V, un Low DropOut tipo AMS1117, pero podemos alimentar el sistema con un LM7805 común, esto se debe a que el primero acepta un voltaje de entrada máximo de 15V contra los 35V del segundo. Si tenemos que alimentar ArduINA226 con el voltaje que tengo que monitorear, es mejor tener valores compatibles. Ademas otra razón es la luz de fondo de la pantalla. El consumo de todo el sistema ronda los 45 mA. El diodo Schottky D1, que tiene una caída de solo 0,2 V, se utiliza para evitar conflictos entre la fuente de alimentación externa y la fuente de alimentación USB.

La figura siguiente muestra el diagrama de cableado :

Lista de componentes

componentedescripcióncomponentedescripción
C1, C2, C4100 nF 50V, condensador cerámicoR4Resistencia de 100 ohmios , ± 5%, 0.5 W
C310 MF 50V, condensador electrolíticoU1Tablero Arduino Nano
C41 MF 25V, condensador cerámicoU2LM7805, regulador de 5V
D11N5819, diodo SchottkyU3Módulo INA226
Rp1Recortadora resistiva de 22k WmonitorLCD de 2 filas x 16 columnas
R1, R210 ohmios , ± 5%, resistencia de 0.25 Wmódulo SDTarjeta micro SD con niveles de 5V
R34.7 kohmios , ± 5%, resistencia de 0.25 WSw1Pulsador normalmente abierto (NO)

Nota : el prefijo de capacidad ‘M’ significa microfaradio (1e-6 F)

 

Cambios en la biblioteca INA226.cpp

Hay algunas bibliotecas para el chip INA226, siendo la biblioteca Korneliusz Jarzebski bastante completa cmbiando la función de calibración porque la corriente de derivación era incorrecta. La función setShuntVoltageLimit también estaba mal.

 

Aquí están las funciones modificadas, para reemplazar las originales:

bool INA226::calibrate(float rShunt, float iMaxExpected){ // MODIFIED by GCAR
    uint16_t calibrationValue;
    float iMaxPossible;
    iMaxPossible = vShuntMax / rShunt;
    currentLSB = iMaxExpected / 32768;// calculate current resolution
    powerLSB = currentLSB * 25;// power resolution
    calibrationValue = (uint16_t)((0.00512) / (currentLSB * rShunt));
    writeRegister16(INA226_REG_CALIBRATION, calibrationValue);
    return true;
}

void INA226::setShuntVoltageLimit(float voltage){// MODIFIED by GCAR
    uint16_t value = voltage/2.5e-6;
    writeRegister16(INA226_REG_ALERTLIMIT, value);
}


El programa

El programa requiere solo dos parámetros que son: rShunt e iMaxExpected . Si desea utilizar una alerta, debe decidir qué señal le interesa. Si elige la corriente máxima, le interesara el voltaje de derivación máximo y el uso las funciones enableShuntOverLimitAlert y setShuntVoltageLimit . 

En el caso de un sistema alimentado por baterías de iones de litio o plomo, sería más útil comprobar la tensión mínima de VBUS con las funciones enableBusUnderLimitAlert y setBusVoltageLimit .

Si no hemos insertado la tarjeta SD o si no es válida, aparecerá el mensaje “¡SD no presente!” 

.Al final de la función setup (), el programa imprime “Push to Start” y espera a que se presione el botón SS, si el SD está presente, comienza a adquirir las medidas en archivo, de lo contrario las imprime solo en la pantalla .

En cuanto a la pantalla LCD, dos filas y 16 columnas son suficientes, por lo que debe administrar bien sus impresiones. Teniendo en cuenta la resolución de las variables a imprimir, estos son sus formatos:

  •  Voltaje del bus: V = xx.xxx (8 caracteres)
  • Corriente de derivación: I = xx.xxx (8 caracteres)
  • Potencia del bus: W = xxx.xx (8 caracteres impresos)

Entonces, en la primera línea veremosel voltaje y la potencia, separados por un espacio y en la segunda la corriente y el número de muestras guardadas en SD (si está insertada), de esta manera:V = xx.xxx W = xxx.xxI = xx.xxx N = xxxxx

Si se produce la alerta cuando se excede el límite máximo de corriente (SOL, sobre límite de voltaje de derivación), aparecerá la impresión “ALERTA SOL” en la segunda línea, en lugar del valor actual y en la primera el voltaje de derivación:Derivación V = xx.xxxxxALERTA SOL

/ * Programar ArduINA226 a corriente, voltaje y potencia
 con módulo Arduino Nano e INA226
utiliza la biblioteca Korneliusz Jarzebski para INA226 (modificado)
guardar las mediciones en SD si está presente
Giovanni Carrera, 14/03/2020
 * /

#include <SPI.h>
#include <LiquidCrystal.h>
#include <Wire.h>
#include <SD.h>
#include <INA226.h>
// Pines de LCD
#define rs 7
#define en 6
#define d4 5
#define d5 4
#define d6 A2
#define d7 A3
#define SSbutton 2
#define SD_CS 10
// inicializar la biblioteca asociando cualquier pin de interfaz LCD necesario
LiquidCrystal lcd (rs, en, d4, d5, d6, d7);
INA226 ina;


char bline [17] = ""; // línea en blanco
const int deltat = 500; // período de muestreo en milisegundos
cmilli largo sin firmar, pmilli;
boolean SDOk = verdadero, FileHeader = verdadero, ACQ = falso;
unsigned int ns = 0;

void setup () {
  // configura el número de columnas y filas de la LCD:
  lcd. comienzo (16, 2);
  pinMode (SSbutton, INPUT); // Botón de inicio / parada
  pinMode (SD_CS, SALIDA); // Selección de chip SD
  // La dirección INA226 predeterminada es 0x40
  ina.begin (0x40);
  lcd.print ("ArduINA226");
  lcd.setCursor (0, 1); // imprime en la segunda fila
  lcd.print ("Monitor de energía");
  // Configurar INA226
  ina.configure (INA226_AVERAGES_1, INA226_BUS_CONV_TIME_1100US, INA226_SHUNT_CONV_TIME_1100US, INA226_MODE_SHUNT_BUS_CONT);
  // Calibrar INA226. Rshunt = 0,004 ohmios, corriente máxima esperada = 20,48 A
  ina.calibrate (0,004, 20,48);
  ina.enableShuntOverLimitAlert (); // habilita la alerta de sobretensión de derivación, corriente por encima del límite
  ina.setShuntVoltageLimit (0.08); // límite de corriente = 20 A para 0,004 ohmios
  ina.setAlertLatch (verdadero);
  if (! SD.begin (SD_CS)) {
    LCDprintLine ("¡SD no está presente!", 1);
    delay (5000);
    SDOk = false;
  }
  LCDprintLine ("Pulsar para iniciar", 1);
  while (digitalRead (SSbutton) == HIGH) {}; // esperar el inicio
  if (SDOk) ACQ = true;
}


void loop () {
  cmilli = milis ();
  si (cmilli - pmilli> deltat) {
    pmilli = cmilli;
    voltios flotantes = ina.readBusVoltage ();
    corriente flotante = ina.readShuntCurrent ();
    LCDprintLine ("V =", 0); // imprimir voltaje del bus
    lcd.print (voltios, 3); 
    float power = ina.readBusPower (); // INA calcula la potencia
// potencia flotante = voltios * corriente; // Arduino calcula la potencia
    lcd.print ("W ="); // potencia de impresión
    lcd.print (potencia, 4);
    flotar Vshunt = ina.readShuntVoltage ();
    String dataString = String (voltios, 3) + ',' + String (actual, 4) + ',' + String (potencia, 4) + ',' + String (Vshunt, 6);
    si (ACQ) {
      Archivo dataFile = SD.open ("powerlog.csv", FILE_WRITE);     
      if (dataFile) {
        if (FileHeader) {// imprimir el encabezado del archivo
          dataFile.print ("Deltat [ms] =");
          dataFile.println (deltat);
          dataFile.println ("Vbus [V], Ishu [A], P [W], Vshu [V]");
          FileHeader = false;
        }
        dataFile.println (dataString);
        ns ++; // número de muestras adquiridas
        dataFile.close ();
        if (digitalRead (SSbutton) == LOW && ns> = 10) {// detente después de al menos 10 muestras
          ACQ = false; // detener la adquisición
          LCDprintLine (Cadena (ns), 1);
          lcd.print ("muestras");
          delay (5000);        
        }
      } else {
        LCDprintLine ("¡No se puede abrir el archivo!", 1);
        ACQ = false;
        delay (5000);
      }
    }
    if (ina.isAlert ()) {
      LCDprintLine ("Shunt V =", 0);
      lcd.print (Vshunt, 5);
      LCDprintLine ("ALERTA SOL", 1);
    }
    else {
      LCDprintLine ("I =", 1);
      lcd.print (actual, 3);
      if(ACQ) {
        dataString = "N =" + String (ns);
        lcd.print (dataString); // número de impresión de la muestra adquirida
      }     
    }
  }
}
/ ************************** Funciones ********************* **** /
void LCDprintLine (texto de cadena, línea de bytes) {
   lcd.setCursor (0, línea);
   lcd.print (bline); // borra la segunda fila
   lcd.setCursor (0, línea);
   lcd.print (texto); // imprimir texto
}

El sistema puede funcionar incluso con periodos delta muy pequeños, pero la visualización, las impresiones y la escritura en archivo llevan algo de tiempo: si queremos reducirlo, por ejemplo a 100 ms, las medidas deben mostrarse cada 5 muestras. 

También se pueden implementar algunos circuitos de hardware, como usar la salida de alerta para desconectar la fuente de alimentación a la carga a través de un relé estático o electromecánico o para hacer un sonido de timbre a través de Arduino.

ArduINa226, si se ha insertado una tarjeta micro SD, transfiere las medidas al archivo “powerlog.csv”.La adquisición comienza presionando el botón ‘SS’ (Start / Stop) que también sirve para finalizar la grabación después de al menos diez muestras.

Las dos primeras líneas contienen el período de muestreo deltat en milisegundos y los nombres de las señales, como se ve en el siguiente ejemplo:Deltat [ms] = 500Vbus [V], Ishu [A], P [W], Vshu [V]10.021,0.0950,0.9531,0.00037510.023,0.0944,0.9531,0.00037710.023,0.0944,0.9531,0.000375……Los datos se agregan al archivo que no se sobrescribe.

Calculo de la potencia

La energía suministrada por la fuente es la integral de la potencia a lo largo del tiempo. El sistema más simple para hacer la integral definida de una función, aunque de manera bastante aproximada, es el de los trapezoides. La curva de potencia se aproxima en muchos trapecios, de los cuales es fácil calcular el área, para dos puntos tenemos:

La integral será igual a la suma de las áreas individuales. 

Este método dará mejores resultados para intervalos de tiempo más pequeños, en nuestro caso:

(t2-t1) = deltat yf (t) es la potencia P (t), entonces

E2 = deltat * (P1 + P2) / 2 + E1 [J]

Donde P1 y P2 son las potencias en dos instantes sucesivos t1, t2 y E1, E2 las energías respectivas expresadas en W × s = julios , que se pueden expresar fácilmente en vatios hora [Wh] dividiendo los julios por 3600

Para calcular la energía con Arduino:

watthour = watthour + (potencia + ppotencia) /2.0*deltath; // entero para trapezoides

ppower = potencia; // luego actualiza el poderDonde potencia es la potencia actual , ppower la del instante anterior ydeltath = deltat / 3.6e6; // período de muestreo en horas (deltat en [ms])

Al principio, debe restablecer la variable watthour  pero el autor no agregoesta función porque no tenía más espacio en la pantalla y preferio calcular la energía en el PC con algoritmos de integración más sofisticados.

Fuente :http://ardupiclab.blogspot.com/2020/03/the-arduina226-power-monitor.html

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