text:
Servicios
Buscar
Idiomas
P. Completa
AN ESTIMATION OF COSTS AND PUBLIC-HEALTH BENEFITS BY THE PM10 MITIGATION IN MEXICALI, BAJA CALIFORNIA, MEXICO
Marco Antonio Reyna Carranza; Efraín Carlos Nieblas Ortíz; Martha Lorena Nava Martínez;
Marco Antonio Reyna Carranza; Efraín Carlos Nieblas Ortíz; Martha Lorena Nava Martínez; Esmeralda Torillo Portilla
AN ESTIMATION OF COSTS AND PUBLIC-HEALTH BENEFITS BY THE PM10 MITIGATION IN MEXICALI, BAJA CALIFORNIA, MEXICO
UNA ESTIMACIÓN DE COSTOS Y BENEFICIOS EN SALUD PÚBLICA POR LA DISMINUCIÓN DEL PM10 EN MEXICALI, BAJA CALIFORNIA, MÉXICO
Revista Internacional de Contaminación Ambiental, vol. 33, no. 1, 2017
resúmenes
secciones
referencias
imágenes

Abstract: This research provides a first approach to the assessment of the health benefits and health costs that would be generated in Mexicali if there were or were not control measures implemented, which could create an impact in the decrease of concentrations of particulate matter ≤ 10 micrometers in aerodynamic diameter (PM10) in the air of the city. Several different scenarios were established to study the projected impact for the 2013─2020 periods. The public health benefits (PHB) are calculated based on the exposure-response function (ERF) estimated using local data, as well as ERF obtained from scientific literature applied previously in local mexican studies. Using a discount rate of 3 %, total costs and total social benefits are brought to present value and are shown as a percentage of the state’s gross domestic product (GDP) reported in 2011. In Mexicali, the study estimates that if there were no measures taken to reduce PM10 pollution in the studied years, there would be a social cost of around $1659 (lower value:$1164, higher value: $2503) millions of dollars for the period (2013─2020). This represents 5.59 % (3.92 %, 8.43 %) of the GDP. If there were measures taken to reduce the PM10 concentrations by about 8 % a year, in order to comply with the federal standards in 2020, there would be savings of around$633 ($444,$955) millions of dollars for the studied period. This would represent 2.13 % (1.50 %, 3.22 %) of the GDP. This study justifies the implementation of control measures for air pollution in Mexicali.

Keywords: air-pollution, morbidity, mortality, economic-impacts.

Resumen: El presente estudio brinda un primer acercamiento a la valoración de los beneficios en salud y el costo social que se generarían en Mexicali si se establecieran o no medidas de control que impactaran en la reducción de las concentraciones del material particulado ≤ 10 micrómetros de diámetro aerodinámico (PM10). Se establecieron distintos escenarios para el estudio de los impactos proyectados para el periodo del 2013 al 2020. Los beneficios en salud pública (BSP) se calcularon tanto con funciones exposición-respuesta (FER) estimadas con datos de la localidad, como con FER obtenidas de la literatura científica que habían sido aplicadas con anterioridad en estudios de localidades mexicanas. Los costos y/o beneficios sociales totales se ubican en un valor presente con tasa de descuento del 3 % y se presentan como porcentajes del producto interno bruto (PIB) estatal reportado en 2011. En Mexicali las proyecciones estiman que de no realizarse acción alguna para disminuir la contaminación por PM10 se produciría un costo social de alrededor de US $1659 (valor bajo:$1164, valor alto: $2503) millones, esto representaría el 5.59 % (3.92 %, 8.43 %) del PIB. De llevarse a cabo medidas que disminuyeran las concentraciones en aproximadamente un 8 % anual para cumplir con la norma en 2020, se tendría un ahorro de alrededor de US$633 ($444,$955) millones, lo que representaría el 2.13 % (1.50 %, 3.22 %) del PIB. El estudio justifica la implementación de medidas de control de la contaminación del aire en Mexicali.

Carátula del artículo

AN ESTIMATION OF COSTS AND PUBLIC-HEALTH BENEFITS BY THE PM10 MITIGATION IN MEXICALI, BAJA CALIFORNIA, MEXICO

UNA ESTIMACIÓN DE COSTOS Y BENEFICIOS EN SALUD PÚBLICA POR LA DISMINUCIÓN DEL PM10 EN MEXICALI, BAJA CALIFORNIA, MÉXICO

Marco Antonio Reyna Carranza
Universidad Autónoma de Baja California, México
Efraín Carlos Nieblas Ortíz
Universidad Autónoma de Baja California, México
Martha Lorena Nava Martínez
Instituto Mexicano del Seguro Social, México
Esmeralda Torillo Portilla
Secretaría de Finanzas del Distrito Federal, México
Revista Internacional de Contaminación Ambiental, vol. 33, no. 1, 2017

Accepted: 01 June 2016

Funding
Funding source: Environmental Protection Agency (EPA) and the Border Environment Cooperation Commission (BECC)
Award recipient: Los autores
INTRODUCTION

The proposal is to study the scope of impact that high levels of air pollution have in the economy and competitiveness of the population in Mexicali, derived from health impact. In order to achieve this, the costs attributed to public health benefits (PHB) are estimated by establishing scenarios in the reduction of particulate matter equal or smaller than 10 micrometers in aerodynamic diameter (PM10), one of the main air pollutants in the city of Mexicali. In fact, the World Health Organization (WHO) has recently indicated that Mexicali is considered the 4th most polluted city in the world (WHO 2011, MEBC 2012) due to particulate matter equal or smaller than 2.5 micrometers in aerodynamic diameter (PM2.5), despite the fact that the National Institute of Ecology and Climate Change (INECC for its acronym in spanish) had previously warned about the serious pollution problem that this city shows (Zuk 2007).

In order to estimate the PHB, the method employed is “health impact assessment (HIA)”, using as a guide the proposal presented by the INECC (2011).

To estimate the economic impact of air pollution, the potential benefits (impacts that could be avoided) that could be achieved in health (i.e. PHB) were the first to measure, due to the action of reducing the level concentration of the air pollutant considered in this study (i.e. PM10). Then, the estimated PHB was economically assessed by the human capital approach for mortality and the cost of illness (CoI) approach for morbidity.

Several studies show that PM10 can be used as a causal pollutant (Schwartz 1992, Ostro 1993, Schwartz 1993, WHO 1996), besides the fact of being a good indicator for the presence of other air pollutants (i.e. a combination of particles from different sources). Therefore, the PHB estimated by PM10 will be understood as being a consequence of the causal relationship between the total mix of this pollutant and their impacts in local health.

Even though the WHO states that all the effects on health as a consequence of the exposure to air pollution are important, it is not practical to include them all in a HIA (WHO 2000). The health impacts considered in this study were selected based on the local data available and on the recommendations made in a forum of experts on binational surveillance of diseases related to air pollution, coordinated by the Universidad Autónoma de Baja California-Mexicali in May 18th, 2006 (SDSU-UABC 2006). The best indicators that would help to measure the effects of air pollution on the local human population were defined. Therefore, the health impacts assessed in this document are: 1) mortality, 2) emergency room visits due to asthma (J45 in ICD-10 codes) and asthmatic status (J46 in ICD-10 codes), 3) hospitalizations due to asthma and asthmatic status, and 4) restricted activity days (RAD).

The aim of this study was to assess, as a first approach, the benefits in terms of Public Health (i.e. PHB for avoided cases) attributable to different PM10 reduction scenarios in Mexicali, Baja California, Mexico.

MATERIALS AND METHODS

To achieve the estimation of the avoided health impacts during the 2013 to 2020 period, data projected by the authors is used from these years. The projections for each variable were made based on their trend, using linear regression models for population projections and linear or nonlinear regression models for other variables which best suit (R2) the historical data for years previous to 2013.

Population data

The total human population data previous to 2013 were taken from the National Population Council (CONAPO for its acronym in spanish of Consejo Nacional de Población; 2010). The economically active population data were taken from the National Institute of Statistics and Geography (INEGI for its acronym in spanish of Instituto Nacional de Estadística y Geografía, years 1990, 1994, 1995, 2000, 2005, 2010), and the Mexicali emissions inventory (ICAR 1999).

Health data
Mortality

The human mortality data previous to 2013 were provided by the Ministry of Health from the Government of the State of Baja California (SA for its acronym in spanish of Secretaría de Salud, pers. comm.). The rates per every 100 000 inhabitants are estimated based on the total mortality count and number of inhabitants per year.

Morbidity
Visits to the emergency room due to asthma and status asthmaticus

The data of new cases of asthma and new cases of status asthmaticus were provided by the SA from the Government of the State of Baja California (pers. comm.). We begin from the assumption that every new case implies a visit to the emergency room. Therefore in the future, the new cases will be referred to as visits to the emergency room due to asthma. The counts per year for the 2013 to 2020 period were estimated by the authors based on the trends from years previous to 2013. It is assumed that asthma hospitalizations also count as visits to the emergency room, so in order to prevent double counts at the time of estimating the avoided cases, net emergency room visits are obtained by subtracting the hospitalizations from emergency room. Then, the net rates of visits to the emergency room per every 100 000 inhabitants are calculated for each year projected.

Hospitalizations due to asthma

Since there were no reliable sources found showing the incidence rates of hospitalization due to asthma, these rates were obtained using data from the Imperial County in California, USA, from the California Department of Public Health (CDPH 2015), California Inpatient Discharge Data (OSHPD 2001, 2004) and the program of the Lucile Packard Foundation for Children’s Health (Kidsdata 2016). It is assumed that the number of hospitalizations related to the emergency room visits due to this pathology on the mexican side, has the same proportion as that on the American side since the Imperial County and the Mexicali Valley are part of the same air basin. The ratio of the number-of-hospitalizations / emergency-room-visits due to asthma that occurs in the Imperial Valley is calculated (i.e. 17 %). This ratio is the average of the ratios occurred from 2005 to 2009. It starts from the assumption that the estimated average ratio will not vary significantly throughout the projected period of years (2013-2020). Then the ratio are multiplied by the number of visits to the emergency room due to asthma occurring in the municipality for each one of the years projected. Finally the hospitalization rate for every 100 000 inhabitants is obtained.

Since there was no sources found showing the RAD rate for Mexicali, this was deducted from the reports in scientific literature (Ostro 1987, WHO 1996). Based on national surveys conducted in the USA it is indicated that the average RAD occurring per person per year is 19, including: work days lost (WDL), days spent in bed, and minor restrictions (i.e. activities are partially restricted due to illness). The number of RAD was taken from the adult population between the ages of 18 and 65. In Mexico, economically active population is considered to be between the ages of 15 and 64. So, for purposes of this study, the average RAD is multiplied by the busy economically active population per year for the projected period (2013─2020) to obtain the total number of RAD in the busy population. From these amounts, the number of hospitalization days and the net number of emergency room visits were subtracted, both, for the occupied population due to the fact that these events also generate restricted activity. Thus, the net RAD per year is obtained.

Pollution data

PM10 data from 1998 to 2002 and from 2009 to 2012 were provided by the Mexican Ministry of Environmental Protection (SPA for its acronym in Spanish of Secretaría de Protección al Ambiente) from the manual environmental monitoring stations (pers. comm.), and data from 2003 to 2008 were downloaded by the authors from the INECC web page (2016). Annual average of 24-h averages were estimated in order to project the annual average values for the 2013 to 2020 period in terms of the trend for years previous to 2013. In order to be consistent with the exposure levels used in the studies reported in the literature which the ERF were extracted, the current HIA uses average exposure concentrations and no other statistics as might be maximum values or percentiles.

Figure 1 displays the behavior of the annual PM10 averages in Mexicali from 1998 to 2012. It shows that the pollution levels maintain a positive trend.

Fig. 1
Annual historical averages of PM10 in Mexicali. Data provided by the Secretaría de Protección al Ambiente (SPA) of the Government of the State of Baja California. a According to the SPA, the value can be sub-estimated due to the fact that the indicator was obtained with information of less than 75 % of the valid data. b Values estimated by the authors with data from the manual stations provided by the SPA. (Source: Pers. comm).

Figure 2 shows the annual averages of PM10 occurred during the period from 2000 to 2012 and those projected until 2020 for Mexicali.

Fig. 2
Annual historical and projected PM10 averages for Mexicali. Data provided by the Secretaría de Protección al Ambiente (SPA) of the Government of the State of Baja California until the year 2012. From year 2013 are projections. For the purpose of such projections, the original values for 2005 and 2006 (i.e. values 210a and 254 of Fig. 1 ) were eliminated since they were extreme severe data and replaced by values estimated by the authors using regression models. Five regression models (i.e. exponential, linear, logarithmic, polynomial and power) were tested and was used that with the highest R2. a According to the SPA of the State Government, the value can be underestimated due to the fact that the indicator was obtained using information that was less than 75 % of the valid data.1 Data estimated by the authors using regression models. 2 Data estimated by the authors using data from the manual stations (COBACH, CONALEP, ITM, PROGRESO and CESPM) provided by the SPA. (Source: pers. comm.).

As shown in figure 2 , PM10 in Mexicali has been showing a systematic increase year by year. If it remains so, it could reach average levels of around 141 µg/m3 of air (95 % IC: 134, 148) in the year 2020 according to the projected trend.

Scenarios

The following are the three scenarios established in the decrease of PM10 pollution levels for Mexicali:

(a) A total decrease per year starting from 2013 and until 2020.

(b) A reduction to the levels indicated by the NOM-025-SSA1-1993 (SSA 1993) of 50 mg/m3 from 2013 to 2020.

(c) A systematic reduction of about 8 % each year from 2013 to reach in 2020 comparable levels to those indicated by the NOM- 025- SSA1-1993 of 50 µg/m3.

Avoided cases

In this present study, avoided cases of mortality and morbidity are calculated combining local ERF with others that are likewise obtained from scientific publications, but the condition being that these should had been previously applied to communities in Mexico.

The number of avoided cases due to changes in pollutant concentration is estimated using the equation (1) (USEPA 2003, 2012), with the assumption that the effects of the pollutant on health are continuous, from the lowest to the highest levels. In other words, the assumption is that there are no threshold effects in the ERF.

(1)

In which:

Hij = Number of cases of impact on health i related to the concentration of the pollutant j.

R i = Basal rate of mortality or morbidity for population N by effect i (cases / inhabitants / year).

Cj = Change of the pollutant concentration j.

N = Exposed population (number of people).

βij = Exposure-response function by effect i due to the exposure to the pollutant j (increase of cases / unit of concentration).

Percentage changes and variability

To visualize variability both in the number of avoided cases as in associated costs for each scenario proposed, their respective low, medium and high values are estimated. These follow similar criteria to that were applied in one of the studies conducted by the INECC (McKinley et al. 2003), because aside from having been done in the country, it addresses aspects that are very similar to those covered in the local study in as much as calculation of the effects of PM10 on the avoided cases of mortality and morbidity.

Mortality

To calculate the low, central and high values the aforementioned INECC study (McKinley et al. 2003) proposes using three ERF that indicate the percentage changes per each 10 μg/m3 of air, which are: 0.5 % for the low value used by Samet et al. (2000), 0.7 % for the central value used by Levy et al. (2000), and 1.4 % for the high value proposed by Evans et al. (2000). For this study, a 0.35 % percentage change is used to calculate the low level. This was taken from a time-series study performed with local data from Mexicali that can be found in scientific publications in Reyna et al. (2012). The 0.5 % from Samet et al. (2000) is used to calculate the central value, and the 0.7 % from Levy et al. (2000) is used to calculate the high value.

Emergency room visits due to asthma

No ERF estimated with local data was found to calculate the avoided cases of morbidity, thus the authors generated them using data from the region. The estimation of local ERF follows a procedure that is similar to that established in the scientific literature (Yang et al. 2009), which is based on auto-regressive semi-parametric models using as a measuring unit weekly counts of the cases, instead of daily counts. This method is well suited as to the limitations of local morbidity data, which represent total counts per week and not total daily counts. The time-series used are from 2005 to 2007. The method used in the estimation of local ERF, implements auto-regressive generalized additive models (GAM) of the Poisson family, where it is first found a baseline model in which the dependent variable are weekly counts of asthma emergency room visits, and as confounding variables are the trend, seasonality, temperature, relative humidity, holidays, vacation days and lagged values of the dependent variable depending on the case of the lag or lags where statistically significant values occur in residual autocorrelation coefficients. The trend, seasonality, temperature and relative humidity are smoothed with cubic spline functions, whereas holidays and vacations are introduced into the model as dummy variables. The degrees of freedom of each spline and the adequacy of the baseline model are determined once the coefficients of the partial autocorrelation of the residuals do not exceed confidence limits (α = 0.05). Afterwards, the pollutant and its lags are added to the baseline model one by one. The statistical significance of the associated regression coefficients is measured and the one showing more significance is selected. This coefficient is used to estimate the exposure-response function of the impact that the PM10 pollutant has on asthma.

The low, central and high levels of emergency room visits for asthma are calculated using weighed averages (WHO 1996): the percentage change mentioned by the INECC (McKinley et al. 2003) of 4 % (95 % CI: 1 %, 7 %) proposed by Schwartz et al. (1993) and the percentage change estimated locally in Mexicali of 2.5 % (95 % CI: 0.43 %, 4.7 %). Each percentage change and their confidence intervals (CI) are weighted with the inverse of their respective variances. Then the average is calculated to obtain the weighted percentage change for every 10 μg/m3 of PM10: 3.05 % (95 % CI: 0.64 %, 5.5 %). And thus, the low level is estimated with the left average confidence interval, the central level with averaged percentage change, and the high level with the right average confidence interval.

Hospitalizations due to asthma

The 2.5 % (95 % CI: 0.43 %, 4.7 %), local percentage change utilized for the asthma hospitalization is the same as that calculated for asthma emergency room visits since it is assumed that the number of hospitalizations in this location depends on the number of emergency room visits for this pathology, which in Mexicali represents 17 % (see “Health Data” section). On the other hand, the INECC study (McKinley et al. 2003) utilizes a 3.02 % (95 % CI: 2.05 %, 4.00 %) percentage increase, stating such is proposed by the World Bank.

For the purpose of the current study these two percentage increments are weighted in the same way that those to emergency room visits as explained above to obtain the average. This new 2.94 % (95 % CI: 1.79 %, 4.12 %) percentage increment is used to calculate the low, central and high levels of hospitalization due to asthma.

Since no local dose-response functions were found to calculate the low, central and high levels of RAD; the decision was made to follow criteria established by WHO (1996). This proposes a percentage increment of 3.0 % (95 % CI: 2.1 %, 4.7 %) per 10 μg/m3 of air of PM10 which is an adaptation of the original percentage increment associate to the effects of PM2.5. The adaptation is made using the relation . For our study the percentage increment is adapted to the relation that was calculated with the actual PM10 and PM2.5 levels in Mexicali. Therefore, the resulting incremental percentage utilized is 3.5 % (95 % CI: 2.4 %, 5.5 %) per 10 μg/m3 of PM10.

Economic assessment
Mortality

In this study, the economic assessment of PHB for mortality is estimated by the human capital approach, with the present value of future income that an individual would perceive if he had not died prematurely, by using the following equation (Sanchez et al. 1998, USEPA 2012):

(2)

Where: PVj is the present value of future income for an individual of age j, Prob[n / j] is the probability that an individual of age j is alive at age n, wn is the average annual labor income by individual of age n, and r is the discount rate for the years after death.

In our study, the survival probabilities (i.e. npx = 1- nqx ) were constructed using the new official mortality charts (i.e. nqx ) (AMA 1999), under the assumption that these probabilities have not changed with time since they were constructed.

The average annual labor income wn is obtained from the general average wages and not from the minimum wage, since it is based on the fact that most of the employed population of Baja California perceives income equal to or above the general average wage. General average wages for Baja California recorded by the Mexican Institute of Social Security (IMSS for its acronym in Spanish of Instituto Mexicano del Seguro Social) in 2011 reached 245.84 mexican pesos per day, this represents 4.11 minimum wage in the geographical area A of the country (source: National Commission on Minimum Wages CONASEMI for its acronym in Spanish of Comisión Nacional de los Salarios Mínimos, and INEGI). So the minimum wage is multiplied by this factor to calculate the general average wage corresponding to each year, assuming that the proportion of the number of minimum wages relative to average wages is maintained approximately the same year by year. Then, the average wages are converted to the corresponding US$parity of each year, according to the historical series obtained from the Bank of Mexico. The minimum wages and parity of the dollar for the years 2013 to 2020 are projected following the tendency of previous years to 2013. It is assumed that the employed population in Mexicali perceives average wages with a similar percentage distribution of that assigned by INEGI in the perception of minimum wages. This scenario means that different proportions of the employed population receive different wages. With the intention to balance this situation, a weighted average wage is estimated in function of the respective employed populations. It is assumed that the employed population that earns general average daily wages keeps a similar distribution to the employed population that earns minimum wages. Then the estimated weighted average daily wages of each year (US$85 in the period of the projected years) are multiplied by 313 officially working days per year in Mexico.

When calculating the present values of future income (PV) for the loss of productivity due to premature death, a discount rate r of 3 % (see ecn 2) is applied to be consistent with other studies reported in the literature (Houtven and Cropper 1996, Sánchez et al. 1998).

Future income calculated by lost productivity are multiplied by the number of avoided death cases respective to every year, to finally obtain the costs attributable to avoided death cases for reducing the levels of PM10, under the proposed scenarios.

Morbidity

To calculate the assessment of morbidity effects (i.e. asthma hospitalizations, emergencies and restricted activity days), in this study, we used the costs associated to disease (i.e. direct cost of treatments and lost salaries due to days not worked), although there are other expenses: preventive or defensive expenses, (i.e. people living in a contaminated area are taking defensive measures to reduce the risk of sickness), also known as health production functions; and the contingent assessment, which is complicated and also expensive because of the implications of the number of gathered information and application of surveys.

As an indirect expense, restricted activity days are assessed, produced by the lost work days, for consequence of sickness, or for taking care of a sick person, or decrease of performance (Ostro 1987).

Costs associated to cases of emergency room visits due to asthma

For the purposes of this study, the associated costs to visits to the emergency room due to asthma are assessed based on the generated cost of the consultation, plus emergency attention, plus the expense of one work day lost. The produced expense of the medical consultation (i.e. US $157.3) is obtained from a study developed in the northeast of the country, by the University Hospital of the Regional Center of Allergy and Immunology of Nuevo Leon (Gallardo et al. 2007). The emergency attention expense is calculated with the estimated average (i.e. US$164.06) in the northeast study (i.e. US $97.33), with the estimated in a study made by the National Institute of Respiratory Diseases (INER for its acronym in Spanish of Instituto Nacional de Enfermedades Respiratorias) (Tapia and Casas 2009) (i.e. US$230.79). The values of consultation and attention from 2013 to 2020, are estimated based on prices obtained from the years 2006 and 2008, therefore a 3.20 % for the average inflation occurring in the 2008-2012 period was added to values from 2013 to 2020. It is assumed that the average inflation will remain without much variation until 2020. The cost of one lost business day is calculated from the average estimated salary. To obtain expenses attributable to the avoided cases by the decrease of PM10 levels, in the contemplated scenarios, the associated costs are multiplied by the avoided emergency cases.

Costs associated to cases related to hospitalization due to asthma

Hospitalization associated expenses are evaluated in this study with the costs generated by a day in bed rest (bed-day), inhalotherapy and the lost work days during the hospital stay. The average cost in Mexico of one bed-day in the hospital due to respiratory diseases in 2009 was US $340.78, inhalotherapy was US$3.33 (LaSalud 2009), while in the study made by the INER (Tapia and Casas 2009) there is an average expense of bed-day of US $333.75 (value obtained with 2007 prices). In the current study, we decided to use the US$340.78 value for bed-day and $3.33 for inhalotherapy. The values of bed-day and inhalotherapy used for the studied period from 2013 to 2020, are estimated based on prices obtained from 2009, therefore it is added a 3.20 % average inflation occurring in the 2008-2012 period. It is assumed that the average inflation will remain without much variation until 2020. The average days of hospital stay due to asthma in Mexicali, are estimated to be 4 disability days (source: General Hospital Pers. comm.). The IMSS is forced to pay for the subsidy disabilities of 60 % of the last contribution salary, starting from the fourth day in disability. (See article 96 of the general law of IMSS and article 42 section II of the Federal Labor Law). This means, that the employee will lose 100 % of the first three days of disability, leaving only 40 % of the remaining days. Thus, corresponding hospitalization expenses are calculated by: cost of bed-day + inhalotherapy expense + cost of: the first three days of hospitalization stay x work day cost + one remaining day of hospitalization stay x 40 % of one work day. Then the avoided hospitalization cases are multiplied by the projected hospitalization expenses to obtain costs attributable to the avoided cases for reducing the PM10 levels in the contemplated scenarios. Costs associated to restricted activity days To calculate the restricted activity days (RAD), the procedure established by Ostro (1994) was used, where it is established that 20 % of RAD pertains to the number of lost work days and the complement (80 % of RAD) is assessed as the third part of the salary. In this manner, the total expense due to daily RAD is calculated with the sum of the 20 % plus a third of the 80 % of the RAD that is a net 46 % of daily productivity. It is assumed that the annual average number of net RAD per employed person in Mexicali is 18.99 (see section “restricted activity days”). Then, the net avoided RAD are multiplied by the total generated cost per daily restricted activity according to the pertaining year; to obtain the costs attributable to avoided RAD for reducing PM10 levels, under one of the established scenarios. RESULTS The 98th percentiles of the 24 hour averages of PM10 levels in Mexicali have historically marked an increasing trend every year, just like the annual averages. In fact, we can see in figure 2, that since year 2000, Mexicali has not been able to comply with the NOM-025-SSA1-1993 (SSA 1993) for the annual PM10 averages. The social costs of the PHB estimated in each of the proposed scenarios are brought to present value using a discount rate of 3 % and presented as a percentage of state GDP reported in 2011. Table I shows the local ERF used in this research along with the affected population and incidence rates. The ERF are shown as percentage increments as increases in the number of units of respective pollutant. TABLE I Percentage changes per 10 mg/m3 of PM10 and their low, central and high levels, respectively 1Net incidence rate calculated after subtracting hospitalizations from emergency room visits. 2Incidence rate estimated using Imperial County, California, USA, data (see section “Hospitalization due to Asthma”). 3The days of Hospital stay and emergency room visits both from the economically active and employed population are also taken as restricted activity days, and are thus subtracted from the total of restricted activity days. 4The coefficient taken from literature (WHO 1996) was obtained by using the employed population of 18 to 65 year olds. For the present study the economically active and employed population is assumed to be in the 15 to 64 year old age range. 5Net rate calculated after subtracting the days of hospital stay and emergency room visits of employed population from the total of the restricted activity days of the employed population. 6Obtained by adapting PM2.5/ PM10of Mexicali = 0.72. Table II shows the low, central and high values of health and economic benefits estimated for each contemplated scenario in the study by effects of the reduction of PM10 in Mexicali. It shows the average avoided cases as well as the accumulated avoided cases for each health impact (2013-2020) with their corresponding expenses in US$ indicated in parentheses. The last line of the chart shows the expenses of each case (i.e. average and accumulated) as a state GDP percentage. The costs are brought to present value with a discount rate of 3 %. The GDP is the one reported in 2011 by INEGI (2012). For example, if PM10 pollution levels in Mexicali followed the projected trend over the studied years, but actions were carried out so that the average concentrations of each year succeed in being on the level of 50 μg/m3 marking the NOM- 025- SSA1-1993 (i.e. scenario b), the city would have a social benefit that would represent the central value of 3.56 % (low value: 2.50 %, high value: 5.37 %) of the state GDP reported in 2011. But, if these actions were about getting reductions of 8 % in the annual concentrations (i.e. scenario c), to reach in 2020 the level established by the norm; then the social benefit would represent a 2.13 % (1.50 %, 3.22 %) of the GDP. If the decision were about not conducting any action to decrease the levels of PM10, then the social expense would represent 5.59 % (3.92 %, 8.43 %) of the GDP.

Table II
Social benefits for reducing PM10 levels under the established scenarios for Mexicali, in the projected period of years (2013─2020).

Scenarios: (a) a total reduction per year starting from 2013 and until 2020. (b) a reduction to achieve every year the levels indicated by the NOM-025-SSA1-1993 of 50 µg/m3, starting from 2013 and until 2020. (c) a reduction of a ~8 % every year starting from 2013 to achieve in 2020 approximate levels as those indicated by the NOM-025-SSA1-1993 of 50 µg/m3.

1Discount rate of 3 %.

GDP: Gross domestic product.

DISCUSSION

We would recommend identifying the specific type of control measures that could be implemented in order to reduce PM10 levels in Mexicali. Perhaps separating by sectors where there should be more actions taken to reduce pollution accordingly to the type of source and pollutant (i.e. stationary sources, area, mobile). In the 2005 Mexicali emissions inventory (ERG 2009), the area sources are the ones contributing the most PM10 and PM2.5, whereas the CO is predominant in mobile sources in roads. There is a “ProAire Mexicali 2011-2020” where there was more than 40 actions proposed aiming to reduce air pollution. It would be advisable to assess the economic benefit that would result from the implementation of each of the proposed actions, in order to prioritize those showing the greatest social benefit impact resulting from the decrease of PM10. The results of the current project in a first approach shed light about the social benefits and costs that would generate, for carrying out or not those actions on achieving reductions in pollutant PM10 according to the 2020 projected scenarios.

If nothing were to be done to reduce PM10 in Mexicali and let it follow the trend projected to 2020, there would be a social cost of around $1659 million dollars. If the decision is to take actions to reduce PM10 levels by about 8 % annually, future projections of 2020 show that the social cost would be reduced to about$1026 million dollars. However, if the decision was to direct actions to reducing PM10 each year to have annual average levels established by the NOM-025-SSA1-1993, then the projections indicate that the social cost would decrease to about $601 million dollars. Then, some obligatory questions to answer would be: What kind of ProAire actions would be the most profitable to make, in terms of social benefit to be achieved according to the proposed scenarios? That is, once the economic benefits that can be obtained depending on the desired scenario are known, the most profitable ProAire actions would be those efficiently impacting the reduction of PM10 and/or other pollutants such as CO; the latter as long as the implementation cost does not exceed the estimated economic benefit in the contemplated scenarios. The results of this study could be underestimated, because the assessment was made using available data from the location, therefore, many aspects such as: the effects of long term exposure to the pollutants, or the inclusion of more diseases, or the use of the approximation “willingness to pay” instead of the approximation “human capital”, were not considered. Neither was considered the financial burden generated from the costs of medication used to control asthma, since it varies because it depends on the quality control of the disease; for example, an asthmatic patient, receiving appropriate care, has an average expense of$635 dollars, whereas the cost for a patient not receiving adequate care may add up to $10 582 dollars (SSA 2000). Another important aspect to mention is that the economic benefits in the current study were assessed under scenarios that consider PM10 levels established in NOM-025-SSA1-1993. However, this standard was replaced with NOM-025-SSA1-2014 on August 20, 2014 (SSA 2014), implying that the public health benefits would rise even more as the permissible limit values for the concentration of PM10 became more strict. The reason why the NOM-025-SSA1-1993 was used rather than the current one, was because the latter was published by the Official Journal of the Federation after this study had been completed. Nevertheless, the results provided by this study as an initial assessment, would justify the implementation of air pollution control measures in Mexicali. It is important to mention that there is a study (Torillo 2008) on the economic assessment of asthma in Mexicali due to the effect of reducing PM10 and O3, where the ERF are estimated through regression models known as error correction models. However, the methodology of these models apparently does not follow the conventional methodology with regards to ERF estimation that is widely used in Health Impact Assessment: there is not a previous baseline model controlling the confounding variables and lags. That is, two confounding variables are included in the model simultaneously (temperature and relative humidity) as well as two pollution variables PM10 and O3, not considering other important control variables such as days of the week, holidays, vacation days, analyzed years, trends, seasonality, among others. The metrics for these variables in the models are not consistent: 8-h mobile averages are used for the O3 variable while mobile averages are used for PM10 and daily averages are used for the other independent variables. The models used are linear, even when the dependent variable is discrete and with a very probable distribution function that is not normal. For the current study it was decided to maintain a conservative scenario, which is the reason why the ERF reported in the theses were not used. Nevertheless, the results from this study would come to reinforce the importance of the implementation of control measures that would aid to decrease the levels of air pollution in Mexicali. CONCLUSSIONS The study evaluates the benefits to health and the social cost for the reduction of the PM10 pollutant in Mexicali, under the three scenarios already described. Scenario A supposes total reduction of PM10 to level zero, meaning the absence of this pollutant in a way that it is presented as the baseline to determine the cost for doing nothing to reduce the concentrations. Scenario B proposes to reach and keep every year, starting 2013, PM10 concentrations at the levels indicated in the NOM-025-SSA1-1993. Both scenarios only establish a frame of reference for the study but its real viability was not evaluated. However, it would be feasible to reach a yearly reduction of 8 %, established by scenario C, (although its viability was not evaluated) by means of the conjunction of the actions established in the ProAire Mexicali 2011-2020. For that it is recommended to elaborate a cost-benefit study of the actions that ProAire establishes, determining the ones with the biggest effectiveness and profitability. However, the effective role of the three levels of government summed to the intergovernmental management and the social co-responsibility, are a requisite to strengthen the public policies that allow the PM10 to reach the levels established in the referred norm. Another important component is the high awareness of the main actors and decision makers; after years of bearing social struggles, health impacts, political costs and loss of economic competitiveness due to the bad air quality of Mexicali, the conditions are given to establish action of reduction of emissions suggested by scenario C. It is recommended the establishment of an intergovernmental organ that involves business sector, academia and civil society; that could be in charge of establishing and to impulse the diary that eventually allows implementing and evaluating the actions of Pro-Air of Mexicali. Supplementary material Alternative link Acknowledgements We thank the Environmental Protection Agency (EPA) and the Border Environment Cooperation Commission (BECC) for the founding of the project. This does not imply support for the findings and conclusions of the project. Also we thank the Autonomous University of Baja California (UABC) for its support to the project. References AMA (1999). Nuevas tablas de mortalidad oficiales. Entre Actuarios. Octubre-diciembre, Boletín de la Asociación Mexicana de Actuarios 4. 6 pp. CDPH (2015). California environmental health tracking program. Asthma data query. Database. California Department of Public Health [online]. http://cehtp.org/page/asthma/query 10/06/2015. CONAPO (2010). Sistema urbano nacional. Indicadores sociodemográficos de los municipios que integran las zonas metropolitanas, 2010. Base de datos. Consejo Nacional de Población [online]. http://www.conapo.gob.mx/work/models/CONAPO/Resource/1212/1/images/BaseDatosCompleta_Indicadores_SUN.xlsx 10/06/2016. ERG (2009). Project No. 3657.00.004.001. Mexicali emissions inventory. Eastern Research Group, Inc. Final Report. Sacramento, California, USA, 132 pp. Evans J., Spengler J., Levy J., Hammitt J., Suh H., Serrano-Trespalacios P., Rojas-Bracho L., Santos-Burgoa C., Rojas-Rodríguez H., Caballero-Ramírez M. and Castillejos M. (2000). Technical Report No. 10. Contaminación atmosférica y salud humana en la Ciudad de México. Massachusetts Institute of Technology-Integrated Program on Urban, Regional and Global Air Pollution. La Jolla, California, USA, 109 pp. Gallardo G., Arias A., González S.N. and Galindo G. (2007). Costos derivados de la atención médica del asma en un grupo de niños del Noreste de México. Rev. Alerg. Mex. 54, 82-85. Houtven G.V. and Cropper M.L. (1996). When is a life too costly to save? The evidence from U.S. environmental regulations. J. Environ. Econ. Manage. 30, 348-368. DOI: 10.1006/jeem.1996.0024. ICAR (1996). Inventario de emisiones de Mexicali. Ingeniería de Control Ambiental y Riesgo Industrial, S. de R.L.M.I. Final Report. Mexicali, Baja California, Mexico, 45 pp. LaSalud (2009). Reduce IMSS costo de hospitalización por afecciones respiratorias con inhaloterapia [online]. http://www.lasalud.mx/permalink/4111.html 16/06/2016. Zuk M., Tzintzun M.G. and Rojas L. (2007). Tercer almanaque de datos y tendencias de la calidad del aire en nueve ciudades mexicanas. IEPSA, Mexico City, Mexico, 116 pp. INECC (2011). Guía para evaluar los impactos en la salud por la instrumentación de medidas de control de la contaminación atmosférica. Manual. Instituto Nacional de Ecología y Cambio Climático. Mexico City, Mexico, 82 pp. INECC (2016). Sistema nacional de información de la calidad del aire. Dirección General del Centro Nacional de Investigación y Capacitación Ambiental. Database. Instituto Nacional de Ecología y Cambio Climático [online]. http://sinaica.inecc.gob.mx/index.html 15/01/2013. INEGI (1990). XI censo general de población y vivienda 1990. Database. Instituto Nacional de Estadística, Geografía e Informática [online]. 10/06/2016. INEGI (1994). Mexicali estado de Baja California. Cuaderno estadístico municipal. Edición 1994. Database. Instituto Nacional de Estadística, Geografía e Informática [online]. 10/06/2016. INEGI (1995). Conteo de población y vivienda 1995. Database. Instituto Nacional de Estadística, Geografía e Informática [online]. http://www.inegi.org.mx/est/contenidos/proyectos/ccpv/cpv1995/default.aspx 10/06/2016. INEGI (2000). XII censo general de población y vivienda 2000. Proceedings. Memoria-Instituto Nacional de Estadística, Geografía e Informática [online]. http://www3.inegi.org.mx/sistemas/biblioteca/ficha.aspx?upc=702825006519 10/06/2016. INEGI (2005). II Conteo de población y vivienda 2005. Database. Instituto Nacional de Estadística, Geografía e Informática [online]. http://www.inegi.org.mx/est/contenidos/proyectos/ccpv/cpv2005/default.aspx 10/06/2016. INEGI (2010). Conteo de población y vivienda 2010. Database. Instituto Nacional de Estadística, Geografía e Informática [online]. http://www.censo2010.org.mx/ 10/06/2016. INEGI (2012). Sistema de cuentas nacionales de México. Producto interno bruto por entidad federativa 2007-2011. Año base 2003. Primera versión. Database. Instituto Nacional de Estadística, Geografía e Informática [online]. http://www3.inegi.org.mx/sistemas/biblioteca/ficha.aspx?upc=702825003971 24/04/2011. Kidsdata (2016). A program of Lucile Packard Foundation for Children’s Health [online]. http://www.kidsdata.org 10/06/2016. Levy J.I., Hammitt J. and Spengler J. (2000). Estimating the mortality impacts of particulate matter: What can be learned from between-study variability? Environ. Health Perspect. 108, 109-117.DOI: 10.1289/ehp.00108109. McKinley G., Zuk M., Hojer M., Avalos M., González I., Hernández M., Iniestra R., Laguna I., Martinez M.A., Osnaya P., Reynales L.M., Valdés R. and Martínez J. (2003). The local benefits of global air pollution control in Mexico City. Instituto Nacional de Ecología and Instituto Nacional de Salud Pública. Final Report. Second phase of the integrated environmental strategies program in Mexico. Mexico City, Mexico, 175 pp. MEBC (2012). Mexicali la cuarta ciudad más contaminada del mundo: Organización Mundial de la Salud. Monitor Económico de Baja California. Año II, No. 456. Jueves 29 de noviembre [online]. http://monitoreconomico.org/ediciones/2012/11%20Noviembre/29%20noviembre%202012.pdf 29/11/2012. Ostro B.D. (1987). Air pollution and morbidity revisited: A specification test. J. Environ. Econ. Manage. 14, 87-98. DOI: 10.1016/0095-0696(87)90008-8. Ostro B.D. (1993). The association of air pollution and mortality: Examining the case for inference. Arch. of Environ. Health 48, 336-342. Ostro B.D. (1994). Policy research working paper 1301. Estimating health effects of air pollutants. A method with an application to Jakarta. World Bank. Washington, D.C., USA, 65 pp. OSHPD (2001). California inpatient discharge data. Database. Office of Statewide Health Planning and Development. California, USA [online]. https://www.hitpages.com/doc/5509496805785600/26#pageTop 07/06 /2016. OSHPD (2004). Imperial County Healthcare Information. California inpatient discharge data. Database. Office of Statewide Health Planning and Development, California, USA [online]. http://www.oshpd.ca.gov/HID/Perspectives/2004/imperial.pdf 07/06/2016. Reyna M.A., Bravo M.E., López R., Nieblas E.C. and Nava M.L. (2012). Relative risk of death from exposure to air pollutants: a short-term (2003-2007) study in Mexicali, Baja California, Mexico. Int. J. Environ. Health Res. 22, 370-86. DOI: 10.1080/09603123.2011.650153. Samet J.M., Zeger S.L., Dominici F., Curriero F., Coursac I., Dockery D.W., Schwartz J. and Zanobetti A. (2000). National morbidity, mortality, and air pollution study. Part II: Morbidity and mortality from air pollution in the United States. Health Effects Institute. Research report. Cambridge, Massachusetts, USA, 94 pp. Sánchez J.M., Valdés S. and Otro B. (1998). Los efectos en salud de la contaminación atmosférica por PM10 en Santiago. Estudios Públicos 69, 125-154. SDSU-UABC (2006). Binational surveillance of air-pollution-related illnesses in imperial county and the municipality of Mexicali. Graduate School of Public Health San Diego State University, San Diego, California, USA. Instituto de Ingeniería Universidad Autónoma de Baja California, Mexicali, Baja California, Mexico. California Center for Border and Regional Economic Studies San Diego State University, Calexico, California, USA. Research Report. Border environment research report (BERR). Mexicali, Baja California, Mexico, 27 pp. Schwartz J. and Dockery D.W. (1992). Particulate air pollution and daily mortality in Steubenville, Ohio. Am. J. Epidemiol. 135, 12-20. Schwartz J., Slater D., Larson T.V., Pierson W.E. and Koenig J.Q. (1993). Particulate air pollution and hospital emergency room visits for asthma in Seattle. Am. Rev. Respir. Dis. 147, 826-831. DOI: 10.1164/ajrccm/147.4.826. SSA (1993). Norma Oficial Mexicana NOM-025-SSA1-1993 y Modificación a la Norma Oficial Mexicana NOM-025-SSA1-1993. Criterios para evaluar la calidad del aire ambiente, con respecto a material particulado. Secretaría de Salud. Diario Oficial de la Federación. 26 de septiembre de 2005. SSA (2000). Entre el 5 y el 10 % de la población padece asma. Secretaría de Salud. Comunicado de prensa No. 537. 27 de diciembre. Mexico City, Mexico [online]. http://www.salud.gob.mx/unidades/dgcs/sala_noticias/comunicados/2000-12-27-537-ASMA.htm 20/03/2013. SSA (2014). Norma Oficial Mexicana NOM-025-SSA1-2014, Salud ambiental. Valores límite permisibles para la concentración de partículas suspendidas PM10 y PM2.5 en el aire ambiente y criterios para su evaluación. Secretaría de Salud. Diario Oficial de la Federación. 04 de agosto de 2014. Tapia A.M.and Casas G.A. (2009). Costos de atención y carga de enfermedad de pacientes asmáticos del Instituto Nacional de Enfermedades Respiratorias. Rev. Inst. Nal. Resp. Mex. 22, 182-189. Torillo E. (2008). Estimación del beneficio económico en materia de Salud por reducciones de Ozono y PM10 en Mexicali: un enfoque de costos evitados. Master Thesis. El Colegio de la Frontera Norte. Tijuana, Baja California, Mexico, 112 pp. USEPA (2003). BenMAP: Environmental benefits mapping and analysis program user’s manual. Abt Associates Inc. Office of Air Quality Planning and Standards. Manual. United States Environmental Protection Agency. Washington, D.C., USA, 415 pp. USEPA (2012). BenMAP: Environmental benefits mapping and analysis program. User’s manual appendices. Abt Associates Inc. Office of Air Quality Planning and Standards. Manual. United States Environmental Protection Agency. Washington, D.C., USA, 271 pp. WHO (1996). WHO/EHG/96.5. A methodology for estimating air pollution health effects. Office of Global and Integrated Environmental Health, Manual. World Health Organization. Geneva, Switzerland, 27 pp. WHO (2000). Quantification of health effects of exposure to air pollution, Regional Office for Europe, World Health Organization. Report. World Health Organozation. Bilthoven, Netherlands, 29 pp. WHO (2011). Urban outdoor air pollution in cities. Department of Public Health and Environment. Database. World Health Organization. [online]. http://www.who.int/phe/health_topics/outdoorair/databases/cities-2011/en/ 11/06/2016. Yang L., Wong C.M., Chan K.P., Chau P.Y., Ou C.Q., Chan K.H. and Malik J.S. (2009). Seasonal effects of influenza on mortality in a subtropical city. BMC Infect. Dis. 9, 133-141. DOI: 10.1186/1471-2334-9-133. Notes Author notes investigador.reyna@gmail.com Fig. 1 Annual historical averages of PM10 in Mexicali. Data provided by the Secretaría de Protección al Ambiente (SPA) of the Government of the State of Baja California. a According to the SPA, the value can be sub-estimated due to the fact that the indicator was obtained with information of less than 75 % of the valid data. b Values estimated by the authors with data from the manual stations provided by the SPA. (Source: Pers. comm). Fig. 2 Annual historical and projected PM10 averages for Mexicali. Data provided by the Secretaría de Protección al Ambiente (SPA) of the Government of the State of Baja California until the year 2012. From year 2013 are projections. For the purpose of such projections, the original values for 2005 and 2006 (i.e. values 210a and 254 of Fig. 1 ) were eliminated since they were extreme severe data and replaced by values estimated by the authors using regression models. Five regression models (i.e. exponential, linear, logarithmic, polynomial and power) were tested and was used that with the highest R2. a According to the SPA of the State Government, the value can be underestimated due to the fact that the indicator was obtained using information that was less than 75 % of the valid data.1 Data estimated by the authors using regression models. 2 Data estimated by the authors using data from the manual stations (COBACH, CONALEP, ITM, PROGRESO and CESPM) provided by the SPA. (Source: pers. comm.). (1) TABLE I Percentage changes per 10 mg/m3 of PM10 and their low, central and high levels, respectively 1Net incidence rate calculated after subtracting hospitalizations from emergency room visits. 2Incidence rate estimated using Imperial County, California, USA, data (see section “Hospitalization due to Asthma”). 3The days of Hospital stay and emergency room visits both from the economically active and employed population are also taken as restricted activity days, and are thus subtracted from the total of restricted activity days. 4The coefficient taken from literature (WHO 1996) was obtained by using the employed population of 18 to 65 year olds. For the present study the economically active and employed population is assumed to be in the 15 to 64 year old age range. 5Net rate calculated after subtracting the days of hospital stay and emergency room visits of employed population from the total of the restricted activity days of the employed population. 6Obtained by adapting PM2.5/ PM10of Mexicali = 0.72. Table II shows the low, central and high values of health and economic benefits estimated for each contemplated scenario in the study by effects of the reduction of PM10 in Mexicali. It shows the average avoided cases as well as the accumulated avoided cases for each health impact (2013-2020) with their corresponding expenses in US$ indicated in parentheses. The last line of the chart shows the expenses of each case (i.e. average and accumulated) as a state GDP percentage. The costs are brought to present value with a discount rate of 3 %. The GDP is the one reported in 2011 by INEGI (2012). For example, if PM10 pollution levels in Mexicali followed the projected trend over the studied years, but actions were carried out so that the average concentrations of each year succeed in being on the level of 50 μg/m3 marking the NOM- 025- SSA1-1993 (i.e. scenario b), the city would have a social benefit that would represent the central value of 3.56 % (low value: 2.50 %, high value: 5.37 %) of the state GDP reported in 2011. But, if these actions were about getting reductions of 8 % in the annual concentrations (i.e. scenario c), to reach in 2020 the level established by the norm; then the social benefit would represent a 2.13 % (1.50 %, 3.22 %) of the GDP. If the decision were about not conducting any action to decrease the levels of PM10, then the social expense would represent 5.59 % (3.92 %, 8.43 %) of the GDP.

Table II
Social benefits for reducing PM10 levels under the established scenarios for Mexicali, in the projected period of years (2013─2020).

Scenarios: (a) a total reduction per year starting from 2013 and until 2020. (b) a reduction to achieve every year the levels indicated by the NOM-025-SSA1-1993 of 50 µg/m3, starting from 2013 and until 2020. (c) a reduction of a ~8 % every year starting from 2013 to achieve in 2020 approximate levels as those indicated by the NOM-025-SSA1-1993 of 50 µg/m3.

1Discount rate of 3 %.

GDP: Gross domestic product.

Buscar:
Contexto
Descargar
Todas
Imágenes


text:
# Autocorrelation - cyclic vs acyclic

In terms of signal analysis (especially frequency decomposition) would it be bad to use cyclic autocorrelation (AC) instead of acyclic/linear AC?

My use case:
I'm trying to characterize events in a sensor stream by periodicity/frequency. At the moment, I'm using DFT spectrum (spectral density) and AC of sliding windows.
On the one handside it's algorithmically way cheaper to calculate the cyclic AC via FFT than acyclic (pad with zeroes and do double length FFT). However, the cyclic AC looks similar (though not equal) to the acyclic for the first N/2 lags and is self-symmetric - thus at most half the information of acyclic AC can be contained.

My motivation:
If the cyclic AC would be sufficient for detection, that would mean better performance compared to linear AC (in terms of memory usage and computation time). I then could also calculate AC inplace (better parallelization potential).

Rephrased question:
Is there any rule when cyclic AC should NOT be used?
Or perhaps, are there known cases, when cyclic AC might actually be better than linear AC?

• The cyclic autocorrelation of a PN sequence has ideal properties; the acyclic autocorrelation does not, except at length $7$ when a couple of PN sequences happen to be Barker sequences and so have acylic autocorrelation that is also ideal. – Dilip Sarwate May 22 '14 at 13:25
If the cyclic autocorrelation output has the properties you need in your application (e.g., it's peak-to-average ratio is high enough for picking out the peaks), then it's probably okay to substitute the cyclical correlation for the linear one. For example, consider a length 7 maximal-length sequence: $$\{+1, -1, -1, +1, +1, +1, -1\}$$ The cyclical and acyclical autocorrelation functions are plotted below.
Notice that both functions have a peak value of 7. The acyclical autocorrelation function has a largest off-peak magnitude of 3, whereas the cyclical autocorrelation function has a largest off-peak magnitude of 1. The peak-to-off-peak amplitude ratios are $7$ and $7/3$ for cyclical vs. acyclical, respectively. This means that the cyclical autocorrelation function is basically three times more peaky, making it much easier to detect. If the extra peak-to-off-peak amplitude ratio is not necessary for your application, it would be okay to correlate using linear correlation instead.

format: plain


text:
Pulsar Characteristic Age

• When astronomers measure the period, or spin-rate of the rotation of a pulsar, they find that pulsars are slowing down, usually at a very consistent rate. This "spin-down" is thought to be due to braking caused by the rotating pulsar's magnetic field, and this information can be used to determine an approximate age for the pulsar, known as the characteristic age.

A radio pulsar's characteristic age τ is usually defined as:

where P is the pulsar's period, and the dot represents the period derivative (the rate the pulsar is slowing). The characteristic age provides an approximate measure of a pulsar's true age, and the calculation is reasonably valid under three assumptions:

• The pulsar's initial spin period was very much smaller than that observed today.
• There is no magnetic field decay.
• The magnetic braking can be approximated by the energy loss a spinning dipolar magnet would experience in a perfect vacuum (in this case, the braking index, n = 3.)

The "correct" formula explicitly includes n and allows for a finite initial spin period P0 and derivative.

It is important to use the same units for period and age and to check the period derivative is dimensionless.

In 1054 AD the pulsar in the Crab supernova remnant was born.
Credit: NASA

Example

In 2007 the Crab pulsar had a period of 0.0331 sec and a period derivative of 4.22×10-13s/s. The characteristic age is around 1240 years. The supernova that produced the pulsar was in 1054 AD, yielding an age of ~950 years.


text:
# Synopsis: Outstanding in the field

#### Anomalous Meissner effect in pnictide superconductors

R. Prozorov, M. A. Tanatar, Bing Shen, Peng Cheng, Hai-Hu Wen, S. L. Bud’ko, and P. C. Canfield

Published November 18, 2010

When a superconductor in a magnetic field is cooled below its transition temperature, it expels the field. The details of how this Meissner effect occurs in real materials depend on the shape and physical properties of the sample, the type of superconductivity it exhibits, and the experimental conditions, but all superconductors are expected to behave as follows: As long as the external magnetic field is below a certain critical value, currents on the surface of the superconductor will form to cancel the field inside, but above this critical field magnetic flux will start to penetrate. The experimental signature of this field dependence is that the magnetization of the superconductor, which opposes the applied field, reaches a maximum at the critical field and then starts to decrease.

However, in a Rapid Communication published in Physical Review B, Ruslan Prozorov and collaborators from Ames Laboratory at Iowa State University, US, and the Institute of Physics of the Chinese Academy of Science in Beijing, show that two iron arsenide (pnictide) superconductors seem to defy this expected behavior. The researchers find that the magnetizations of $Ba(Fe0.926Co0.074)2As2$ and $Ba0.6K0.4Fe2As2$ continually increase in an approximately linear fashion, without reaching a maximum, even when the applied field far exceeds the estimated critical field for either material. Based on their results, Prozorov et al. suggest that the magnetic field suppresses magnetic scattering, yielding more resilient Cooper pairing in iron arsenide compounds than in other superconductors, but more experiments are needed to support or refute this proposal. – Matthew Eager


text:
# Victor Torres Perez: Combinatorial Principles without MA

Dear all,

The seminar meets on Wednesday October 3rd at 11:00 in the Institute of
Mathematics CAS, Zitna 25, seminar room, 3rd floor, front building.

Program: Victor Torres Perez — Combinatorial Principles without MA

Rado’s Conjecture (RC) in the formulation of Todorcevic is the statement
that every tree T that is not decomposable into countably many
antichains contains a subtree of cardinality ℵ_1 with the same property.
Todorcevic has shown the consistency of this statement relative to the
consistency of the existence of a strongly compact cardinal.
Todorcevic also showed that RC implies the Singular Cardinal Hypothesis,
a strong form of Chang’s Conjecture, the continuum is at most ℵ_2, the
negation of Box(θ) for every regular θ ≥ ω_2, etc. These implications
are very similar to the ones obtained from traditional forcing axioms
such as MM or PFA. However, RC is incompatible even with MA(ℵ_1).
In this talk we will take the opportunity to give an overview of our
results with different coauthors obtained in the last few years together
with recent ones, involving RC, certain weak square principles and
instances of tree properties. These new implications seem to continue
suggesting that RC is a good alternative to forcing axioms. We will
discuss to which extent this may hold true and where we can find some
limitations. We will end the talk with some open
problems and possible new directions.
For example, we will also discuss some recent results regarding the
P-ideal dichotomy (which can be consistent with the negation of MA(ℵ_1))
and square principles.

Best,
David

This site uses Akismet to reduce spam. Learn how your comment data is processed.


text:
# Text Classification with Naive Bayes in numpy

Goal: step by step build a model to classify text documents from newsgroups20 into their according categories.

You can download the accompanying jupyter notebook for this exercise from here and use the attached environment.yml to reproduce the used conda environment..

In this exercise we want to get to know a first classifier which is commonly referred to as “naive bayes”. But don’t get discouraged by the word “naive” as it refers to the circumstance that it uses some unrealistic but practical assumptions.

As the name suggests, in the Naive Bayes Classifier (NB) we make use of Bayes’ theorem to build a classifier. We are going to use the NB to map text documents to class labels $c \in \mathcal{C} = \{1, 2, \dots, M\}$. You could think of the example of “spam classification” from previous exercises in which we thought about a system which takes text of a mail and automatically tells whether it is spam (class “1”) or not spam (class “0”). Here, we look at a larger space of class labels ($M > 2$) and thus we’re looking at a multinomial case. In particular, we want to look at a dataset called “newsgroups20” with 20 possible classes.

Newsgroups20 can be found at http://qwone.com/~jason/20Newsgroups/ and is a popular entry classification problem in the text domain (the research field can be put into natural language processing). It contains roughly 20,000 documents from newsgroups evenly partitioned over 20 types of newsgroups. You can download the dataset on your system in your local notebook path by executing the following:

# shell scripts for downloading the data and placing it in a corresponding directory
!mkdir newsgroups
!curl -o newsgroups/news.tar.gz "http://qwone.com/~jason/20Newsgroups/20news-bydate.tar.gz"
# extract the files
!gzip -d < newsgroups/news.tar.gz | tar xf - --directory newsgroups

  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current
100 13.7M  100 13.7M    0     0  2948k      0  0:00:04  0:00:04 --:--:-- 3400k

Let’s have a look into what we’ve downloaded: we have a directory with two sub-directories “20news-bydate-test” and “20news-bydate-train”. Both contain for each newsgroup (or class) another subdirectory. So let’s first have a look into how they are named and if there are actual 20 of them:

# Reference the paths to our dataset
base_path_newsgroups = "./newsgroups/"
path_newsgroups_train = os.path.join(base_path_newsgroups, "20news-bydate-train")
path_newsgroups_test = os.path.join(base_path_newsgroups, "20news-bydate-test")

subdirectory_names = [file_name for file_name in os.listdir(path_newsgroups_train)]
print("We have <%s> categories/subdirs in our training folder" % len(subdirectory_names))
print("Named categories:")
print("\t- "+",\n\t- ".join(sorted(subdirectory_names)))

We have <20> categories/subdirs in our training folder
Named categories:
- alt.atheism,
- comp.graphics,
- comp.os.ms-windows.misc,
- comp.sys.ibm.pc.hardware,
- comp.sys.mac.hardware,
- comp.windows.x,
- misc.forsale,
- rec.autos,
- rec.motorcycles,
- rec.sport.baseball,
- rec.sport.hockey,
- sci.crypt,
- sci.electronics,
- sci.med,
- sci.space,
- soc.religion.christian,
- talk.politics.guns,
- talk.politics.mideast,
- talk.politics.misc,
- talk.religion.misc

Understanding what the output of a classifier (a class label) looks like is now quite simple: the classifier is a function returning us one of the 20 possible class labels. How we represent a particular class is sometimes different depending on what kind of classifier we look at. We already got to learn about the trivial and the random classifier. So let’s first count how many examplary documents we have available per class / category (btw those terms are often used interchangably).

categories = sorted(subdirectory_names)
row_format ="{:>5}" * (len(categories)+1)
print(row_format.format("Cat", *range(len(categories))))
samples_per_cat = [len(list(os.listdir(os.path.join(path_newsgroups_train, cat_name)))) for cat_name in categories]
print(row_format.format("#", *samples_per_cat))

  Cat    0    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16   17   18   19
#  480  584  591  590  578  593  585  594  598  597  600  595  591  594  593  599  546  564  465  377

Now we create a custom “document” and we will use that document and pass it to 1) a trivial classifier, 2) a random classifier and later 3) to our naive bayes classifier.

sample_document = "some document content"

We can now “build” our trivial classifier (as a callable function). It just returns the class with the most occurences (have a look at above categories and notice, that the eleventh category “rec.sport.hockey” contains most examples).

trivial_classifier = lambda document = None: np.argmax(samples_per_cat)
trivial_classifier(sample_document)

10

categories[trivial_classifier(sample_document)]

'rec.sport.hockey'

And similar to that we “build” a random classifier. Both of them are ignorant to the actual input:

total_samples = sum(samples_per_cat)
row_format ="{:>5}" * (len(categories)+1)
print(row_format.format("Cat", *range(len(categories))))
frequency_per_category = [samples_per_cat[cat]/total_samples for cat in range(len(categories))]
print(row_format.format("Frq", *[np.round(f, 2) for f in frequency_per_category]))

  Cat    0    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16   17   18   19
Frq 0.04 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.04 0.03

total_samples

11314

random_classifier = lambda document = None: np.argmax(np.random.multinomial(1, frequency_per_category))
random_classifier(sample_document)

1

You could now evaluate your classifier in terms of class-wise precision, recall and calculate their macro- and micro-averaged aggregations as an exercise. In case of the random classifier, you need to once run it over all documents of the test set and then evaluate the generated assignments.

The output of a classifier could also be represented slightly different as a “one-hot”-encoding. This is often quite suitable as we can better continue further calculations with it. Also, it has a nice probabilistic interpretation in our next section.

def one_hot(cat, num_cat):
assert cat <= num_cat
rep = np.zeros(num_cat)
rep[cat] = 1
return rep

print(one_hot(trivial_classifier(sample_document), len(categories)))
print(one_hot(random_classifier(sample_document), len(categories)))

[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]

trivial_classifier(sample_document)

10

But to get to Naive Bayes we now need to clarify how we represent the input to our classifier. Currently it is a string but we need to “learn” and “calculate” something based on this input. Using a “one-hot”-encoding is the simplest approach applied to that issue of getting a numeric representation of our input: we can represent each possible word (term) as a dimension (or category) and then represent a word as one-hot-encoded vector or a document as a matrix of such vectors (or as the sum of its one-hot-terms). Example: if we only consider the words “apple”, “peach” and “orange” then we can represent the first as $(1, 0, 0)$, the second as $(0, 1, 0)$ and the third one as $(0, 0, 1)$. A document containing “apple” and “orange” could then be a matrix of two rows with those vectors $\begin{pmatrix}1 & 0 & 0\\0 & 0 & 1\end{pmatrix}$ but there are also other representations possible.

# Naive Bayes

Each text document is represented as a “bag of words” - a feature vector $\mathbf{x} = (x_1, \dots, x_p)$ that has as many entries as there are words in our vocabulary $\mathcal{V}$, p = $|\mathcal{V}|$. Each dimension $i$ in $\mathbf{x}$ is associated with one word in the vocabulary. A value $x_i = 1$ if the word associated with dimension $i$ occurs in the document, and $x_i = 0$ if the word does not occur in the document.

Always keep some simple examples in mind: we can have $V = \{“apple”, “peache”, “orange”\}$ with $|V| = 3$ and we get $\mathbf{x} = (1, 0, 1)$ for a document containing “apple” and “orange”.

Before we can define the probabilistic model, we define the following random variables:

• $X_1, \dots, X_p$: the $p$ binary indicator random variables associated with attributes/words
• $C \in \{1, 2, \dots, M\}$: the random variable associated with the category of the document / class of the document

Given a dataset of text documents and class labels, our goal is to estimate:

$$P[C | X_1, \dots, X_p]$$

Then we can classify documents by simply finding the class $c$ that maximizes this conditional probability:

$$h(\mathbf{x}) = \argmax_{c \in \{1, 2, \dots, M\}} P[C=c | X_1=x_1, \dots, X_p=x_p]$$

We can’t estimate this quantity directly from our dataset, but by using Bayes theorem we can write:

$$P[C | X_1, \dots, X_p] = \frac{P[X_1, \dots, X_p | C] \cdot P[C]}{P[X_1, \dots, X_p]}$$

That’s already helpful. We can estimate $P[C]$ by calculating the relative frequencies of classes over all documents (regardless of their content). Remember, we already did this above: for the first category “alt.atheism” we get:

frequency_per_category[0]

0.04242531377054976

However, the joint distribution $P[X_1, \dots, X_p | C]$ is a little annoying. Even for tiny vocabulary sizes (p $\approx 100$) we would have to store and estimate $2^{100}$-1 parameters.
The Naive Bayes Classifier is called “Naive” because it makes one very strong assumption about the observation variables $X_1, \dots, X_p$: All $X_i$ are conditionally independent given the class $C$:

$$P[X_1, \dots, X_p | C] \cdot P[C] = P[C] \cdot \prod_{i=1}^p P[X_i | C]$$

Sanity Check: Let's say we have documents from a newsgroup in which people discuss about politics. When we accept this independence assumption we acknowledge: Knowing that a document contains "Donald" does not alter our belief about whether this same document would also contain the word "Trump". $P[X_{'Trump'} = 1 | C = 2] = P[X_{'Trump'} = 1 | X_{'Donald'} = 1, C = 2]$. These independence assumptions will most likely be violated by any newsgroup. However, they lead to a very simple model and, in practice, the classifer yields very accurate predictions anyway.

Now the calculations are straight-forward. The terms $P[X_i | C]$ can be estimated by calculating the relative frequency of word $i$ per class $c$, resulting in a table of $p \cdot M$ conditional probabilities. And finally, we can make predictions by calculating:

$$h(\mathbf{x}) = \argmax_{c \in \{1, 2, \dots, M\}} P[C] \cdot \prod_{i=1}^p P[X_i | C]$$

NOTE: We completely ignored the denominator $P[X_1, \dots, X_p]$ in the maximization. This ignorance is justified by the fact that the denominator is a function of the data and the data is fixed. Since the data is fixed, the denominator is a constant. As we are only interested in the class $c$ that maximizes the above quantity, the constant denominator is irrelevant. It would be a different story if we were interested in the class probabilities, then certainly would we have to take $P[X_1, \dots, X_p]$ into account. Otherwise $P[C | \cdot]$ wouldn’t be a proper probability mass function (pmf).

## Naive Bayes Pseudocode

Now we want to implement Naive Bayes and learn how to translate parts of the mathematical expressions into computational estimates of e.g. $P[C]$ or $P[X_i|C]$.

#### TrainMultiNomialNB($\mathbb C$,$\mathbb D$)

$V \leftarrow extractVocabulary(\mathbb D)$
$N \leftarrow countDocs(\mathbb D)$
for $c \in \mathbb C$:
$N_c \leftarrow countDocsInClass(\mathbb D, c)$
$prior[c] \leftarrow \frac{N_c}{N}$
$text_c \leftarrow concatenateTextOfAllDocsInClass(\mathbb D, c)$
for $t \in V$:
$T_{ct} \leftarrow countTokensOfTerm(text_c,t)$
for $t \in V$:
$condprob[c][t] \leftarrow \frac{T_{ct} + 1}{\sum_{t’}(T_{ct’} + 1)}$
return $V,prior,condprob$

#### ApplyMultinomialNB($\mathbb C,V,prior,condprob,d$)

$W \leftarrow extractTokensFromDoc(V,d)$
for $c \in \mathbb C$:
$score[c] \leftarrow log(prior[c])$
for $t \in W$:
$score[c] += log(condprob[c][t])$
return $argmax_{c \in \mathbb C} score[c]$

# Implementing Multinomial Naive Bayes

What we’re now going to do is to fill below class template. We’re going to use a class as an instantiated object of it will first be “trained” on a given dataset (the newsgroups20). During this process we obtain for example how many categories we need and how large $p$ (the number of terms) is. At the end we will have a function (method of the object) which will return a class-prediction for a given document.

class NaiveBayesClassifier:
def train(self, path):
pass  # This has to be implemented

def test(self, path):
pass  # This has to be implemented

def predict(self, document):
pass  # This has to be implemented

Let’s first read a single sample of one category of our newsgroups training corpus:

exemplary_category = "comp.graphics"
examplary_category_doc = "37915"
exemplary_path = os.path.join(path_newsgroups_train, exemplary_category, examplary_category_doc)
exemplary_path

'./newsgroups/20news-bydate-train/comp.graphics/37915'

with codecs.open(exemplary_path, encoding='latin1') as file:

['From: kph2q@onyx.cs.Virginia.EDU (Kenneth Hinckley)\n',
'Subject:   VOICE INPUT -- vendor information needed\n',
'Organization: University of Virginia\n',
'Lines: 27\n',
'\n',
'\n',
'Hello,\n',
'     I am looking to add voice input capability to a user interface I am\n',
'developing on an HP730 (UNIX) workstation.  I would greatly appreciate \n',
'information anyone would care to offer about voice input systems that are \n',
'easily accessible from the UNIX environment. \n',
'\n',
'     The names or adresses of applicable vendors, as well as any \n',
'\n',
'     Please respond via email; I will post a summary if there is \n',
'sufficient interest.\n',
'\n',
'\n',
'Thanks,\n',
'Ken\n',
'\n',
'\n',
"P.S.  I have found several impressive systems for IBM PC's, but I would \n",
'like to avoid the hassle of purchasing and maintaining a separate PC if \n',
'at all possible.\n',
'\n',
'-------------------------------------------------------------------------------\n',
'Ken Hinckley (kph2q@virginia.edu)\n',
'University of Virginia \n',
'Neurosurgical Visualization Laboratory\n',
'-------------------------------------------------------------------------------\n']

def tokenize_str(doc):
return re.findall(r'\b\w\w+\b', doc) # return all words with #characters > 1

tokenize_str("I constructed this example.")

['constructed', 'this', 'example']

def tokenize_file(doc_file):
with codecs.open(doc_file, encoding='latin1') as doc:

print(tokenize_file(exemplary_path)[:10])

['hello', 'am', 'looking', 'to', 'add', 'voice', 'input', 'capability', 'to', 'user']

### Iterating over the training set to “build the vocabulary”

We need to iterate over our whole dataset and read it into memory in a way that we can use its information later to compute conditional probabilities. Let’s have a first look over this loop: we are iterating over all categories (outer loop) and all documents of each category (inner loop). As a neat trick to look into these dynamics, we add zip(a_list, range(5)) to iterate only over the first five categories and within each category over the first five documents and see how many terms are in there:

for class_name, _class_idx in zip(os.listdir(path_newsgroups_train), range(5)):
print(class_name)
path_class = os.path.join(path_newsgroups_train, class_name)
for doc_name, _idx in zip(os.listdir(path_class), range(3)):
terms = tokenize_file(os.path.join(path_class, doc_name))
print("\tdoc<%s> #terms=%s, First five terms:" % (doc_name, len(terms)), terms[:5])

comp.sys.ibm.pc.hardware
doc<60529> #terms=118, First five terms: ['in', 'article', '1993apr19', '132748', '18044']
doc<60402> #terms=105, First five terms: ['okay', 'here', 'is', 'my', 'configuration']
doc<60421> #terms=48, First five terms: ['have', '486', 'machine', 'with', 'drive']
rec.autos
doc<103011> #terms=96, First five terms: ['monthly', 'posting', 'regarding', 'the', 'buick']
doc<102874> #terms=189, First five terms: ['in', 'article', 'c5jnk3', 'jkt', 'news']
doc<103406> #terms=36, First five terms: ['in', 'article', 'apr16', '215151', '28035']
comp.windows.x
doc<67169> #terms=128, First five terms: ['usually', 'when', 'start', 'up', 'an']
doc<67113> #terms=285, First five terms: ['the', 'multi', 'lingual', 'archives', 'at']
doc<66416> #terms=133, First five terms: ['does', 'anyone', 'have', 'any', 'good']
sci.space
doc<61050> #terms=357, First five terms: ['in', 'article', '1993apr21', '212202', 'aurora']
doc<61236> #terms=70, First five terms: ['can', 'tell', 'you', 'that', 'when']
doc<61055> #terms=99, First five terms: ['related', 'question', 'which', 'haven', 'given']
talk.politics.misc
doc<178364> #terms=244, First five terms: ['in', 'article', '1993apr17', '111713', '4063']
doc<178413> #terms=355, First five terms: ['in', 'article', '16apr199317110543', 'rigel', 'tamu']
doc<178587> #terms=554, First five terms: ['in', 'article', '1993apr15', '170731', '8797']

Now let’s collect that information into a pythonic structure in which we have direct access to

• a vocabulary which contains each term such as {“hello”: 2, “world”: 3}
• the size of the vocabulary
• the number of total documents
• the number of documents per category
• the number of total terms per category
• the terms within each category (after looking into all its documents)
min_count = 1
vocabulary = {}
num_docs = 0
classes = {}
priors = {}
conditionals = {}

num_docs = 0
for class_name in os.listdir(path_newsgroups_train):
classes[class_name] = {"doc_counts": 0, "term_counts": 0, "terms": {}}
print(class_name)
path_class = os.path.join(path_newsgroups_train, class_name)
for doc_name in os.listdir(path_class):
terms = tokenize_file(os.path.join(path_class, doc_name))

num_docs += 1
classes[class_name]["doc_counts"] += 1

# build vocabulary and count terms
for term in terms:
classes[class_name]["term_counts"] += 1
if not term in vocabulary:
vocabulary[term] = 1
classes[class_name]["terms"][term] = 1
else:
vocabulary[term] += 1
if not term in classes[class_name]["terms"]:
classes[class_name]["terms"][term] = 1
else:
classes[class_name]["terms"][term] += 1

# remove terms with frequency < min_count
vocabulary = {k:v for k,v in vocabulary.items() if v > min_count}

print("Term counts per class:")
print({name: classes[name]["term_counts"] for name in classes})
print("Vocabulary size:", len(vocabulary))

comp.sys.ibm.pc.hardware
rec.autos
comp.windows.x
sci.space
talk.politics.misc
alt.atheism
sci.med
sci.crypt
rec.motorcycles
rec.sport.hockey
rec.sport.baseball
misc.forsale
comp.os.ms-windows.misc
comp.graphics
comp.sys.mac.hardware
talk.politics.guns
talk.politics.mideast
talk.religion.misc
sci.electronics
soc.religion.christian
Term counts per class:
{'comp.sys.ibm.pc.hardware': 104937, 'rec.autos': 118241, 'comp.windows.x': 159094, 'sci.space': 161284, 'talk.politics.misc': 190031, 'alt.atheism': 150831, 'sci.med': 160614, 'sci.crypt': 208234, 'rec.motorcycles': 107499, 'rec.sport.hockey': 157104, 'rec.sport.baseball': 115492, 'misc.forsale': 68699, 'comp.os.ms-windows.misc': 264777, 'comp.graphics': 116182, 'comp.sys.mac.hardware': 90823, 'talk.politics.guns': 185157, 'talk.politics.mideast': 260529, 'talk.religion.misc': 121243, 'sci.electronics': 107404, 'soc.religion.christian': 205727}
Vocabulary size: 60738

### Calculating the prior probabilities and the conditional probabilities for each term

To now come up with the prior probability of the category $P(C)$ and the conditional probabilities $P(X_i|C)$ we estimate the first one simply by taking its relative document frequency. However, we apply a neat trick by working in logarithmic space (this is why in probability theory and machine learning you will often see the term log-probability (or logprob)). This costs us additional effort when computing the prior and conditional probabilities in our training phase but later we only have to sum them up - which is equivalent to multiplying them: $log(P(X_1|C))+log(P(X_2|C)) = log(P(X_1|C)*P(X_2|C))$

categories[0], classes[categories[0]]["doc_counts"]

('alt.atheism', 480)

num_docs

11314

$P(C=0)$

classes[categories[0]]["doc_counts"]/num_docs

0.04242531377054976

$log(A/B) = log(A)-log(B)$

math.log(classes[categories[0]]["doc_counts"])-math.log(num_docs)

6.173786103901937

math.log(1/len(vocabulary.values()))+math.log(5/len(vocabulary.values()))

-20.419211709219514

See how for this exemplary category we can calculate the prior $P(C)$ in “normal” space and in log space:

classes[categories[0]]["doc_counts"]/num_docs, math.log(classes[categories[0]]['doc_counts']) - math.log(num_docs)

(0.04242531377054976, -3.160010072001165)

len(vocabulary.values())

60738

categories[0]

'alt.atheism'

list(classes[categories[0]]['terms'].keys())[8], list(classes[categories[0]]['terms'].values())[8]

('chimps', 15)

sum(classes[categories[0]]['terms'].values())

150831

Now our training phase consists in calculating this prior probability for all categories and then calculate the conditionals by counting all unique terms in all documents of the particular category. Note, that we apply another technique here: laplacian smoothing by assuming that each term occurs once by default. What would happen if we would not apply it?

We could possibly have no occurences of a term within a document and therefore we would divide by zero or try to evaluate the logarithm of zero.

logspace_num_docs = math.log(num_docs)
for cn in classes:
# calculate priors: doc frequency
priors[cn] = math.log(classes[cn]['doc_counts']) - logspace_num_docs

# calculate conditionals
conditionals[cn] = {}
cdict = classes[cn]['terms']
terms_in_class = sum(cdict.values())

for term in vocabulary:
t_ct = 1.  # We assume each term to occur in each document at least once - smoothing!
t_ct += cdict[term] if term in cdict else 0.
conditionals[cn][term] = math.log(t_ct) - math.log(terms_in_class + len(vocabulary))

print(math.e**conditionals[categories[0]]["disturbing"], math.e**conditionals["sci.crypt"]["disturbing"])
print(conditionals[categories[0]]["disturbing"], conditionals["sci.crypt"]["disturbing"])

3.62350349305737e-06 3.599506867559147e-05
-12.52806918439661 -10.232128610124521

### Calculating a score based on the conditional probabilities

print("Testing <%s>" % exemplary_path)
token_list = tokenize_file(exemplary_path)
token_list[:5]

Testing <./newsgroups/20news-bydate-train/comp.graphics/37915>

conditionals[categories[0]]['hello']

-11.163694176730045

np.array(frequency_per_category)

array([0.04242531, 0.05161747, 0.05223617, 0.05214778, 0.05108715,
0.05241294, 0.05170585, 0.05250133, 0.05285487, 0.05276648,
0.05303164, 0.05258971, 0.05223617, 0.05250133, 0.05241294,
0.05294326, 0.04825879, 0.04984974, 0.04109952, 0.03332155])

np.array([priors[p] for p in priors])

array([-2.95367364, -2.94691686, -2.94860178, -2.94860178, -3.19175877,
-3.16001007, -2.94691686, -2.94523477, -2.94020542, -2.93686652,
-2.94187906, -2.96218433, -2.95198016, -2.96389519, -2.97422231,
-3.0311772 , -2.99874192, -3.40155099, -2.95198016, -2.93853458])

scores = {}
for class_num, class_name in enumerate(classes):
scores[class_name] = priors[class_name]
for term in token_list: # ['hello', 'am', 'looking', ..]
if term in vocabulary: # ['hello', 'am', ..]
scores[class_name] += conditionals[class_name][term]
print(scores)
print("Result:")
print(max(scores, key=scores.get))

{'comp.sys.ibm.pc.hardware': -833.0541069962882, 'rec.autos': -862.3330497046505, 'comp.windows.x': -817.2169384851644, 'sci.space': -849.5240956512641, 'talk.politics.misc': -871.7723219710546, 'alt.atheism': -882.9604871670316, 'sci.med': -843.1841184508672, 'sci.crypt': -835.3925630963009, 'rec.motorcycles': -877.4131260209258, 'rec.sport.hockey': -922.3793867458648, 'rec.sport.baseball': -902.1850471341099, 'misc.forsale': -874.320329015849, 'comp.os.ms-windows.misc': -915.1018517113815, 'comp.graphics': -810.6187571263675, 'comp.sys.mac.hardware': -843.6084671998246, 'talk.politics.guns': -889.6462333560933, 'talk.politics.mideast': -899.2959096703631, 'talk.religion.misc': -893.1896022184637, 'sci.electronics': -835.1713581601585, 'soc.religion.christian': -879.8719468915807}
Result:
comp.graphics

### Putting it all together

class NaiveBayesClassifier:
def __init__(self, min_count=1):
self._min_count = min_count
self._vocabulary = {}
self._num_docs = 0
self._classes = {}
self._priors = {}
self._conditionals = {}

def _build_vocab(self, path_base):
self._num_docs = 0
for class_name in os.listdir(path_base):
self._classes[class_name] = {"doc_counts": 0, "term_counts": 0, "terms": {}}
print(class_name)
path_class = os.path.join(path_base, class_name)
for doc_name in os.listdir(path_class):
terms = tokenize_file(os.path.join(path_class, doc_name))

self._num_docs += 1
self._classes[class_name]["doc_counts"] += 1

# build vocabulary and count terms
for term in terms:
self._classes[class_name]["term_counts"] += 1
if not term in self._vocabulary:
self._vocabulary[term] = 1
self._classes[class_name]["terms"][term] = 1
else:
self._vocabulary[term] += 1
if not term in self._classes[class_name]["terms"]:
self._classes[class_name]["terms"][term] = 1
else:
self._classes[class_name]["terms"][term] += 1
# remove terms with frequency < min_count
self._vocabulary = {k:v for k,v in self._vocabulary.items() if v > self._min_count}

def train(self, path):
self._build_vocab(path)

logspace_num_docs = math.log(self._num_docs)
for cn in self._classes:
# calculate priors: doc frequency
self._priors[cn] = math.log(self._classes[cn]['doc_counts']) - logspace_num_docs

# calculate conditionals
self._conditionals[cn] = {}
cdict = self._classes[cn]['terms']
terms_in_class = sum(cdict.values())

for term in self._vocabulary:
t_ct = 1.  # We assume each term to occur in each document at least once - smoothing!
t_ct += cdict[term] if term in cdict else 0.
self._conditionals[cn][term] = math.log(t_ct) - math.log(terms_in_class + len(self._vocabulary))

def test(self, path_test):
true_y = []
pred_y = []

for c in self._classes:
print(c)
for f in os.listdir(os.path.join(path_test, c)):
doc_path = os.path.join(path_test, c, f)
_, result_class = self.predict_path(doc_path)
pred_y.append(result_class)
true_y.append(c)
return true_y, pred_y

def predict_path(self, path: str):
token_list = tokenize_file(path)
return self._predict(token_list)

def predict(self, document: str):
token_list = tokenize_str(document)
return self._predict(token_list)

def _predict(self, token_list: list):
scores = {}
for c in self._classes:
scores[c] = self._priors[c]
for term in token_list:
if term in self._vocabulary:
scores[c] += self._conditionals[c][term]
return scores, max(scores, key=scores.get)

clf_naivebayes = NaiveBayesClassifier()
clf_naivebayes.train(path_newsgroups_train)

comp.sys.ibm.pc.hardware
rec.autos
comp.windows.x
sci.space
talk.politics.misc
alt.atheism
sci.med
sci.crypt
rec.motorcycles
rec.sport.hockey
rec.sport.baseball
misc.forsale
comp.os.ms-windows.misc
comp.graphics
comp.sys.mac.hardware
talk.politics.guns
talk.politics.mideast
talk.religion.misc
sci.electronics
soc.religion.christian

ty, py = clf_naivebayes.test(path_newsgroups_test)

comp.sys.ibm.pc.hardware
rec.autos
comp.windows.x
sci.space
talk.politics.misc
alt.atheism
sci.med
sci.crypt
rec.motorcycles
rec.sport.hockey
rec.sport.baseball
misc.forsale
comp.os.ms-windows.misc
comp.graphics
comp.sys.mac.hardware
talk.politics.guns
talk.politics.mideast
talk.religion.misc
sci.electronics
soc.religion.christian

from sklearn import metrics

print('Accuracy: %2.2f %%' % (100. * metrics.accuracy_score(ty, py)))

Accuracy: 75.73 %

### Other utility functions

#
# Code from https://scikit-learn.org/stable/auto_examples/model_selection/plot_confusion_matrix.html#sphx-glr-auto-examples-model-selection-plot-confusion-matrix-py
#

def plot_confusion_matrix(cm, classes,
normalize=False,
title='Confusion matrix',
cmap=plt.cm.Blues):
"""
This function prints and plots the confusion matrix.
Normalization can be applied by setting normalize=True.
"""
if normalize:
cm = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis]
print("Normalized confusion matrix")
else:
print('Confusion matrix, without normalization')

#print(cm)

plt.imshow(cm, interpolation='nearest', cmap=cmap)
plt.title(title)
plt.colorbar()
tick_marks = np.arange(len(classes))
plt.xticks(tick_marks, classes, rotation=90)
plt.yticks(tick_marks, classes)

fmt = '.2f' if normalize else 'd'
thresh = cm.max() / 2.
for i, j in itertools.product(range(cm.shape[0]), range(cm.shape[1])):
plt.text(j, i, format(cm[i, j], fmt),
horizontalalignment="center",
color="white" if cm[i, j] > thresh else "black")

plt.ylabel('True label')
plt.xlabel('Predicted label')
plt.tight_layout()

# Compute confusion matrix
cnf_matrix = metrics.confusion_matrix(ty, py)
np.set_printoptions(precision=2)

# Plot non-normalized confusion matrix
plt.figure(figsize=(12,12))
plot_confusion_matrix(cnf_matrix, classes=clf_naivebayes._classes.keys(), title='Confusion matrix, without normalization')
plt.show()

Confusion matrix, without normalization

# SKlearn implementation of NaiveBayes

from sklearn.datasets import load_files
from sklearn import feature_extraction
from sklearn import metrics
from sklearn.naive_bayes import MultinomialNB

twenty_train = load_files(path_newsgroups_train, encoding='latin1')

vectorizer = feature_extraction.text.CountVectorizer()#(max_features=len(clf.vocabulary))
train_X, train_y = vectorizer.fit_transform(twenty_train.data), twenty_train.target
test_X, test_y = vectorizer.transform(twenty_test.data), twenty_test.target
print(train_X.shape, train_y.shape, test_X.shape, test_y.shape)

(11314, 130107) (11314,) (7532, 130107) (7532,)

clf = MultinomialNB()
clf.fit(train_X, train_y)
pred_y = clf.predict(test_X)

print('Accuracy: %2.2f %%' % (100. * metrics.accuracy_score(test_y, pred_y)))

Accuracy: 77.28 %

### Can we do better?

from sklearn.pipeline import Pipeline
from sklearn.feature_extraction.text import CountVectorizer
from sklearn.feature_extraction.text import TfidfTransformer
from sklearn.model_selection import GridSearchCV

# Warning: this needs some computation time!
clf = Pipeline([('vect', CountVectorizer()),
('tfidf', TfidfTransformer()),
('clf', MultinomialNB()),
])

parameters = {'vect__ngram_range': [(1, 1), (1, 2)],
'tfidf__use_idf': (True, False),
'clf__alpha': (1.0, 0.1),
}

gs_clf = GridSearchCV(clf, parameters, n_jobs=-1)
gs_clf.fit(twenty_train.data, twenty_train.target)
pred = gs_clf.predict(twenty_test.data)

print('Accuracy: %2.2f %%' % (100. * metrics.accuracy_score(twenty_test.target, pred)))

for param_name in sorted(parameters.keys()):
print(param_name,":", gs_clf.best_params_[param_name])

Accuracy: 81.78 %
clf__alpha : 0.1
tfidf__use_idf : True
vect__ngram_range : (1, 2)

# Warning: this needs some computation time!
from sklearn.linear_model import SGDClassifier
#SGDClassifier with hinge loss gives a linear SVM

clf = Pipeline([('vect', CountVectorizer()),
('tfidf', TfidfTransformer()),
('clf', SGDClassifier(loss='hinge', penalty='l2', alpha=1e-3, random_state=42)),
])

parameters = {'vect__ngram_range': [(1, 1), (1, 2)],
'tfidf__use_idf': (True, False),
'clf__alpha': (1e-2, 1e-3),
}

gs_clf = GridSearchCV(clf, parameters, n_jobs=-1)
gs_clf.fit(twenty_train.data, twenty_train.target)
pred = gs_clf.predict(twenty_test.data)

print('Accuracy: %2.2f %%' % (100. * metrics.accuracy_score(twenty_test.target, pred)))

for param_name in sorted(parameters.keys()):
print(param_name,":", gs_clf.best_params_[param_name])

## Authorship

This introduction was partly created by

in context of the Advanced Data Science lecture at the University of Passau and with examples taken from the scikit-learn documentation.