Regresion No Parametrica en R

REGRESION
NO PARAMETRICA
EN R
Trabajo Fin de Master
Master en Estadstica Aplicada
Autora: Nisa Boukichou Abdelkader

Tutora: Mara Dolores Martnez Miranda
Indice general
Pr
ologo
1. Introducci
on
1.1. Antecedentes . . . . . . . . . . . . . . . . . . . . . .
1.2. Estimacion del modelo de Regresion No Parametrica
1.3. Regresion Polinomial Local . . . . . . . . . . . . . . .
1.4. Metodos de seleccion de la complejidad del modelo .
1.5. Extension multivariante . . . . . . . . . . . . . . . .
1.5.1. El problema de la dimensionalidad . . . . . .
1.5.2. Modelos aditivos no parametricos . . . . . . .
1.6. Seleccion del parametro ancho de banda . . . . . . .
2. Software disponible en R
2.1. Introduccion . . . . . . .
2.2. Libro KernSmooth . . . .
2.3. Libro Locpol . . . . . . .
2.4. Libro Locfit . . . . . . .
2.5. Libro sm . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
5
5
6
10
13
14
16
16
17
.
.
.
.
.
19
19
20
30
63
68
3. Aplicaci
on pr
actica
87
3.1. Estudio con datos reales . . . . . . . . . . . . . . . . . . . . . . . . . 87
3.2. Estudio con datos simulados . . . . . . . . . . . . . . . . . . . . . . . 94
Bibliografa
99
iii
Pr
ologo
Dado el rapido avance que ha experimentado la Estadstica Computacional en
las u
ltimas decadas, hoy en da podemos agradecerle el desarrollo de diversos campos dentro de la Estadstica, que eran impensables dado que requeran costosos
procedimientos de calculo. Un ejemplo de este tipo lo constituyen los enfoques no
parametricos del Analisis de Regresion.
Las tecnicas de Regresion No Parametrica logran una mejor adaptacion a los
datos disponibles, mediante la obtencion de estimaciones mas proximas a la curva
de regresion subyacente. Esto es posible usando la informacion suministrada directamente desde los datos, sin formular rgidos modelos parametricos.
En este trabajo nuestro objetivo ha sido el de explorar las tecnicas de regresion
no parametrica mas habituales y las capacidades que R incorpora actualmente para
su aplicacion practica. En este sentido el trabajo se ha estructurado en tres captulos.
El primero tiene como finalidad establecer los elementos teoricos fundamentales de
la regresion no parametrica, desde la propia formulacion del modelo. De este modo
para un problema general de regresion se definen dos vas se solucion. Una sera la
regresion parametrica o clasica que presenta la ventaja de ser mas sencilla y menos
costosa desde el punto de vista computacional, pero que suele ser muy poco flexible
y de difcil adaptacion en situaciones complejas. Paralelamente y no necesariamente
en contraposicion (puesto que ambas pueden ir de la mano) la denominada regresion
no parametrica. De esta u
ltima destacamos fundamentalmente su flexibilidad, ya que
permite una mejor adaptacion a diversas situaciones y problemas, si bien requiere
de un elevado coste computacional y una mayor complejidad desde el punto de vista
teorico.
Una vez definido el contexto general y establecidas las caractersticas particulares
que perfilan el problema de regresion no parametrica frente a los planteamientos
clasicos parametricos, se procede a analizar algunas de las mas relevantes tecnicas de
este tipo. El tratamiento que se ha hecho de dichos metodos en este trabajo, ha sido
dirigido fundamentalmente hacia la practica y concreto la practica con el software
R. De este modo no se ha profundizado en aspectos teoricos de complejidad como
son los estudios de tipo asintotico. Bajo tal perspectiva se han explorado metodos
univariantes y multivariantes, perfilandose los denominados metodos de regresion
polinomial local como una buena solucion, dadas sus buenas propiedades teoricas y
1
PROLOGO
sus deseables caractersticas de interpretabilidad y sencillez en la practica.

De forma sencilla se plantea tambien en este primer captulo el conocido problema
de la dimensionalidad. Desde dicha motivacion se introducen metodos que permiten
salvar dicho problema, como son los modelos de regresion aditivos no parametricos. Dichos modelos se caracterizan fundamentalmente porque la naturaleza de los
efectos de las variables explicativas sobre la variable de respuesta se considera de
forma individual. Esto obviamente permite ganar en simplicidad y tambien en interpretabilidad.
Asociado a los metodos de regresion no parametrica (univariantes o multivariantes) se introduce uno de los problemas tecnicos cruciales en la practica, la eleccion
del parametro de suavizado o ancho de banda que define la complejidad del modelo.
Desde el punto de vista teorico se formula el problema de seleccion y se perfilan los
distintos metodos dise
nados para su seleccion automatica. En concreto se distingue
entre los metodos basados en la metodologa plug-in, los basados en el criterio de
validacion cruzada (cross-validation) y los procedimientos basados en Bootstrap.
Una vez expuestos los elementos teoricos necesarios comienza el captulo dos,
donde se hace un estudio profundo de los aspectos computacionales asociados a
dichos metodos. El software analizado es el entorno de analisis y programacion estadstica R y en contreto algunos de los libros especficos de funciones, actualmente
disponibles en la web http://cran.es.r-project.org, para la aplicacion practica
de los metodos de regresion no parametricos. Nuestra atencion se ha centrado fundamentalmente en los libros kernSmooth, locpol, locfit y sm, si bien existen
funciones disponibles en otros libros (stats, monreg, lokern, loess, lowess, np,
psplines, etc.). Asociado a alguna aplicacion de datos concreta se ilustrara tambien
el uso de alguna de estas funciones adicionales.
De estos libros de R se ha hecho una descripcion casi exhaustiva, obviando solo
en algunos casos partes que no corresponden a los objetivos concretos de este trabajo
y mas concretamente de los metodos teoricos desarrollados en el captulo primero.
Hemos de destacar que todo el trabajo aqu desarrollado esta sujeto a la necesaria
y continua actualizacion, dado el rapido avance en esta materia computacional. En
este sentido para cada libro se ha especificado la version utilizada incorporando
dentro de la u
ltima (hasta septiembre de 2009) version de R.
Para finalizar este trabajo se desarrollan en el captulo tres algunas aplicaciones
practicas. Nos hemos centrado en modelos de regresion univariante, haciendo uso de
algunas de las funciones analizadas en el captulo dos. Se han ilustrado los metodos
de regresion no parametrica para distintos conjuntos de datos. Hemos querido realizar dichas ilustraciones usando datos reales (disponibles en libros de R) y tambien
mediante ejercicios de simulacion. Con esto pretendemos iniciarnos en el metodo
habitualmente usado en la investigacion para la validacion practica de las nuevas
metodologas propuestas. En la resolucion de estos ejercicios practicos hemos utilizado varias funciones disponibles para el mismo problema. Esto nos ha permitido

INDICE GENERAL
realizar conclusiones acerca de dichos procedimientos.

A modo de conclusion podemos decir que este trabajo nos ha permitido un
acercamiento a los metodos de regresion no parametrica mas habituales con un
enfoque eminentemente practico. Las ilustraciones desarrolladas con datos permiten
explorar el modo en que actualmente se puede trabajar en R, para dar soluciones
al problema de regresion no parametrico. Los metodos de regresion polinomial local
ofrecen una sencilla va de solucion, suficientemente documentada e implementada
en diversos libros de R. No obstante sigue siendo necesaria la incorporacion de nuevas
funciones, que permitan la implementacion de novedosas tecnicas, que ocupan las
publicaciones actuales en dicha materia. Estamos seguros de que es cuestion de
no mucho tiempo y entre nuestros objetivos esta el de intentar colaborar en dicha
tarea. El inicio de este proposito se traduce en el trabajo que actualmente tiene en
sus manos.
Nisa Boukichou Abdelkader
Granada, 27 de septiembre de 2009
INDICE GENERAL
Captulo 1
Introducci
on
1.1.
Antecedentes
Consideremos en primer lugar el planteamiento de un problema de regresion,

habitualmente especificado como sigue:
Sea un conjunto de n observaciones, {(Xi , Yi ), i = 1, . . . , n}, de una variable
aleatoria bidimensional, (X, Y ), satisfaciendo el modelo,
Yi = m(Xi ) + i ,
i = 1, . . . , n,
(1.1)
donde
los residuos i son variables aleatorias independientes (vv.aa.ii.) con media
cero y varianza 2 (Xi ),
y la funcion m es desconocida y se define como la funcion de regresion, m(x) =
E[Y |X = x].
Un planteamiento tal corresponde a un modelo de regresion de tipo univariante,
esto es, con tan solo una variable explicativa, basado en un dise
no aleatorio, donde
las observaciones constituyen una muestra aleatoria de la poblacion (X,Y). Ademas,
estamos considerando una situacion general de heterocedasticidad, es decir, las
varianzas de los errores se suponen distintas.
Con tal planteamiento el interes se centra principalmente en los tres objetivos
siguientes:
1. Explorar y representar la relacion existente entre la variable explicativa, X, y
la variable Y .
2. Predecir el comportamiento de Y para valores de la variable X a
un no observados.
5
CAP
ITULO 1. INTRODUCCION
3. Estimar las derivadas de m, as como algunas propiedades relacionadas.

Para alcanzar tales objetivos se puede optar por dos aproximaciones diferentes:
asumir alg
un modelo parametrico para la funcion de regresion o no imponer tal restriccion, asumiendo tan solo que la funcion presenta algunas propiedades deseables
generalmente relacionadas con la derivabilidad.
La primera aproximacion correspondera a lo que se denomina Regresi
on Param
etrica, y supone que la funcion de regresion desconocida, m, pertenece a alguna
familia parametrica de funciones, m {m | }. De esta manera se procede a la
seleccion de alguna de ellas seg
un alg
un criterio previamente especificado, como
puede ser el criterio de mnimos cuadrados. Por contraposicion a la segunda
aproximacion se le denomina Regresi
on No Param
etrica, y la u
nica restriccion
que se le impone a la funcion m, es que sea suave, entendiendo esta suavidad en
terminos de derivabilidad. En este u
ltimo sentido a estos metodos tambien se les
suele denominar metodos de suavizamiento o suavizado.
De partida podemos destacar de las dos aproximaciones una serie de ventajas e
inconvenientes:
- Los modelos de Regresi
on Parametrica en general son sencillos y de calculo
rapido. Ademas, el presuponer que tienen una forma definida y parametrica
para la funcion de regresion, hace que se tengan garantizadas de antemano
las propiedades de las estimaciones resultantes. Sin embargo, se caracterizan
por ser metodos poco flexibles y de difcil adaptacion en diversas situaciones
reales.
- El planteamiento no parametrico por el contrario, permite una mayor flexibilidad y por eso, se considera una de las mejores herramientas de tipo exploratorio. No obstante, el no realizar muchas hipotesis sobre la estructura subyacente
en los datos se traduce en soluciones mucho mas complejas y con un mayor
coste computacional.
En este trabajo nos centraremos en la teora de Regresi
on No Parametrica y de
hecho sera la aproximacion que desarrollaremos con mas detalle en lo que sigue.
1.2.
Estimaci
on del modelo de Regresi
on No Param
etrica
La teora y los metodos de suavizamiento o regresion no parametrica han cobrado

un gran auge en las u
ltimas decadas unido al avance en materia computacional. Una
revison de los mismos se puede encontrar por ejemplo en los libros de Wand y
Jones (1995), Fan y Gijbels (1996) y Loader (1999). El creciente interes por estas
DEL MODELO DE REGRESION

NO PARAMETRICA
1.2. ESTIMACION
metodologas ha tenido dos razones principales: la primera, que los estadsticos se

dieron cuenta de que planteamientos puramente parametricos no aportaban la flexibilidad necesaria para la estimacion de las curvas que aparecan en la practica. Y
la segunda razon estaba ligada al avance de la informatica y al desarrollo de un
hardware que posibilitaba el costoso calculo de esos estimadores no parametricos.
Los primeros estimadores de regresion no parametrica propuestos fueron los sencillos
estimadores de tipo n
ucleo de Nadaraya (1964) y Watson (1964). Dichos estimadores
se han ido refinando y perfeccionando dentro de los denominados metodos de regresion polinomial local, convirtendose en uno de los metodos mas empleados por
diversos analistas en la actualidad.
Y centrando la atencion en la adecuada eleccion de los parametros decisivos en
el buen comportamiento de los estimadores de regresion resultantes, como son: el
ancho de banda
o parametro de suavizado y el grado de los ajustes polinomiales
locales. A continuacion, en Eubank (1988) se puede encontrar el siguiente postulado
que dice lo siguiente:
Supongamos que la funcion de regresi
on desconocida m, es suave entonces, podemos esperar que las observaciones tomadas en puntos proximos a uno dado x puedan
darnos informacion del vector de dicha funcion en x.
Con este planteamiento y sobre el diagrama de dispersion de los datos, se trata
de definir unas bandas centradas en cada punto y de un ancho h, y calcular la estimacion en el punto utilizando tan solo las observaciones que caen dentro. Ademas,
para obtener la curva de regresion estimada, la banda definida por el parametro
h recorrera todo el diagrama de dispersion de izquierda a derecha. No obstante
tambien, al parametro h se le suele llamar ancho de banda o parametro de suavizado.
Figura 1.1: Ancho de Banda
CAP
Como podemos observar en este grafico (Figura 1.1), en primer lugar hemos
definido unas bandas centradas en los puntos x = 25 y x = 30 con un ancho h en
cada una de ellos. Y despues, hemos estimado las curvas de regresion en cada uno
de los puntos utilizando solamente las observaciones que caen dentro de cada una
de las bandas.
Por tanto, siguiendo esta idea se han desarrollado diversas tecnicas, entre las que
destacaremos estos cuatro tipos:
1. Estimadores tipo n
ucleo: realizan un promedio de las observaciones que
caen en cada banda. Se definen como:
c(x) =
m
n
X
Wi (x)Yi
(N adaraya W atson, 1964)
i=1
donde:

Xi x
Wi (x)
= n1 h1 k
h
2. Regresi
on Polinomial Local: realiza un ajuste polinomial con las observaciones que caen en la banda. Se define como:
mn
n
X
i=1
Y
i
p
X
j=0
j (Xi x)j h1 k
Xi x
h
(Cleveland, 1979)
3. Suavizamiento por Splines: se define como la solucion a un problema de

mnimos cuadrados penalizados. Y se calcula de la siguiente forma:
min n1
n
X
i=1
{m(Xi ) Yi }2 +
(m(r) (x))2 dx ( > 0)
4. Estimadores basados en Desarrollos en Serie Ortogonal o los que estan

tan de moda denominados Wavelets, se definen como:
DEL MODELO DE REGRESION

NO PARAMETRICA
1.2. ESTIMACION
cN (x) =
m
N
X
j=1
bj qj (x)
donde:
bj = n1
n
X
qj (Xi )Yi
i=1
Ruppert, Wand y Carroll (2003) establecen los siguientes factores a tener en

cuenta a la hora de evaluar y decidir que tipo de estimador no parametrico utilizar
en la practica:
1. Conveniencia. Esta disponible en el software de uso habitual o favorito por
el analista?
2. Facilidad para su implementacion. Si no esta disponible directamente, es facil
implementarlo en el lenguaje de programacion habitual del analista?
3. Flexibilidad. Se trata de ver si el estimador es capaz de explicar un amplio
abanico de tipos de relaciones que pueden existir entre las variables de interes.
4. Simplicidad y sencillez. Es intuitivo? esto es, si es facil entender como el
metodo act
ua sobre los datos para dar respuestas.
5. Tratabilidad: es sencillo estudiar las propiedades matematicas del estimador?.
6. Fiabilidad. Si se tiene garanta de que las respuestas proporcionadas por el
estimador son verdad.
7. Eficiencia. Si el estimador hace un uso de los datos eficiente.
8. Posibilidades de extension a otras situaciones o problemas mas complicados.
Este trabajo se desarrolla desde una perspectiva computacional y practica, en
este sentido ponemos especial atencion en los puntos 1, 2, 3 y 8 anteriores. As,
desde tales premisas, los metodos de regresion polinomial local constituyen una
adecuada eleccion. De hecho se trata de un procedimiento sencillo y muy intuitivo,
que esta implementado y ampliamente documentado en programas estadsticos y
lenguajes de programacion de uso extendido. Destacan implementaciones en S-Plus,
MatLab y en R atraves de diversos libros especficos.
A continuacion describiremos con mas detalle el procedimiento de regresion polinomial local desde el punto de vista teorico.
10
CAP
1.3.
Regresi
on Polinomial Local
Si suponemos que la funcion de regresion,m tiene p derivadas en un punto x0 ,

entonces va el teorema de Taylor tenemos una aproximacion de este tipo para los
valores en un entorno de x0 .
m(x) m(x0 ) + m (x0 )(x x0 ) +
m (x0 )
m(p) (x0 )
(x x0 )2 + ... +
(x x0 )p
2!
p!
Luego, esto justifica que se puede aproximar localmente m por funciones polin
omicas de grado p.
Pp (x) =
p
X
j (x x0 )j
j=0
As, se obtienen estimaciones de los coeficientes bj con j = 0, . . . , p y entonces,

observando la expresion, vemos que la estimacion del termino independiente 0
sera un estimador de m en x0 y el resto de coeficientes bj proporcionaran estimaciones de sus derivadas.
Por eso, con el fin de estimar m localmente mediante polinomios de grado p
consideraremos un problema de mnimos cuadrados ponderados:
min
n
X
Yi
i=1
donde:
p
X
j=0
j (Xi x0 )j kh (Xi x0 )
h es un parametro denominado ancho de banda o par

ametro de suavizado
que controla las observaciones que caen en cada entorno.
Kh (u) = h1 K( uh ), donde la funcion K(), se denomina funci
on n
ucleo.
Dicha funcion define las ponderaciones que se asignan a cada observacion en
el entorno local considerado. Habitualmente se supone una densidad simetrica
y con soporte compacto.
y p es el grado del ajuste polinomial local.
Ademas,como casos particulares se puede obtener el conocido estimador n
ucleo
de Nadaraya-Watson, que supone realizar ajustes polinomiales locales de grado cero,
y tambien cuando el ajuste polinomial es de grado uno, se obtiene el denominado estimador lineal local. Si bien los ajustes constantes has sido estudiados y ampliamente
11
POLINOMIAL LOCAL
1.3. REGRESION
utilizados por teoricos y analistas de datos, es el ajuste lineal el que ha mostrado

ser mas conveniente y el mas usado actualmente en la practica (en este sentido es
interesante la discusion del artculo de Cleveland (1979).
Como anteriormente hemos mencionado la definicion del estimador esta determinado por tres parametros: el ancho de banda, h, la funcion n
ucleo y el grado
p.
El ancho de banda se define como un parametro positivo cuyo rango en principio sera la amplitud del intervalo de estimacion. No obstante la eleccion de dicho
parametro constituye uno de los aspectos cruciales del procedimiento de estimacion.
Las buenas propiedades del estimador resultante dependera en gran medida de la
eleccion que se haga de dicho parametro. En este sentido existe una compensacion
(trade-off en ingles) entre flexiblidad y complejidad que es contralada a traves de dicho parametro. As, si dicho parametro se toma muy peque
no, tan solo observaciones
muy proximas al punto de estimacion intervendran en el calculo del estimador, describiendo muy bien comportamientos locales pero obtendremos una curva estimada
muy variable. Si por el contrario, dicho parametro se toma muy grande, las estimaciones en cada punto se veran afectadas por observaciones en puntos alejados de
forma que difcilmente se podran recoger los comportamientos locales, dando lugar
a grandes sesgos y como consecuencia obtendremos poca variabilidad. Notese por
tanto que una practica con anchos de banda muy grandes se traduce en procedimientos computacionalmente costosos.
As, desde un punto de vista teorico, la eleccion del ancho de banda tendra que
buscar una adecuada compensacion entre sesgo y varianza. Y desde un punto de
vista practico supondra la eleccion de la complejidad del modelo.
Por otro lado, con la eleccion del grado de los ajustes polinomiales nos ocurre
algo similar. Es decir, la utilizacion de ajustes de grado cero o uno nos daran estimaciones con poca variabilidad, muy suaves, pero con sesgos muy elevados. Por
el contrario, cuando los ajustes son con grados mayores (dos o tres) nos permitiran
mayor adaptabilidad, o sea, menores sesgos pero, obtendremos mayor varianza.
Luego, del mismo modo que para la eleccion del ancho de banda es necesario una
compensacion entre sesgo y varianza, para elegir p se buscara tambien la compensacion optima entre sesgo y varianza.
Y finalmente, la funcion n
ucleo, K, define la forma de los pesos que se asocian a
cada observacion dentro de la banda definida y por tanto, determina su importancia
en el calculo de la estimacion. Para dicha funcion es habitual el uso de alguna de las
siguientes densidades:
1. Triangular: k(u) = (1 |u|)1|u|1
2. Epanechnikov: k(u) = 43 (1 u2 )1|u|1
3. Biponderado: k(u) =
15
(1
16
u2 )2 1|u|1
12
CAP
4. Gaussiano: k(u) = (2)1 exp( x2 )1|u|1

En la practica, la eleccion del n
ucleo no afecta al buen comportamiento de las
estimaciones resultantes, por lo que considerar una funcion u otra atiende fundamentalmente a razones tecnicas.
Por otro lado, tenemos que el ancho de banda y el grado del polinomio necesitan
una compensacion entre el sesgo y la varianza, por tanto, se definen a continuacion
criterios para medir dicho compromiso y asimismo, la consecuente bondad de las
estimaciones.
Luego, entre las diversas medidas que se pueden definir, hemos considerado las
siguientes:
- El Error Cuadratico Medio Condicional en su version local (MSE):
h
dh (x0 )) = EY /X (m
dh (x0 ) m(x0 ))2
MSE(m
- El Error Cuadratico Medio Condicional en su version global (MISE):

dh ) = EY /X
MISE(m
Z
dh (x) m(x))2 dx
(m
Por tanto, como sabemos el error cuadratico medio se puede descomponer en

un termino de sesgo al cuadrado mas uno de varianza, por lo que podemos asumir
que este criterio permitira la compensacion entre sesgo y varianza anteriormente
mencionado.
n
o2
dh (x0 )) = EY /X [m
dh (x0 )] m(x0 )
M SE(m
|
{z
sesgo2
dh (x0 ))
+ V arY /X (m
|
{z
varianza
En cuanto, a la eleccion entre un criterio u otro dependera de que busquemos

optimalidad en cada punto (MSE) o bien, que busquemos la optimalidad conjunta
en todo el intervalo de estimacion (MISE).
A continuacion, veremos varios procedimientos de estimacion de h o tambien
denominados selectores del ancho de banda o metodos se seleccion de la complejidad
del modelo.

DE LA COMPLEJIDAD DEL MODELO
1.4. METODOS
DE SELECCION
1.4.
13
M
etodos de selecci
on de la complejidad del
modelo
A) M
etodos Plug-in
La idea que subyace en este tipo de tecnicas es que a partir de la expresion
optima teorica del ancho de banda, se proponen estimadores de lo desconocido
que se incrustaran en dicha expresion teorica. No obstante, aunque estos procedimientos comparten una misma filosofa, podemos hacer una clasificacion
de los mismos en dos categoras:
- Plug-in Asint
otico:
Parte de expresiones asintoticas teoricas del error (MSE o MISE) y estima
lo desconocido. Y uno de los inconvenientes que tiene es que es necesario
realizar hipotesis bastante restrictivas sobre la funcion de regresion desconocida, lo cual limita su campo de aplicacion en la practica. Ademas,
es valida para tama
nos muestrales grandes porque en caso, de tama
nos
muestrales peque
nos dicho procedimiento no es efectivo.
- Plug-in no Asint
otico:
Debido al problema que presenta el plug-in asintotico en muestras peque
nas,
se toma como punto de partida expresiones desarrolladas de los errores
teoricos de naturaleza no asintotica, incrustando en ellas las estimaciones
propuestas para los terminos desconocidos. Por otro lado, dependiendo
de si utilizan iteraciones para llegar a la soluci
on
optima, se pueden
distinguir entre:
- Plug-in iterativo:
A
naden una sucesion de iteraciones con el fin de conseguir una progresiva
mejora del ancho de banda. Su inconveniente es que aumenta el coste
computacional.
- Plug-in no iterativo
o directo:
Son los que no utilizan iteraciones para la obtencion del selector del
parametro.
B) Otros m
etodos de selecci
on autom
atica basados en los datos:
Se basan en la minimizacion de criterios de error relacionados con con los
anteriores MSE y MISE pero conocidos, como es la suma residual de cuadrados
para obtener la seleccion del parametro. Ademas, dentro del mismo destacamos
algunas de las elecciones mas habituales, como son:
14
CAP
La validacion cruzada por mnimos cuadrados.

La validacion cruzada generalizada.
El criterio de informacion de Akaike.
El error de prediccion finito.
El selector de Shibata.
La T de Rice.
C) Procedimientos basados en aproximaciones Bootstrap
La metodologa bootstrap tiene como proposito ganar informacion acerca de
la distribucion de un estimador. Sin embargo, en regresion no parametrica
la metodologa bootstrap es utilizada fundamentalmente para dos tareas: la
primera es la de elegir el parametro de suavizamiento o ancho de banda y la
segunda es la de construir intervalos de confianza para la curva de regresion.
1.5.
Extensi
on multivariante
Consideramos ahora la extension del modelo (1.1) al caso en que se considera

mas de una covariable. De este modo sea Xi un vector de D covariables y el modelo
general de regresion multivariante heterocedastico dado por:
Yi = m(Xi ) + i
(1.2)
donde la funcion (x) = V ar[Y |X = x] es finita, y los residuos, i , son variables independientes e identicamente distribuidas, con media cero, varianza 2 (Xi )
y son independientes de los vectores aleatorios, Xi . La funcion de regresion, m(x) =
E[Y |X = x], con x = (x1 , . . . , xD )T .
El modelo de regresion lineal multivariante supone que la relacion entre la variable de respuesta Y y cada una de las variables independientes es lineal. A veces,
es evidente que esta relacion no es lineal, por lo que hay que considerar modelos
que sean mas flexibles. Las tecnicas de regresion no parametricas responden a esta
flexibilidad ya que no imponen condiciones sobre la forma de la funcion D-variante,
m(x).
En esta situacion, los estimadores mas comunes para m(x) son versiones multivariantes de los estimadores tipo n
ucleo (como los polinomiales locales descritos
anteriormente) o splines de suavizamiento. Ruppert y Wand (1994) introducen la
extension multivariante del estimador polinomial local. A continuacion describimos
el estimador polinomial local para el caso general de grado p y D > 1. Consideramos
el siguiente problema de mnimos cuadrados ponderados:
15
MULTIVARIANTE
1.5. EXTENSION
mn
n
X
Yi
i=1
p
X
l1 ,...,ld
D
Y
j=1
L=1 l1 +...+ld =L
(Xi xj )lj KH (Xi x)
donde: = {l1 ,...,lD : l1 + . . . + lD = L} y L = {0, . . . , p} es un vector de coeficientes, H es una matriz de dimension D D simetrica, definida positiva; K() es
una funcion n
ucleo no negativa D-variante y KH (u) = |H|1/2 K(H 1/2 u).
A la matriz H se le denomina matriz ancho de banda, dado que es la extension
multivariante del parametro ancho de banda univariante. Si denotamos por bj , j =
0, . . . , p, a las soluciones del problema anterior entonces, usando el desarrollo de
Taylor, el estimador polinomico local vendra dado por la primera de ellas, esto es:
cp (x) = b0,...,0
m
El problema en forma matricial y su solucion, se definen mediante:

mn (Y Xx )T Wx,H (Y Xx )
cp (x) = eT1 XTx Wx,H Xx

m
1
XTx Wx,H Y
donde e1 es un vector (pD + 1), con un 1 en la primera posicion y 0 en el resto,

Y = (Y1 , . . . , Yn )T ,
1 (X1 x)T ((X1 x)p )T
.
..
..
...
Xx =
.
.
..
1 (Xn x)T ((Xn x)p )T
y
Wx,H = diag {KH (X1 x), . . . , KH (Xn x)}
Casos particulares:
p = 0, este coincide con la version multivariante del estimador n
ucleo de
Nadaraya-Watson:
Pn
KH (Xi x)Yi
i=1 KH (Xi x)
i=1
c0 (x) = P
m
n
16
CAP
p = 1, es el estimador lineal local multivariante (Ruppert y Wand, 1994):

c1 (x) = eT1 XTx,H Wx,H Xx

m
1.5.1.
1
XTx Wx,H Y
El problema de la dimensionalidad
Aunque la generalizacion al caso multidimensional de la mayora de tecnicas

de suavizamiento es posible, aparece un problema importante conocido como: el
problema de la dimensionalidad, (en ingles the curse of dimensionality, Bellman,1961).
Este problema se refiere al hecho de que cuando estamos estimando, considerando
un entorno con un n
umero fijo de datos, y tenemos una superficie de gran dimension,
dicho entorno puede ser demasiado grande como para ser llamado local, es decir,
si un entorno local contiene 10 datos de cada variable, entonces el correspondiente
entorno local D dimensional contiene 10D datos.
Como consecuencia se necesitan conjuntos de datos mucho mas grandes incluso
cuando D no es muy elevado, y en la practica puede que tales conjuntos no esten
disponibles.
Otro problema que presentan los estimadores multivariantes es la falta de interpretabilidad, ya que sera difcil de visualizarlos graficamente. No es posible representar superficies para D > 2.
Tambien resulta un inconveniente el excesivo coste computacional de las versiones
multivariantes, que requieren un gran n
umero de operaciones. Esto hace que en la
practica los suavizadores multidimensionales solo se apliquen hasta dimensiones 2
o 3.
Los problemas anteriores nos llevan a plantear modelos alternativos que eviten
estos inconvenientes, como los que se comentan en este trabajo y que reciben el
nombre de modelos aditivos.
Los modelos aditivos se presentan como una herramienta u
til para el analisis de
datos. Estos modelos mantienen una importante caracterstica de interpretacion de
los modelos lineales, al tener representada cada variable de forma separada. As la
naturaleza de los efectos de una variable sobre la variable de respuesta no depende
de los valores de las otras variables.
Los modelos aditivos fueron formulados por Friedman y Stuetzle, (1981) y constituan el centro del algoritmo ACE de Breiman y Friedman(1985).
1.5.2.
Modelos aditivos no param

etricos
Propuestos por Hastie y Tibshirani (1990), son una tecnica de regresion no

parametrica multivariante muy utilizada. Frente a otras tecnicas de regresion no
17
DEL PARAMETRO
1.6. SELECCION
ANCHO DE BANDA
parametrica multivariante, los modelos aditivos pueden interpretar facilmente los

efectos individuales de las variables independientes sobre la variable dependiente,
sin tener en cuenta que el n
umero de variables sea elevado. Si partimos del modelo
de regresion m
ultiple heterocedastico (1.2), el modelo aditivo supone que la funcion,
m, se puede escribir como la suma de funciones univariantes de cada una de las
variables independientes, que se suponen suaves; es decir:
h
m(x) = E Y |X = (x1 , . . . , xD )T = + m1 (x1 ) + . . . + mD (xD )

Para estimar este modelo se han propuesto en los u
ltimos a
nos diferentes metodos, entre ellos destacamos el algoritmo backfitting, el metodo de integraci
on marginal
y el smooth backfitting (Nielsen and Sperlich, 2005), como los mas importantes, y a
partir de estos se han propuesto distintas modificaciones pero, todas ellas con el fin
de mejorar las estimaciones obtenidas.
El algoritmo backfitting, propuesto por Buja, Hastie y Tibshirani (1989) proporciona un metodo iterativo para estimar las componentes unidimensionales. Sin
embargo, el hecho de ser un proceso iterativo dificulta el estudio de las propiedades
asintoticas de las estimaciones. Y por eso, esta dificultad del estudio de las propiedades
del algoritmo backfitting es lo que lleva a otros autores a plantearse metodos de estimacion alternativos.
De esta forma, Linton y Nielsen(1995) proponen el m
etodo de integraci
on
marginal para un modelo aditivo bivariante, este es un metodo que estima directamente cada una de las componentes univariantes y permite el estudio de sus
propiedades, y fue generalizado para un modelo multivariante (D > 2) por Hengartner (1996) y Kim, Linton y Hengartner (1997).
1.6.
Selecci
on del par
ametro ancho de banda
La seleccion de una matriz de anchos de banda en el caso general multivariante,

H, a partir de los datos, ha sido mucho menos considerado en la literatura que el
problema en el caso univariante, en el que interviene un u
nico ancho de banda.
En el caso de los modelos aditivos, el parametro es un vector que define el grado
de suavizamiento empleado en cada componente del modelo y puede considerar un
ancho de banda diferente para cada una de las componentes aditivas o un u
nico
parametro com
un para todas ellas.
Al igual que en el caso univariante se pueden utilizar diversos metodos de seleccion del parametro ancho de banda. Los mas destacados son el metodo plug-in y
el metodo de validacion cruzada. Para una revision de estos metodos se pueden ver
en el libro de Hastie y Tibshirani (1990) y en el de Martnez Miranda et al., (2008).
18
CAP
El metodo de validacion cruzada de este modelo se puede encontrar en el libro

de Kim, Linton y Hengartner, 1999; en el de Kauermann y Opsomer, 2003 y en el
de Nielsen y Sperlich, 2005.
Y para los metodos plug-in ver los libros de Opsomer y Ruppert, 1997; el de
Severance-Lossin y Sperlich, 1997 y el de Mammen y Park, 2005.
Captulo 2
Software disponible en R
2.1.
Introducci
on
El programa R es un entorno de analisis y programacion estadstico que forma

parte del proyecto de software libre GNU General Public Licence. R esta disponible
en la direccion http://www.r-project.org. El proyecto R comenzo en 1995 por
un grupo de estadsticos de la universidad de Auckland, dirigidos por Ross Ihaka
y Robert Gentleman. R esta basado en el lenguaje de programacion S, dise
nado
especficamente para la programacion de tareas estadsticas en los a
nos 80 por los
Laboratorios Bell AT&T. El lenguaje S se considera un lenguaje de programacion
estadstica orientado a objetos de alto nivel.
Frente a otros lenguajes de programacion, R es sencillo, intuitivo y eficiente
ya que se trata de un lenguaje interpretado (a diferencia de otros como Fortran,
C++, Visual Basic, etc.). Como programa de analisis estadstico, R-base permite
realizar tareas estadsticas sencillas habituales y ademas permite extensiones que
implementan tecnicas estadsticas avanzadas. De este modo se cubre las necesidades
de cualquier analista, tanto en el ambito de la estadstica profesional como en el de
la investigacion estadstica.
R consta de un sistema base pero la mayora de las funciones estadsticas vienen
agrupadas en distintos libros (o bibliotecas del ingles packages) que se incorporan
de forma opcional. Para los metodos de regresion no parametrica existen funciones
disponibles en el libro basico stats, no obstante el uso mas adecuado de dichos metodos puede conseguirse a traves de funciones incorporadas en varios libros adicionales
y actualmente disponibles en la web. Entre estos libros destacan kernSmooth,
locpol, np, locfit,lokern, monreg, loess, sm, lowess etc.
En este captulo haremos una descripcion casi exhaustiva de algunos de estos
libros. En concreto nos centraremos en los libros kernSmooth, locpol, locfit y sm.
Se ha obviado solo en algunos casos, partes que no corresponden a los objetivos
concretos de este trabajo y mas concretamente de los metodos teoricos desarrollados
19
20
CAP
ITULO 2. SOFTWARE DISPONIBLE EN R
en el captulo primero.
Hemos de destacar que todo el trabajo aqui desarrollado esta sujeto a la necesaria
y continua actualizacion, dado el rapido avance en esta materia computacional. En
este sentido los libros descritos se ha actualizado solo en la fecha que se indica en
sus correspondientes apartados.
2.2.
Libro KernSmooth
En este libro llamado KernSmooth se trabajara con funciones para el suavizado

del n
ucleo y la estimacion de densidades correspondientes al texto de Wand y Jones
(1995). Dicho libro corresponde a la actualizacion del 24-05-2009, version 2.23-2 y
esta sujeta a cualquier modificacion a da de hoy porque esta en constante actualizacion.
Ademas, nos vamos a encontrar en este libro diferentes funciones con las que
podremos trabajar el tema que nos ocupa como es la regresion no parametrica. Los
detalles de cada una de las funciones se veran en cada seccion de las mismas con
mas detenimiento pero, basicamente tratan de calcular el parametro de seleccion, el
estimador de densidad, el calculo de los bins y la funcion de regresion.
Funci
on bkde
Descripci
on: devuelve las coordenadas x e y de un Binned estimado de la
densidad del n
ucleo para la funcion de probabilidad de los datos.
Funci
on:
bkde(x, kernel = "normal", canonical = FALSE, bandwidth,
gridsize = 401L, range.x, truncate = TRUE)
donde:
El n
ucleo kernel puede ser: normal, box (rectangular), epanech (Beta(2,2)),
biweight (Beta(3,3)), triweight (Beta(3,3)).
Valor devuelto por la funci
on: esta funcion devuelve un grafico de la funcion
de probabilidad.
Ejemplo de uso:
> data(geyser, package="MASS")
> x <- geyser$duration
21
2.2. LIBRO KERNSMOOTH
0.3
0.0
0.1
0.2
est$y
0.4
0.5
> est <- bkde(x, bandwidth=0.25)

> plot(est, type="l")
est$x
Figura 2.1: KernSmooth de bkde
Funci
on bkde2D
Descripci
on: devuelve la red de puntos y la matriz de densidad estimada inducida. El n
ucleo es la densidad normal bivariada.
Funci
on:
bkde2D(x, bandwidth, gridsize = c(51L, 51L), range.x, truncate = TRUE)
on: esta funcion devuelve dos graficos, el primero
es el contorno de las variables x1 y x2 y las estimaciones de fb. Y el segundo grafico,
es la perpestiva de las estimaciones de fb.
Ejemplo de uso:
> data(geyser, package="MASS")
> x <- cbind(geyser$duration, geyser$waiting)
22
CAP
100
120
> est <- bkde2D(x, bandwidth=c(0.7, 7))

> contour(est$x1, est$x2, est$fhat)
> persp(est$fhat)
0.001
0.003
0.007
80
0.006
0.006
0.004
60
0.002
0.007
40
0.005
Figura 2.2: KernSmooth de bkde2D1
Z
est$fhat
Figura 2.3: KernSmooth de bkde2D2
Funci
on bkfe
Descripci
on: calcula una estimacion aproximada (binned ) para la estimacion
tipo n
ucleo de una funcion de densidad. El n
ucleo es la densidad normal estandar.
23
Funci
on:
bkfe(x, drv, bandwidth, gridsize = 401L, range.x, binned = FALSE,
truncate = TRUE)
donde:
El drv: es el orden de la derivada en la funcion de densidad.
on: esta funcion lo que nos devuelve es un valor
aproximado de la estimacion de la densidad.
Ejemplo de uso:
>
>
>
>
data(geyser, package="MASS")
x <- geyser$duration
est <- bkfe(x, drv=4, bandwidth=0.3)
est
[1] 33.04156
Funci
on dpih
Descripci
on: Se utiliza directamente el metodo plug-in para seleccionar el tama
no
del bin de un estimador de tipo histograma.
NOTA: Plug-in puede ser asintotico y no asintotico, y si se usa iteraciones para
llegar a la solucion optima pueden ser plug-in iterativo y plug-in no iterativo (directo).
Funci
on:
dpih(x, scalest = "minim", level = 2L, gridsize = 401L,
range.x = range(x), truncate = TRUE)
donde:
El minim es el mnimo de la desviacion tpica estandar (stdev ) y del rango
intercuartilico (iqr ) dividido por 1349.
24
CAP

on: esta funcion devuelve un valor y con este se
construye el vector de los bins que posteriormente se utilizaran para construir el
histograma.
Ejemplo de uso:
>
>
>
>
>
h <- dpih(x)
bins <- seq(min(x)-0.1, max(x)+0.1+h, by=h)
hist(x, breaks=bins)
> h <- dpih(x)

> h
[1] 0.2088382
> bins <- seq(min(x)-0.1, max(x)+0.1+h, by=h)
> bins
0.9421715
2.4040386
3.8659057
5.3277728
1.1510096
2.6128767
4.0747438
5.5366110
1.3598478 1.5686859 1.7775241 1.9863623

2.8217149 3.0305530 3.2393912 3.4482294
4.2835820 4.4924202 4.7012583 4.9100965
5.7454491
30
20
Frequency
40
50
60
Histogram of x
10
0.7333333
2.1952004
3.6570675
5.1189346
[1]
[8]
[15]
[22]
Figura 2.4: KernSmooth de dpih
25
Funci
on dpik
Descripci
on: Se utiliza directamente el metodo plug-in para seleccionar el ancho
de banda de la densidad estimada tipo n
ucleo.
Funci
on:
dpik(x, scalest = "minim", level = 2L, kernel = "normal",
canonical = FALSE, gridsize = 401L, range.x = range(x),
truncate = TRUE)
donde:
El minim es el mnimo de la desviacion tpica estandar, stdev y del rango
intercuartlico, iqr dividido por 1.349.
on: esta funcion devuelve el valor del ancho de
banda y despues, lo utiliza con la funcion bkde para calcular los bins estimados de la
densidad del n
ucleo para finalmente representarlos graficamente dicha estimacion.
Ejemplo de uso:
>
>
>
>
h <- dpik(x)
h
[1] 0.1434183
> est <- bkde(x, bandwidth=h)
> plot(est,type="l")
26
0.4
0.0
0.2
est$y
0.6
CAP
est$x
Figura 2.5: KernSmooth de dpik
Funci
on dpill
Descripci
on: Se utiliza el metodo plug-in para la seleccion del ancho de banda
de un estimador lineal local con n
ucleo gaussiano, descrito por Ruppert, Sheather y
Wand (1995).
Funci
on:
dpill(x, y, blockmax = 5, divisor = 20, trim = 0.01, proptrun = 0.05,
gridsize = 401L, range.x, truncate = TRUE)
donde:
El blockmax es el n
umero maximo de bloques de los datos para la construccion
de un parametro inicial estimado.
El trim es la proporcion de la muestra de cada uno de los extremos recortados
en la direccion x antes de aplicar el metodo plug-in.
El proptrun es la proporcion de la serie de x en cada extremo truncado en la
funcion estimada.
27

on: esta funcion devuelve el valor del ancho de
banda seg
un el metodo descrito por Ruppert, Sheather y Wand (1995). Este valor
se puede utilizar para calcular el estimador de regresion lineal local.
Ejemplo de uso:
> local({pkg <- select.list(sort(.packages(all.available = TRUE)))
+ if(nchar(pkg)) library(pkg, character.only=TRUE)})
KernSmooth 2.23 loaded
Copyright M. P. Wand 1997-2009
>
>
>
>
>
>
data(geyser, package = "MASS")

y <- geyser$waiting
plot(x, y)
h <- dpill(x, y)
h
[1] 0.2342897
50
60
70
80
90
100
110
> fit <- locpoly(x, y, bandwidth = h)

> lines(fit)
Figura 2.6: KernSmooth de dpill
28
CAP
Funci
on locpoly
Descripci
on: Estima una funcion de densidad de probabilidad, una funcion de
regresion o sus derivadas utilizando polinomios locales.
Funci
on:
locpoly(x, y, drv = 0L, degree, kernel = "normal",
bandwidth, gridsize = 401L, bwdisc = 25,
range.x, binned = FALSE, truncate = TRUE)
donde:
El drv es el orden de la derivada para ser estimada.
El bwdisc es el n
umero de anchos de banda equiespaciados logaritmicamente
utilizados para acelerar el calculo.

on: esta funcion devuelve un objeto con los resultados de la estimacion de la funcion de densidad o la funcion de regresion o sus
derivadas utilizando polinomios locales. A partir de este objeto se puede representar
el ajuste realizado usando el metodo plot.
Ejemplo de uso:
>
>
>
>
>
data(geyser, package = "MASS")

# local linear density estimate
est <- locpoly(x, bandwidth = 0.25)
plot(est, type = "l")
> # local linear regression estimate

> y <- geyser$waiting
> plot(x, y)
> fit <- locpoly(x, y, bandwidth = 0.25)
> lines(fit)
29
0.3
0.0
0.1
0.2
est$y
0.4
0.5
est$x
50
60
70
80
90
100
110
Figura 2.7: KernSmooth de locpoly
Figura 2.8: KernSmooth de locpoly
30
CAP
2.3.
Libro Locpol
En este libro se recogen funciones que calculan estimadores de tipo polinomial

local, seg
un se describen en Fan y Gijbels (1996) o Wand y Jones (1995). Dicho libro
corresponde a la actualizacion del 17-04-2009, version es del 0.2-0.
Dentro de este libro se recogen funciones para calcular cantidades u
tiles asociadas
a los n
ucleos, el estimador de densidad, la estimacion polinomial local as como
selectores del parametro ancho de banda de tipo plug-in y validacion cruzada.
Funci
on compKernVals
Descripci
on: permiten calcular valores u
tiles relacionados con el n
ucleo.
Funciones:
computeRK(kernel, lower=dom(kernel)[[1]], upper=dom(kernel)[[2]],
subdivisions = 25)
computeK4(kernel, lower=dom(kernel)[[1]], upper=dom(kernel)[[2]],
subdivisions = 25)
computeMu(i, kernel, lower=dom(kernel)[[1]], upper=dom(kernel)[[2]],
subdivisions = 25)
computeMu0(kernel, lower=dom(kernel)[[1]], upper=dom(kernel)[[2]],
subdivisions = 25)
Kconvol(kernel,lower=dom(kernel)[[1]],upper=dom(kernel)[[2]],
subdivisions = 25)
donde:
lower, upper: son los lmites de integracion inferior y superior.
Valor devuelto por las funciones: estas funciones devuelven un valor numerico
que son del tipo:
computeK4 : es el cuarto fin de autoconvolucion de K.
computeRK : es el segundo fin de autoconvolucion de K.
31
2.3. LIBRO LOCPOL
computeMu0 : es el integrante de K.
computeMu2 : es el segundo momento de orden K.
computeMu: es el i-esimo momento de orden K.
Kconvol : es la autoconvolucion de K.
Ejemplo de uso:
>
+
+
+
+
+
+
+
+
+
+
+
>
g <- function(kernels)
{
mu0 <- sapply(kernels,function(x) computeMu0(x,))
mu0.ok <- sapply(kernels,mu0K)
mu2 <- sapply(kernels,function(x) computeMu(2,x))
Rk.ok <- sapply(kernels,RK)
RK <- sapply(kernels,function(x) computeRK(x))
K4 <- sapply(kernels,function(x) computeK4(x))
res <- data.frame(mu0,mu0.ok,mu2,mu2.ok,RK,Rk.ok,K4)
res
}
g(kernels=c(EpaK,gaussK,TriweigK,TrianK))
1
2
3
4
mu0 mu0.ok
mu2
mu2.ok
RK
Rk.ok
1
1 0.2000000 0.200000 0.6000000 0.6000000
1
1 1.0000000 1.000000 0.2820948 0.2820948
1
1 0.1111111 0.111111 0.8158508 0.8158510
1
1 0.1666667 0.166667 0.6666667 0.6666670
K4
0.4337657
0.1994711
0.5879019
0.4793655
## Sacamos los valores de cada funci

on individualmente pero
## son los mismos que se observan en la tabla anterior.
> EpaK
function (x)
ifelse(abs(x) <= 1, 3/4 * (1 - x^2), 0)
attr(,"RK")
[1] 0.6
attr(,"RdK")
[1] 1.5
32
CAP
attr(,"mu0K")
[1] 1
attr(,"mu2K")
[1] 0.2
attr(,"K4")
[1] 0.4337657
attr(,"dom")
[1] -1 1
> gaussK
function (x)
dnorm(x, 0, 1)
attr(,"RK")
[1] 0.2820948
attr(,"RdK")
[1] 0.1410474
attr(,"mu0K")
[1] 1
attr(,"mu2K")
[1] 1
attr(,"K4")
[1] 0.1994711
attr(,"dom")
[1] -Inf Inf
> TriweigK
function (x)
ifelse(abs(x) <= 1, 35/32 * (1 - x^2)^3, 0)
attr(,"RK")
[1] 0.815851
attr(,"RdK")
[1] 3.18182
attr(,"mu0K")
[1] 1
attr(,"mu2K")
[1] 0.111111
attr(,"K4")
[1] 0.5879012
attr(,"dom")
2.3. LIBRO LOCPOL
[1] -1
33
> TrianK
function (x)
ifelse(abs(x) <= 1, (1 - abs(x)), 0)
attr(,"RK")
[1] 0.666667
attr(,"RdK")
[1] 2
attr(,"mu0K")
[1] 1
attr(,"mu2K")
[1] 0.166667
attr(,"K4")
[1] 0.4793729
attr(,"dom")
[1] -1 1
Funci
on denCVBwSelC
Descripci
on: calcula el selector mediante validacion cruzada del ancho de banda
para el estimador de la densidad de Parsen-Rosenblatt.
Funci
on:
denCVBwSelC(x, kernel = gaussK, weig = rep(1, length(x)),
interval = .lokestOptInt)
donde:
weig: es un vector de pesos para las observaciones.
on: esta funcion devuelve un valor numerico con
el ancho de banda.
Ejemplo de uso:
> stdy <- function(size=100,rVar=rnorm,dVar=dnorm,kernel=gaussK,x=NULL)
+ {
34
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
CAP
if( is.null(x) ) x <- rVar(size)

Tc <- system.time( dbwc <- denCVBwSelC(x,kernel) )[3]
ucvT <- system.time( ucvBw <- bw.ucv(x,lower=0.00001,upper=2.0) )[3]
nrdT <- system.time( nrdBw <- bw.nrd(x) )[3]
{
xeval <- seq( min(x)+dbwc , max(x)-dbwc ,length=50)
hist(x,probability=TRUE)
lines(xeval,trueDen <- dVar(xeval),col="red")
lines(density(x),col="cyan")
xevalDenc <- PRDenEstC(x,xeval,dbwc,kernel)
dencMSE <- mean( (trueDen-xevalDenc)^2 )
xevalucvDen <- PRDenEstC(x,xeval,ucvBw,kernel)
ucvMSE <- mean( (trueDen-xevalucvDen)^2 )
xevalDenNrd <- PRDenEstC(x,xeval,nrdBw,kernel)
nrdMSE <- mean( (trueDen-xevalDenNrd)^2 )
lines(xevalDenc,col="green")
lines(xevalucvDen,col="blue")
lines(xevalDenNrd,col="grey")
}
return( cbind( bwVal=c(evalC=dbwc,ucvBw=ucvBw,nrdBw=nrdBw),
mse=c(dencMSE,ucvMSE,nrdMSE),
time=c(Tc,ucvT,nrdT) ) )
}
> stdy(100,kernel=gaussK)
bwVal
mse time
evalC 0.1806189 0.5476062 6.92
ucvBw 0.1846771 0.5474982 0.06
nrdBw 0.3460004 0.5448560 0.01
> stdy(100,rVar=rexp,dVar=dexp,kernel=gaussK)
bwVal
mse time
evalC 0.1197885 4.646692 7.59
ucvBw 0.1209562 4.646678 0.00
nrdBw 0.3183121 4.647766 0.00
> stdy(200,rVar=rexp,dVar=dexp,kernel=gaussK)
35
2.3. LIBRO LOCPOL
bwVal
mse time
evalC 0.06563341 8.532414 33.04
ucvBw 0.06655133 8.532399 0.01
nrdBw 0.28667317 8.533728 0.00
## check stdy with other kernel, distributions
0.0
0.1
0.2
Density
0.3
0.4
0.5
Histogram of x
Figura 2.9: Locpol de denCVBwSelC1
0.4
0.0
0.2
Density
0.6
0.8
Histogram of x
36
CAP
0.4
0.3
0.0
0.1
0.2
Density
0.5
0.6
Histogram of x
Funci
on equivKernel
Descripci
on: calcula el n
ucleo equivalente asociado a la estimacion polinomica
local de la funcion de regresion o sus derivadas.
Funci
on:
equivKernel(kernel,nu,deg,lower=dom(kernel)[[1]],upper=dom(kernel)[[2]],
subdivisions=25)
donde:
nu: es el orden de derivadas a estimar.
on: esta funcion devuelve un vector cuyas componentes son el n
ucleo equivalente utilizado para calcular el estimador polinomial
local de las derivadas de la funcion de regresion.
Ejemplo de uso:
>
>
>
>
curve(EpaK(x),-3,3,ylim=c(-.5,1))
f <- equivKernel(EpaK,0,3)
curve(f(x),-3,3,add=TRUE,col="blue")
curve(gaussK(x),-3,3,add=TRUE)
37
2.3. LIBRO LOCPOL
> f <- equivKernel(gaussK,0,3)

> curve(f(x),-3,3,add=TRUE,col="blue")
> ## Draw several Equivalent local polynomial kernels
0.5
0.0
EpaK(x)
0.5
1.0
curve(EpaK(x),-3,3,ylim=c(-.5,1))
for(p in 1:5){
curve(equivKernel(gaussK,0,p)(x),-3,3,add=TRUE)
}
0.0
EpaK(x)
0.5
1.0
Figura 2.12: Locpol de equivKernel
0.5
>
>
+
+
38
0.5
0.0
EpaK(x)
0.5
1.0
CAP
0.5
0.0
EpaK(x)
0.5
1.0
2.3. LIBRO LOCPOL
39
Funci
on KernelChars
Descripci
on: para un determinado n
ucleo devuelve algunos de los valores numericos mas com
unmente utilizados para funciones relacionadas con ellos.
Funciones: aqu tenemos las mismas que en el punto 2.2.1.
-- RK(K)
-- RdK(K)
-- mu2K(K)
-- mu0K(K)
-- K4(K)
-- dom(K)
Valor devuelto por las funciones: estas funciones devuelven un valor numerico.
Ejemplo de uso:
>
+
+
+
+
+
+
+
+
+
+
+
>
g <- function(kernels)
{
mu0 <- sapply(kernels,function(x) computeMu0(x,))
mu2 <- sapply(kernels,function(x) computeMu(2,x))
Rk.ok <- sapply(kernels,RK)
RK <- sapply(kernels,function(x) computeRK(x))
K4 <- sapply(kernels,function(x) computeK4(x))
res <- data.frame(mu0,mu0.ok,mu2,mu2.ok,RK,Rk.ok,K4)
res
}
g(kernels=c(EpaK,gaussK,TriweigK,TrianK))
1
2
mu0 mu0.ok
mu2
mu2.ok
RK
Rk.ok
K4
1
1 0.2000000 0.200000 0.6000000 0.6000000 0.4337657
1
1 1.0000000 1.000000 0.2820948 0.2820948 0.1994711
40
3
4
CAP
1
1
1 0.1111111 0.111111 0.8158508 0.8158510 0.5879019

1 0.1666667 0.166667 0.6666667 0.6666670 0.4793655
Funci
on kernelCte
Descripci
on: estos son los valores en funcion del n
ucleo y del grado del polinomio local que se utilizan en la seleccion del ancho de banda, tal como se propone
en Fan y Gijbels, 1996.
Funciones:
cteNuK(nu,p,kernel,lower=dom(kernel)[[1]],upper=dom(kernel)[[2]],
subdivisions= 25)
adjNuK(nu,p,kernel,lower=dom(kernel)[[1]],upper=dom(kernel)[[2]],
subdivisions= 25)
donde:
p: es el grado del polinomio local a estimar.
Valor devuelto por las funciones: ambas funciones devuelven valores numericos.
Funci
on Kernels
Descripci
on: definicion de los n
ucleos utilizados en la estimacion polinomial
local.
Funciones:
------
CosK(x)
EpaK(x)
Epa2K(x)
gaussK(x)
...
2.3. LIBRO LOCPOL
41
Funci
on locCteWeights
Descripci
on: pesos asociados al estimador polinomial local de grado cero (o
grado uno, locLinWeightsC.
Funciones:
locCteWeightsC(x, xeval, bw, kernel, weig = rep(1, length(x)))
locLinWeightsC(x, xeval, bw, kernel, weig = rep(1, length(x)))
locWeightsEval(lpweig, y)
locWeightsEvalC(lpweig, y)
donde:
xeval: es un vector con los puntos de evaluacion de los pesos.
lpweig: son pesos (X T W X) 1X T W evaluados en la matriz xeval.
Valor devuelto por las funciones:
locCteWeightsC y locLinWeightsC devuelven una lista con dos componentes
que son:
den : estimacion de (n bw f (x))p + 1.
locWeig : (X T W X) 1X T W evaluado en la matriz xeval.
locWeightsEvalC y locWeightsEval devuelven un vector con la estimacion. Y
realizan el producto matricial entre locWeig y y, para obtener la estimacion
dada por los puntos de xeval.
Ejemplo de uso:
> size <- 200
> sigma <- 0.25
42
>
>
>
>
>
>
>
>
>
>
>
>
CAP
deg <- 1
kernel <- EpaK
bw <- .25
xeval <- 0:100/100
regFun <- function(x) x^3
x <- runif(size)
y <- regFun(x) + rnorm(x, sd = sigma)
d <- data.frame(x, y)
lcw <- locCteWeightsC(d$x, xeval, bw, kernel)$locWeig
lce <- locWeightsEval(lcw, y)
lceB <- locCteSmootherC(d$x, d$y, xeval, bw, kernel)$beta0
mean((lce-lceB)^2)
[1] 1.582489e-32
>
>
>
>
llw <- locLinWeightsC(d$x, xeval, bw, kernel)$locWeig

lle <- locWeightsEval(llw, y)
lleB <- locLinSmootherC(d$x, d$y, xeval, bw, kernel)$beta0
mean((lle-lleB)^2)
[1] 5.319263e-32
Funci
on locpol (caso 1)
Descripci
on: formula para el calculo asociado la estimacion polinomial local.
Funci
on:
locpol(formula,data,weig=rep(1,nrow(data)),bw=NULL,kernel=EpaK,deg=1,
xeval=NULL,xevalLen=100)
confInterval(x)
## S3 method for class locpol:
residuals(object,...)
fitted(object,deg=0,...)
2.3. LIBRO LOCPOL
43
summary(object,...)
print(x,...)
plot(x,...)
donde:
xevalLen: es la longitud de xeval si este no se especifica.
on: devuelve una lista que contiene entre otros
componentes los siguientes:
mf: modelo de marco para los datos y la formula.
data: hoja de datos con datos.
weig: vector de pesos para cada observaciones.
xeval: vector de puntos de evaluacion.
bw: suavizado de parametros, ancho de banda.
kernel: n
ucleo usado, ver n
ucleos.
kName: nombre del n
ucleo, una cadena con el nombre del n
ucleo.
deg: estimacion del grado (p) del polinomio local.
X,Y: nombre de los datos de la respuesta y de la covarianza. Tambien, se
utilizan en lpF it para el nombre de los datos ajustados.
residuals: residuos del polinomio local ajustado.
lpFit: hoja de datos con el ajuste polinomico local, que contiene la covariable,
la respuesta, la derivada estimada, la densidad estimada X y la estimacion de
la varianza.
44
CAP
Ejemplo de uso:
> N <- 250

> xeval <- 0:100/100
## ex1
d <- data.frame(x = runif(N))
d$y <- d$x^2 - d$x + 1 + rnorm(N, sd = 0.1)
r <- locpol(y~x,d)
plot(r)
## ex2
d <- data.frame(x = runif(N))
d$y <- d$x^2 - d$x + 1 + (1+d$x)*rnorm(N, sd = 0.1)
r <- locpol(y~x,d)
plot(r)
## length biased data !!
d <- data.frame(x = runif(10*N))
d$y <- d$x^2 - d$x + 1 + (rexp(10*N,rate=4)-.25)
posy <- d$y[ whichYPos <- which(d$y>0) ];
d <- d[sample(whichYPos, N,prob=posy,replace=FALSE),]
rBiased <- locpol(y~x,d)
r <- locpol(y~x,d,weig=1/d$y)
plot(d)
points(r$lpFit[,r$X],r$lpFit[,r$Y],type="l",col="blue")
points(rBiased$lpFit[,rBiased$X],rBiased$lpFit[,rBiased$Y],type="l")
curve(x^2 - x + 1,add=TRUE,col="red")
95% Conf. Int. for xPoints
1.1
+
+
+
+
++
++
+
1.0
+
++
+
+
+
+
+
+
+ +
+
0.9
++
+
++
+
+
+
+
+
0.8
+ +
+
+ +
+
+
+
++
+
+
+
+
+ +
+
+
+ +
+
+ ++ +
++
+
+ +
+
+
+ ++
++ + +
+
+
0.6
+
+
+
+
+
++
+
++
+
+
++
+
+
+
+
+
+
+
+
+++
++
+
+
++
+
++ + +
+
+
++
++
+
++ +
++
+
+
+
+
+
+
+
+
++
++
+
+
+
+ +
+
++ + +
+
+
+
++
+
+
+ +
++
+
++
+ +
++
+
+
+
+
+
+
+ +
+
+
+
+
+
+
+
+
++
+
+ +
+
+
+
+ +
0.7
+
+
+
+++
++
+
+
++
+
+
+
0.5
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
0.0
0.2
0.4
0.6
0.8
1.0
Figura 2.16: Locpol de locpol1
45
2.3. LIBRO LOCPOL
95% Conf. Int. for xPoints
1.4
1.2
++
+
+
+
+
+
+
+ +
++
1.0
0.8
+
+
+
+
+
+
+
++
+
+
++
+
+
+ +
+ +
+
+
+ +
+
+
+ + +
+
+
+
+
+ +
+ ++
++
+ + +
+ + ++ +
++ +
+
+
+
++
+
+
+
+
+
+ ++
+
++
+ +
+ +
+ +
+
+
+ +
+
+
+
+
+ +
+
+
+
+
+
+
++
++
+
+
+
+
+ +
+ +
++
+
+
+
+
+
+
+ ++
+
+
+
+
+
++
+
+ +
+ +
++
+
+
+ +
+
++
+
+
++
+
+
+ +
++
+
+ +
+ + ++ +
++
+
+
+
+
0.6
+
+
+
+ +
+
+
++ ++
+ +
+
+
++
+
+
+
+
++ +
+
++ +
+
+
+
+
+
+
+
+
+
++
0.4
+
+
0.0
0.2
0.4
0.6
0.8
1.0
0.5
1.0
1.5
2.0
0.0
0.2
0.4
0.6
0.8
1.0
0.5
1.0
1.5
2.0
0.0
0.2
0.4
0.6
0.8
1.0
46
CAP
Funci
on locpol (caso 2)
Descripci
on: calcula directamente la estimacion polinomial local de la funcion
de regresion.
Funci
on:
locCteSmootherC(x, y, xeval, bw, kernel, weig = rep(1, length(y)))
locLinSmootherC(x, y, xeval, bw, kernel, weig = rep(1, length(y)))
locCuadSmootherC(x, y, xeval, bw, kernel, weig = rep(1, length(y)))
locPolSmootherC(x, y, xeval, bw, deg, kernel, DET = FALSE,
weig = rep(1, length(y)))
looLocPolSmootherC(x, y, bw, deg, kernel, weig = rep(1, length(y)),
DET = FALSE)
donde:
DET: es el Booleano que se pide para calcular el calculo del determinante si
la matriz es X T W X.
Valor devuelto por las funciones: devuelven una hoja de datos cuyos componentes dan los puntos de la evaluacion, el estimador de regresion para la funcion
m(x) y sus derivadas en cada punto, y la estimacion de la densidad marginal de x
para el grado p + 1. Estos componentes estan dados por:
x: evaluacion de puntos.
beta0, beta1, beta2,...: estimacion de la i-esima derivada de la funcion de regresion.
den: estimacion de (n bw f (x))p + 1.
Ejemplo de uso:
> N <- 100
> xeval <- 0:10/10
> d <- data.frame(x = runif(N))
47
2.3. LIBRO LOCPOL
>
>
>
>
d
bw <- 0.125
fx <- xeval^2 - xeval + 1
fx
[1] 1.00 0.91 0.84 0.79 0.76 0.75 0.76 0.79 0.84 0.91 1.00
> ## Non random
> d$y <- d$x^2 - d$x + 1
> d$y
[1]
[8]
[15]
[22]
[29]
[36]
[43]
[50]
[57]
[64]
[71]
[78]
[85]
[92]
[99]
0.7612712
0.7856257
0.9431480
0.7895587
0.7903537
0.9376779
0.7538522
0.7935019
0.9235209
0.9384282
0.7550881
0.9153462
0.7521302
0.7874924
0.7850518
0.7559905
0.9741636
0.9800950
0.7531085
0.7556701
0.7546431
0.7552131
0.7500432
0.9114959
0.7897166
0.9455074
0.8034103
0.9085868
0.8297123
0.7722768
0.7522013
0.9471402
0.9716219
0.8553279
0.7860397
0.8301950
0.7965644
0.8074889
0.7629582
0.7922242
0.9149992
0.8810769
0.9313345
0.7909301
0.8848619
0.9154741
0.8589399
0.7591280
0.8094275
0.7556629
0.8778957
0.8835923
0.8350024
0.7500432
0.9838134
0.7504405
0.7854696
0.7507330
0.8160407
0.8486005
0.7505613
0.9631150
0.7536389
0.8893358
0.8928852
0.7947588
0.8640522
0.7503208
0.8332296
0.9838679
0.7500239
0.7631920
0.7507679
0.7595106
0.7940482
0.9631849
0.7867965
0.9162690
0.7509552
0.7504915
0.7598914
0.9851169
0.7501226
0.7711603
0.9594925
0.9656503
> cuest <- locCuadSmootherC(d$x, d$y ,xeval, bw, EpaK)

> cuest
1
2
3
4
5
6
7
8
9
10
11
x beta0
beta1 beta2
den
0.0 1.00 -1.00000e+00
2 0.02372248
0.1 0.91 -8.00000e-01
2 8.02909644
0.2 0.84 -6.00000e-01
2 23.31658125
0.3 0.79 -4.00000e-01
2 11.78540653
0.4 0.76 -2.00000e-01
2 38.67861685
0.5 0.75 1.93709e-16
2 35.01524608
0.6 0.76 2.00000e-01
2 28.50127243
0.7 0.79 4.00000e-01
2 6.56710512
0.8 0.84 6.00000e-01
2 19.20445392
0.9 0.91 8.00000e-01
2 8.74252631
1.0 1.00 1.00000e+00
2 0.01881323
0.7964555
0.7562525
0.7525843
0.8603828
0.8181295
0.7812961
0.8604513
0.7627245
0.8609709
0.7852533
0.8569621
0.8152093
0.8544892
0.8223636
48
CAP
> lpest2 <- locPolSmootherC(d$x, d$y , xeval, bw, 2, EpaK)

> lpest2
x beta0
beta1 beta2
1 0.0 1.00 -1.00000e+00
2
2 0.1 0.91 -8.00000e-01
2
3 0.2 0.84 -6.00000e-01
2
4 0.3 0.79 -4.00000e-01
2
5 0.4 0.76 -2.00000e-01
2
6 0.5 0.75 9.48996e-16
2
7 0.6 0.76 2.00000e-01
2
8 0.7 0.79 4.00000e-01
2
9 0.8 0.84 6.00000e-01
2
10 0.9 0.91 8.00000e-01
2
11 1.0 1.00 1.00000e+00
2
> print(cbind(x = xeval, fx, cuad0 = cuest$beta0, lp0 = lpest2$beta0,
cuad1 = cuest$beta1, lp1=lpest2$beta1))
[1,]
[2,]
[3,]
[4,]
[5,]
[6,]
[7,]
[8,]
[9,]
[10,]
[11,]
x
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
fx cuad0 lp0
cuad1
lp1
1.00 1.00 1.00 -1.00000e+00 -1.00000e+00
0.91 0.91 0.91 -8.00000e-01 -8.00000e-01
0.84 0.84 0.84 -6.00000e-01 -6.00000e-01
0.79 0.79 0.79 -4.00000e-01 -4.00000e-01
0.76 0.76 0.76 -2.00000e-01 -2.00000e-01
0.75 0.75 0.75 1.93709e-16 9.48996e-16
0.76 0.76 0.76 2.00000e-01 2.00000e-01
0.79 0.79 0.79 4.00000e-01 4.00000e-01
0.84 0.84 0.84 6.00000e-01 6.00000e-01
0.91 0.91 0.91 8.00000e-01 8.00000e-01
1.00 1.00 1.00 1.00000e+00 1.00000e+00
> d$y <- d$x^2 - d$x + 1 + rnorm(d$x, sd = 0.1)

> d$y
[1]
[8]
[15]
[22]
[29]
[36]
[43]
0.6723636
0.8505989
0.9601393
0.8048367
0.7696472
0.8876989
0.6733276
0.8448013
1.0953380
1.0679246
0.6544773
0.7438772
0.9026008
0.6779314
0.6243922
0.9551620
1.1700375
0.9297912
0.8958657
0.8190183
0.7221547
0.8162543
0.8478272
0.7825367
0.7927135
0.8276265
0.9027349
0.8123858
0.7436355
0.9471580
0.6269945
1.0375346
0.7252702
0.6894363
0.8134340
0.6866786
0.8222345
0.7835987
1.1104370
0.8883434
0.9984543
0.8745750
0.7271204
0.7517749
0.7351159
0.7743744
0.6869968
0.8354297
0.8664553
49
2.3. LIBRO LOCPOL
[50]
[57]
[64]
[71]
[78]
[85]
[92]
[99]
0.8047627
0.8460743
0.8374226
0.7706758
0.9301271
0.8154858
0.7293871
0.8457371
0.8801028
1.0463443
0.7110618
0.9308003
0.8300982
1.0819505
0.8928125
0.8609747
0.6856138
0.6829823
0.8485211
0.9222196
0.8937102
0.7909390
0.8333529
0.9211008
0.8607412
0.7864358
1.0207491
0.5815393
0.9901345
0.8743793
0.8376500
0.9401312
0.7215284
0.8385355
0.9806565
0.5212752
0.6063388
0.7492864
0.7307411
0.9580189
0.7530619
0.6800853
1.0151708
0.8478267
> cuest <- locCuadSmootherC(d$x,d$y , xeval, bw, EpaK)

> cuest
x
1 0.0
2 0.1
3 0.2
4 0.3
5 0.4
6 0.5
7 0.6
8 0.7
9 0.8
10 0.9
11 1.0
beta0
1.1696954
0.8791549
0.8136628
0.7869036
0.7310876
0.7329695
0.7619366
0.8089766
0.8248833
0.8967816
1.0870074
beta1
-4.82196709
-1.52413568
-0.54206496
-0.73315011
-0.22747304
-0.06783130
0.51090384
0.07681645
0.62602660
0.85090352
5.53515036
beta2
41.562041
24.858571
6.356705
-9.195076
7.882866
4.586677
2.573981
-5.305648
7.108892
8.014109
86.530261
den
0.02372248
8.02909644
23.31658125
11.78540653
38.67861685
35.01524608
28.50127243
6.56710512
19.20445392
8.74252631
0.01881323
> lpest2 <- locPolSmootherC(d$x,d$y , xeval, bw, 2, EpaK)

> lpest2
x
1 0.0
2 0.1
3 0.2
4 0.3
5 0.4
6 0.5
7 0.6
8 0.7
9 0.8
10 0.9
11 1.0
beta0
1.1696954
0.8791549
0.8136628
0.7869036
0.7310876
0.7329695
0.7619366
0.8089766
0.8248833
0.8967816
1.0870074
beta1
-4.82196709
-1.52413568
-0.54206496
-0.73315011
-0.22747304
-0.06783130
0.51090384
0.07681645
0.62602660
0.85090352
5.53515036
beta2
41.562041
24.858571
6.356705
-9.195076
7.882866
4.586677
2.573981
-5.305648
7.108892
8.014109
86.530261
> lpest3 <- locPolSmootherC(d$x,d$y , xeval, bw, 3, EpaK)
0.5999255
0.6731402
0.6532467
0.7516319
0.7677295
1.0160035
0.8197416
50
CAP
> lpest3
x
1 0.0
2 0.1
3 0.2
4 0.3
5 0.4
6 0.5
7 0.6
8 0.7
9 0.8
10 0.9
11 1.0
beta0
0.5804087
0.8791243
0.8144113
0.7906162
0.7218013
0.7316335
0.7684981
0.8052345
0.8250382
0.8984336
1.0834538
beta1
beta2
beta3
39.0917900 -1761.581230 32073.140200
-1.5169100
24.921877
-8.022927
-0.7355235
5.935609
148.187404
0.1022015
-8.019882 -742.843339
1.1560714
5.543538 -1033.073360
-0.5044145
5.830562
410.750916
0.8838808
1.012286 -278.143321
-1.1792922
1.192569 1065.659358
0.6029547
7.136487
18.045586
1.0882845
6.575599 -266.262808
5.2542255
75.285771 -191.756862
> cbind(x = xeval, fx, cuad0 = cuest$beta0, lp20 = lpest2$beta0,

+ lp30 = lpest3$beta0, cuad1 = cuest$beta1, lp21 = lpest2$beta1,
+ lp31 = lpest3$beta1)
[1,]
[2,]
[3,]
[4,]
[5,]
[6,]
[7,]
[8,]
[9,]
[10,]
[11,]
x
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
fx
1.00
0.91
0.84
0.79
0.76
0.75
0.76
0.79
0.84
0.91
1.00
cuad0
1.1696954
0.8791549
0.8136628
0.7869036
0.7310876
0.7329695
0.7619366
0.8089766
0.8248833
0.8967816
1.0870074
lp20
1.1696954
0.8791549
0.8136628
0.7869036
0.7310876
0.7329695
0.7619366
0.8089766
0.8248833
0.8967816
1.0870074
lp30
0.5804087
0.8791243
0.8144113
0.7906162
0.7218013
0.7316335
0.7684981
0.8052345
0.8250382
0.8984336
1.0834538
cuad1
-4.82196709
-1.52413568
-0.54206496
-0.73315011
-0.22747304
-0.06783130
0.51090384
0.07681645
0.62602660
0.85090352
5.53515036
lp21
-4.82196709
-1.52413568
-0.54206496
-0.73315011
-0.22747304
-0.06783130
0.51090384
0.07681645
0.62602660
0.85090352
5.53515036
lp31
39.0917900
-1.5169100
-0.7355235
0.1022015
1.1560714
-0.5044145
0.8838808
-1.1792922
0.6029547
1.0882845
5.2542255
> print(cbind(x = xeval, fx, cuad0 = cuest$beta0, lp20 = lpest2$beta0,

+ lp30 = lpest3$beta0, cuad1 = cuest$beta1, lp21 = lpest2$beta1,
+ lp31 = lpest3$beta1))
Funci
on pluginBw
Descripci
on: implementa un selector del ancho de banda de tipo plug-in para
la funcion de regresion. El metodo se describe en el texto de Fan y Gijbels (1996)
2.3. LIBRO LOCPOL
51
paginas 110-112, y utiliza estimaciones piloto seg

un al regla del pulgar (Rule of
thumb.
Funci
on:
pluginBw(x, y, deg, kernel, weig = rep(1,length(y)))
on: esta funcion devuelve un valor numerico.
Ejemplo de uso:
>
>
>
>
>
>
>
>
>
size <- 200

sigma <- 0.25
deg <- 1
kernel <- EpaK
xeval <- 0:100/100
x <- runif(size)
> cvBwSel <- regCVBwSelC(d$x,d$y, deg, kernel, interval = c(0, 0.25))

> cvBwSel
[1] 0.1962416
> thBwSel <- thumbBw(d$x, d$y, deg, kernel)
> thBwSel
[1] 0.1993822
> piBwSel <- pluginBw(d$x, d$y, deg, kernel)
> piBwSel
[1] 0.1268736
> est <- function(bw, dat, x) return(locPolSmootherC(dat$x,dat$y, x,
+ bw, deg, kernel)$beta0)
> ise <- function(val, est) return(sum((val - est)^2 * xeval[[2]]))
52
CAP
> plot(d$x, d$y)

> trueVal <- regFun(xeval)
> trueVal
[1]
[9]
[17]
[25]
[33]
[41]
[49]
[57]
[65]
[73]
[81]
[89]
[97]
>
>
>
>
>
>
>
>
>
>
+
>
>
0.000000
0.000512
0.004096
0.013824
0.032768
0.064000
0.110592
0.175616
0.262144
0.373248
0.512000
0.681472
0.884736
0.000001
0.000729
0.004913
0.015625
0.035937
0.068921
0.117649
0.185193
0.274625
0.389017
0.531441
0.704969
0.912673
0.000008
0.001000
0.005832
0.017576
0.039304
0.074088
0.125000
0.195112
0.287496
0.405224
0.551368
0.729000
0.941192
0.000027
0.001331
0.006859
0.019683
0.042875
0.079507
0.132651
0.205379
0.300763
0.421875
0.571787
0.753571
0.970299
0.000064
0.001728
0.008000
0.021952
0.046656
0.085184
0.140608
0.216000
0.314432
0.438976
0.592704
0.778688
1.000000
0.000125
0.002197
0.009261
0.024389
0.050653
0.091125
0.148877
0.226981
0.328509
0.456533
0.614125
0.804357
lines(xeval, trueVal, col = "red")

xevalRes <- est(cvBwSel, d, xeval)
cvIse <- ise(trueVal, xevalRes)
lines(xeval, xevalRes, col = "blue")
xevalRes <- est(thBwSel, d, xeval)
thIse <- ise(trueVal, xevalRes)
xevalRes <- est(piBwSel, d, xeval)
piIse <- ise(trueVal, xevalRes)
lines(xeval, xevalRes, col = "blue", lty = "dashed")
res <- rbind( bw = c(cvBwSel, thBwSel, piBwSel),
ise = c(cvIse, thIse, piIse) )
colnames(res) <- c("CV", "th", "PI")
res
CV
th
PI
bw 0.196241577 0.199382210 0.126873646
ise 0.001127976 0.001100554 0.002027269
0.000216
0.002744
0.010648
0.027000
0.054872
0.097336
0.157464
0.238328
0.343000
0.474552
0.636056
0.830584
0.000343
0.003375
0.012167
0.029791
0.059319
0.103823
0.166375
0.250047
0.357911
0.493039
0.658503
0.857375
53
0.5
0.5
0.0
d$y
1.0
2.3. LIBRO LOCPOL
0.0
0.2
0.4
0.6
0.8
1.0
d$x
0.5
0.0
0.5
d$y
1.0
Figura 2.20: Locpol de pluginBw1
0.0
0.2
0.4
0.6
0.8
1.0
d$x
Figura 2.21: Locpol de pluginBw2
54
CAP
Funci
on PRDenEstC
Descripci
on: estimador Parsen-Rosenblat de densidad univariante.
Funci
on:
PRDenEstC(x, xeval, bw, kernel, weig = rep(1, length(x)))
on: esta funcion devuelve una hoja de datos de
la forma: (x, den).
Ejemplo de uso:
>
>
>
>
>
N <- 100
x <- runif(N)
xeval <- 0:10/10
b0.125 <- PRDenEstC(x, xeval, 0.125, EpaK)
b0.125
1
2
3
4
5
6
7
8
9
10
11
x
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
den
0.3762050
0.8861180
0.9478700
0.9940558
1.2625850
1.2979144
1.1465889
1.2102660
0.8060191
0.7589499
0.3192066
> b0.05 <- PRDenEstC(x, xeval, 0.05, EpaK)

> b0.05
1
2
3
4
x
0.0
0.1
0.2
0.3
den
0.4603164
1.0064206
0.5727967
1.1694444
2.3. LIBRO LOCPOL
5
6
7
8
9
10
11
0.4
0.5
0.6
0.7
0.8
0.9
1.0
55
1.0153751
1.4813057
0.9270388
1.4077981
0.5424289
1.1261434
0.2701843
> cbind(x = xeval, fx = 1, b0.125 = b0.125$den, b0.05 = b0.05$den)
[1,]
[2,]
[3,]
[4,]
[5,]
[6,]
[7,]
[8,]
[9,]
[10,]
[11,]
x fx
b0.125
b0.05
0.0 1 0.3762050 0.4603164
0.1 1 0.8861180 1.0064206
0.2 1 0.9478700 0.5727967
0.3 1 0.9940558 1.1694444
0.4 1 1.2625850 1.0153751
0.5 1 1.2979144 1.4813057
0.6 1 1.1465889 0.9270388
0.7 1 1.2102660 1.4077981
0.8 1 0.8060191 0.5424289
0.9 1 0.7589499 1.1261434
1.0 1 0.3192066 0.2701843
Funci
on regCVBwSelC
Descripci
on: implementa el metodo por validacion cruzada para la seleccion
del ancho de banda del estimador de regresion polinomial local.
Funci
on:
regCVBwSelC(x, y, deg, kernel=gaussK,weig=rep(1,length(y)),
interval=.lokestOptInt)
on: esta funcion devuelve un valor numerico, el
ancho de banda seleccionado.
Ejemplo de uso:
> size <- 200
> sigma <- 0.25
56
>
>
>
>
>
>
>
CAP
deg <- 1
kernel <- EpaK
xeval <- 0:100/100
x <- runif(size)
> cvBwSel <- regCVBwSelC(d$x,d$y, deg, kernel, interval = c(0, 0.25))

> cvBwSel
[1] 0.1962416
> thBwSel
[1] 0.1993822
> piBwSel
[1] 0.1268736
>
+
>
>
>
>
est <- function(bw, dat, x) return(locPolSmootherC(dat$x,dat$y, x,

bw, deg, kernel)$beta0)
ise <- function(val, est) return(sum((val - est)^2 * xeval[[2]]))
plot(d$x, d$y)
trueVal <- regFun(xeval)
trueVal
[1]
[9]
[17]
[25]
[33]
[41]
[49]
[57]
[65]
[73]
[81]
[89]
0.000000
0.000512
0.004096
0.013824
0.032768
0.064000
0.110592
0.175616
0.262144
0.373248
0.512000
0.681472
0.000001
0.000729
0.004913
0.015625
0.035937
0.068921
0.117649
0.185193
0.274625
0.389017
0.531441
0.704969
0.000008
0.001000
0.005832
0.017576
0.039304
0.074088
0.125000
0.195112
0.287496
0.405224
0.551368
0.729000
0.000027
0.001331
0.006859
0.019683
0.042875
0.079507
0.132651
0.205379
0.300763
0.421875
0.571787
0.753571
0.000064
0.001728
0.008000
0.021952
0.046656
0.085184
0.140608
0.216000
0.314432
0.438976
0.592704
0.778688
0.000125
0.002197
0.009261
0.024389
0.050653
0.091125
0.148877
0.226981
0.328509
0.456533
0.614125
0.804357
0.000216
0.002744
0.010648
0.027000
0.054872
0.097336
0.157464
0.238328
0.343000
0.474552
0.636056
0.830584
0.000343
0.003375
0.012167
0.029791
0.059319
0.103823
0.166375
0.250047
0.357911
0.493039
0.658503
0.857375
57
2.3. LIBRO LOCPOL
[97] 0.884736 0.912673 0.941192 0.970299 1.000000
>
>
>
>
>
>
>
>
>
>
+
>
>

res
0.5
0.5
0.0
d$y
1.0
CV
th
PI
bw 0.196241577 0.199382210 0.126873646
ise 0.001127976 0.001100554 0.002027269
0.0
0.2
0.4
0.6
0.8
1.0
d$x
Figura 2.22: Locpol de regCVBwSelC
58
CAP
Funci
on selKernel
Descripci
on: utilizacion de los atributos del n
ucleo para seleccionar n
ucleos.
Esta funcion se utiliza principalmente para fines internos.
Funci
on:
selKernel(kernel)
on: esta funcion devuelve un entero que es u
nico
para cada n
ucleo.
Funci
on simpleSmoothers
Descripci
on: calcula un suavizador elemental de las respuestas, definido como
el promedio ponderado por el n
ucleo de las respuestas observadas. La version simpleSqSmootherC es igual solo que promediando las respuestas al cuadrado.
Funci
on:
simpleSmootherC(x, y, xeval, bw, kernel, weig = rep(1, length(y)))
simpleSqSmootherC(x, y, xeval, bw, kernel)
Valor devuelto por las funciones: ambas funciones devuelven una hoja de
datos con:
x: puntos de evaluacion de x.
reg: valores de suavizado en los puntos de x.
Ejemplo de uso:
>
>
>
>
size <- 1000

x <- runif(100)
bw <- 0.125
kernel <- EpaK
2.3. LIBRO LOCPOL
>
>
>
>
>
59
xeval <- 1:9/10

y <- rep(1,100)
## x kern. aver. should give density f(x)
prDen <- PRDenEstC(x,xeval,bw,kernel)$den
prDen
[1] 1.0888622 0.8859709 0.8785186 1.0515013 0.9550039 1.3555096 1.1455056

[8] 0.5906089 0.7851241
> ssDen <- simpleSmootherC(x,y,xeval,bw,kernel)$reg
> ssDen
[1] 1.0888622 0.8859709 0.8785186 1.0515013 0.9550039 1.3555096 1.1455056
[8] 0.5906089 0.7851241
> all(abs(prDen-ssDen)<1e-15)
[1] TRUE
> ## x kern. aver. should be f(x)*R2(K) aprox.
> s2Den <- simpleSqSmootherC(x,y,xeval,bw,kernel)$reg
> s2Den
[1] 0.6397285 0.5396532 0.5139364 0.6437790 0.5363702 0.8570298 0.6951637
[8] 0.3604343 0.4755898
> summary( abs(prDen*RK(kernel)-s2Den) )
Min. 1st Qu.
Median
Mean 3rd Qu.
Max.
0.004515 0.007860 0.012880 0.016280 0.013590 0.043720
> summary( abs(1*RK(kernel)-s2Den) )
Min. 1st Qu. Median
Mean 3rd Qu.
Max.
0.03973 0.06035 0.08606 0.11220 0.12440 0.25700
>
>
+
+
+
## x kern. aver. should be f(x)*R2(K) aprox.

for(n in c(1000,1e4,1e5))
{
s2D <- simpleSqSmootherC(runif(n),rep(1,n),xeval,bw,kernel)$reg
cat("\n",n,"\n")
60
CAP
+ print( summary( abs(1*RK(kernel)-s2D) ) )

+ }
1000
Min.
1st Qu.
Median
Mean
3rd Qu.
Max.
0.0009869 0.0150400 0.0264600 0.0371300 0.0441600 0.1071000
10000
Min. 1st Qu.
Median
Mean 3rd Qu.
Max.
0.001878 0.007167 0.010570 0.013100 0.021160 0.024560
1e+05
Min.
1st Qu.
Median
Mean
3rd Qu.
Max.
0.0008332 0.0015290 0.0049780 0.0051580 0.0092530 0.0117500
Funci
on thumbBw
Descripci
on: implementa la regla del pulgar descrita en el texto de Fan y Gijbels
(1996) para la seleccion del ancho de banda en regresion polinomial local.
Funci
on:
thumbBw(x, y, deg, kernel, weig = rep(1, length(y)))
compDerEst(x, y, p, weig = rep(1, length(y)))
Valor devuelto por las funciones: la funcion thumbBw devuelve un u
nico
valor numerico, mientras que compDerEst devuelve una hoja de datos cuyos componentes son los siguientes:
x: valores de x.
y: valores de y.
res: residuos de la estimacion parametrica.
61
2.3. LIBRO LOCPOL
der: derivada estimada de los valores de x.

Ejemplo de uso:
>
>
>
>
>
>
>
>
>
>
>
size <- 200

sigma <- 0.25
deg <- 1
kernel <- EpaK
xeval <- 0:100/100
x <- runif(size)
cvBwSel <- regCVBwSelC(d$x,d$y, deg, kernel, interval = c(0, 0.25))
cvBwSel
[1] 0.1962416
> thBwSel
[1] 0.1993822
> piBwSel
[1] 0.1268736
>
+
>
>
>
>
est <- function(bw, dat, x) return(locPolSmootherC(dat$x,dat$y, x,

bw, deg, kernel)$beta0)
ise <- function(val, est) return(sum((val - est)^2 * xeval[[2]]))
plot(d$x, d$y)
trueVal <- regFun(xeval)
trueVal
[1]
[9]
[17]
[25]
[33]
[41]
0.000000
0.000512
0.004096
0.013824
0.032768
0.064000
0.000001
0.000729
0.004913
0.015625
0.035937
0.068921
0.000008
0.001000
0.005832
0.017576
0.039304
0.074088
0.000027
0.001331
0.006859
0.019683
0.042875
0.079507
0.000064
0.001728
0.008000
0.021952
0.046656
0.085184
0.000125
0.002197
0.009261
0.024389
0.050653
0.091125
0.000216
0.002744
0.010648
0.027000
0.054872
0.097336
0.000343
0.003375
0.012167
0.029791
0.059319
0.103823
62
CAP
[49]
[57]
[65]
[73]
[81]
[89]
[97]
>
>
>
>
>
>
>
>
>
>
+
>
>
0.110592
0.175616
0.262144
0.373248
0.512000
0.681472
0.884736
0.117649
0.185193
0.274625
0.389017
0.531441
0.704969
0.912673
0.125000
0.195112
0.287496
0.405224
0.551368
0.729000
0.941192
0.132651
0.205379
0.300763
0.421875
0.571787
0.753571
0.970299
0.140608
0.216000
0.314432
0.438976
0.592704
0.778688
1.000000
0.148877
0.226981
0.328509
0.456533
0.614125
0.804357

res
0.5
0.5
0.0
d$y
1.0
CV
th
PI
bw 0.196241577 0.199382210 0.126873646
ise 0.001127976 0.001100554 0.002027269
0.0
0.2
0.4
0.6
0.8
1.0
d$x
Figura 2.23: Locpol de thumbBw
0.157464
0.238328
0.343000
0.474552
0.636056
0.830584
0.166375
0.250047
0.357911
0.493039
0.658503
0.857375
2.4. LIBRO LOCFIT
2.4.
63
Libro Locfit
En este libro llamado Locfit se implementan funciones asociadas al calculo de

la regresion polinomial local, la verosimilitud y la densidad estimada. Los metodos
implementados se pueden consultar en los textos de Cook y Weisberg (1994), Loader
(1999). Dicho libro corresponde a la actualizacion del 17-04-2009, version 1.5-4.
El libro contiene muchas funciones, en este trabajo nos centramos tan solo en
aquellas asociadas al estimador de la densidad y de regresion de tipo polinomial local.
As como aquellas que permiten la seleccion del parametro de suavizado asociado.
Funciones que calculan la densidad estimada

Funcion density.lf : proporciona un interfaz para el calculo de la funcion de
densidad.
Uso de la funci
on:
density.lf(x, n = 50, window = "gaussian", width, from, to,
cut = if(iwindow == 4.) 0.75 else 0.5,
ev = lfgrid(mg = n, ll = from, ur = to),
deg = 0, family = "density", link = "ident", ...)
donde:
window: es el tipo de ventana que se utiliza para la estimacion.
from: es el lmite inferior para el dominio de estimacion.
to: es el lmite superior para el dominio de estimacion.
cut: es la expansion de los controles por defecto del dominio.
on: esta funcion devuelve una lista con los
componentes de x (puntos de evaluacion) e y (densidad estimada).
Ejemplo de uso:
> data(geyser)
> density.lf(geyser, window="tria")
64
CAP
$x
[1]
[8]
[15]
[22]
[29]
[36]
[43]
[50]
1.564723
2.060516
2.556310
3.052103
3.547897
4.043690
4.539484
5.035277
$y
[1]
[6]
[11]
[16]
[21]
[26]
[31]
[36]
[41]
[46]
1.182697e-16
5.571365e-01
1.963319e-01
1.046841e-02
9.107183e-02
1.022381e-01
5.441002e-01
8.292375e-01
4.034270e-01
1.409679e-01
1.635550
2.131344
2.627137
3.122931
3.618724
4.114518
4.610311
1.706378
2.202172
2.697965
3.193759
3.689552
4.185346
4.681139
1.890007e-01
3.739426e-01
1.132892e-01
7.019261e-02
4.838509e-02
1.911011e-01
5.917874e-01
6.788447e-01
3.677768e-01
9.322311e-02
1.777206
2.272999
2.768793
3.264586
3.760380
4.256173
4.751967
1.848033
2.343827
2.839620
3.335414
3.831207
4.327001
4.822794
4.623314e-01
1.256785e-01
3.419681e-02
4.762960e-02
9.322311e-02
3.630067e-01
3.663234e-01
4.238963e-01
5.616119e-01
9.322311e-02
1.918861
2.414654
2.910448
3.406241
3.902035
4.397828
4.893622
7.174016e-01
3.316671e-03
7.653125e-02
2.512419e-02
9.322311e-02
3.338104e-01
3.868570e-01
4.480279e-01
7.364176e-01
5.972420e-02
1.989689
2.485482
2.981276
3.477069
3.972863
4.468656
4.964450
8.081279e-01
8.002072e-02
4.129095e-02
1.445726e-01
6.792613e-02
2.442562e-01
6.537848e-01
4.545461e-01
3.963838e-01
5.913483e-17
> density(geyser, window="tria")

Call:
density.default(x = geyser, window = "tria")
Data: geyser (107
x
Min.
:0.5668
1st Qu.:1.9334
Median :3.3000
Mean
:3.3000
3rd Qu.:4.6666
Max.
:6.0332
obs.);
Bandwidth bw = 0.3677
y
Min.
:0.0000
1st Qu.:0.0605
Median :0.1560
Mean
:0.1828
3rd Qu.:0.2568
Max.
:0.4834
Funcion gam.lf : es una llamada a la funcion locfit por lf() utilizada en los
terminos del modelo aditivo. Normalmente no es llamada directamente por los
usuarios.
Uso de la funci
on:
gam.lf(x, y, w, xeval, ...)
2.4. LIBRO LOCFIT
65
donde:
w: son los pesos a priori.
Funcion locfit: formula base del libro para el ajuste de la regresion local y
para los modelos de verosimilitud.
Uso de la funci
on:
locfit(formula, data=sys.frame(sys.parent()), weights=1, cens=0, base=0,
subset, geth=FALSE, ..., lfproc=locfit.raw)
donde:
formula: es la formula del modelo; por ejmplo, y~lp(x) para un modelo de
regresion; ~lp(x) para un modelo de densidad estimada. Es recomendable, usar para lp() los RHS, especialmente cuando por defecto no se
utilizan los parametros de suavizacion.
subset: son las observaciones de subconjuntos en la hoja de datos.
geth: no utilizar; siempre se pone geth = F ALSE.
lfproc: es una de las funciones de procesamiento para calcular el ajuste
local. El valor predeterminado es locfit.raw(). Otras opciones incluyen
locfit.robust(), locfit.censor() y locfit.quasi().
on: esta funcion devuelve un objeto con clase
locfit y un conjunto estandar de metodos para impresion, graficos, etc.
Ejemplo de uso:
>
>
>
>
# fit and plot a univariate local regression

data(ethanol, package="locfit")
fit <- locfit(NOx ~ E, data=ethanol)
plot(fit, get.data=TRUE)
> # a bivariate local regression with smaller smoothing parameter

> fit <- locfit(NOx~lp(E,C,nn=0.5,scale=0), data=ethanol)
> plot(fit)
66
CAP
# density estimation
data(geyser, data="locfit")
fit <- locfit( ~ lp(geyser, nn=0.1, h=0.8))
plot(fit,get.data=TRUE)
o o
o
o
o
o
oo
o
o
o
o
o
o
ooo
o
o o
o oo
o o
o
o
oo
o
oo
NOx
o
o o
o
o
o
o
o
o
o
o
o
oo
oo
o oo
o
o
oo
o
o
o
o
o
ooo o o
o oo
o oo o
o
oo o
o
o
o oo
oo
0.6
0.7
0.8
0.9
1.0
1.1
1.2
16
1.5
2.5
18
Figura 2.24: Locfit de locfit1
0.5
10
12
14
0.5
0.8
0.7
2.5
1.5
0.5
0.6
3.5
0.9
1.0
1.1
1.2
0.4
0.3
0.1
0.2
density
0.5
0.6
0.7
0.0
>
>
>
>
2.0
2.5
3.0
3.5
4.0
4.5
5.0
geyser
2.4. LIBRO LOCFIT
67
Funciones que calculan el par

ametro de suavizado
Funcion sjpi: calcula un ancho de banda a traves del metodo plug-in propuesto
por Sheather y Jones (1991) para la estimacion de una densidad.
Uso de la funci
on:
sjpi(x, a)
donde:
a: es un vector de anchos de banda piloto.
on: esta funcion devuelve una matriz con cuatro
columnas; el n
umero de filas es igual a la longitud de a. La primera columna
es el plug-in de ancho de banda seleccionado. La segunda columna es el ancho
de banda piloto a. La tercera columna es el ancho de banda piloto de acuerdo
con la supuesta relacion de Sheather-Jones. La cuarta columna es un calculo
intermedio.
Ejemplo de uso:
>
>
>
>
data(geyser)
gf <- 2.5
a <- seq(0.05, 0.7, length=100)
z <- sjpi(geyser, a)
>
>
>
+
# the plug-in curve. Multiplying by gf=2.5 corresponds

# to Locfits standard scaling for the Gaussian kernel.
plot(gf*z[, 2], gf*z[, 1], type = "l", xlab = "Pilot Bandwidth k",
ylab = "Bandwidth h")
>
>
>
+
# Add the assumed curve.

lines(gf * z[, 3], gf * z[, 1], lty = 2)
legend(gf*0.05, gf*0.4, lty = 1:2,
legend = c("Plug-in", "SJ assumed"))
68
1.0
CAP
0.6
0.2
0.4
Bandwidth h
0.8
Plugin
SJ assumed
0.5
1.0
1.5
Pilot Bandwidth k
Figura 2.27: Locfit de sjpi
2.5.
Libro sm
En este libro incluye metodos de regresion no parametrica y estimacion n

ucleo
de la densidad descritos en el texto de Bowman y Azzalini (1997). Dicho libro corresponde a la actualizacion del 19-04-2009, version es del 2.2-3.
El libro contiene m
ultiples funciones que no seran descritas aqui. De hecho nos
centramos tan solo en aquellas asociadas al estimador de la densidad y de regresion
de tipo polinomial local, en concreto aquellas que permiten la seleccion del parametro
de suavizado asociado.
Estimacion n
ucleo de la funci
on de densidad
La funcion sm.density calcula la estimacion de la densidad de los datos en
una, dos o tres dimensiones. En dos dimensiones, una variedad de representaciones
graficas pueden ser seleccionadas, y en tres dimensiones, puede representarse una
superficie de contorno.
uso de la funci
on:
sm.density(x, h, model = "none", weights = NA, group=NA, ...)
on: esta funcion devuelve una lista que contiene los
69
2.5. LIBRO SM
valores de la estimacion de la densidad en los puntos de evaluacion, el parametro de

suavizacion, los pesos del parametro de suavizacion y los pesos del n
ucleo. Tambien,
para los datos de uno o dos dimensiones, se proporciona el error estandar de la
estimacion y los extremos sueprior e inferior de una banda de la variabilidad.
Ejemplo de uso:
>
>
>
>
# A one-dimensional example
y <- rnorm(50)
sm.density(y, model = "Normal")
# sm.density(y, panel = TRUE)
> # A two-dimensional example

> y <- cbind(rnorm(50), rnorm(50))
> sm.density(y, display = "image")
0.4
0.3
0.2
0.0
0.1
Probability density function
0.5
> # A three-dimensional example

> y <- cbind(rnorm(50), rnorm(50), rnorm(50))
> sm.density(y)
y[2]
Figura 2.28: sm de sm-density1
y[1]
70
CAP
C
alculo del par
ametro de suavizado
La funcion h.select selecciona un parametro de suavizacion para la densidad
estimada en una o dos dimensiones y para la regresion no parametrica con una
o dos covariables.
Uso de la funci
on:
h.select(x, y = NA, weights = NA, group = NA, ...)
on: esta funcion devuelve el valor del parametro
de suavizacion seleccionado.
Ejemplo de uso:
> x <- rnorm(50)
> h.select(x)
[1] 0.4766912
> h.select(x, method = "sj")
[1] 0.4164856
71
2.5. LIBRO SM
> x <- matrix(rnorm(100), ncol = 2)

> h.select(x)
[1] 0.5315978 0.4961265
> sm.density(x, method = "cv")
> x <- rnorm(50)
> y <- x^2 + rnorm(50)
> h.select(x, y)
[1] 0.3869072
> sm.regression(x, y, method = "aicc")
> x <- matrix(rnorm(100), ncol = 2)
> y <- x[,1]^2 + x[,2]^2 + rnorm(50)
> h.select(x, y, method = "cv", structure.2d = "common")
[1] 0.4695942 0.4695942
> sm.regression(x, y, df = 8)
Density fun
0.10
0.05
ction
0.00
2
1
]
x[2
1
1
0
]
x[1
Figura 2.31: sm de sm1
72
CAP
4
2
2
1
x[2
]
0
0
1
1
1
[1]
2.5. LIBRO SM
73
La funcion hcv utiliza la tecnica de validacion cruzada para seleccionar un

parametro de suavizacion adecuado para la construccion de la densidad estimada o la curva de regresion no parametrica en una o dos dimensiones.
Uso de la funci
on:
hcv(x, y = NA, hstart = NA, hend = NA, ...)
donde:
hstart: es el valor mas peque
no de los puntos de la red para ser utilizado en
una b
usqueda en la red inicial para el valor del parametro de suavizacion.
hend: es el mayor valor de los puntos de la red para ser utilizado en una
b
usqueda en la red inicial para el valor del parametro de suavizacion.
de suavizacion que minimiza el criterio de validacion cruzada en la cuadrcula
seleccionada.
Ejemplo de uso:
>
>
>
>
>
# Density estimation
x <- rnorm(50)
par(mfrow=c(1,2))
h.cv <- hcv(x, display="lines", ngrid=32)
sm.density(x, h=hcv(x))
>
>
>
>
>
>
# Nonparametric regression
x <- seq(0, 1, length = 50)
y <- rnorm(50, sin(2 * pi * x), 0.2)
par(mfrow=c(1,2))
h.cv <- hcv(x, y, display="lines", ngrid=32)
sm.regression(x, y, h=hcv(x, y))
74
0.2
0.4
h
0.15
0.6
0.8
0.25
0.0
0.8
CAP
0.4
0.05
Probability density function

y
0.24
0.26
0.28
0.6
0.4
0.2
0.0
1.0
0.5
0.0
0.5
1.0
CV
CV
0.30
0.32
0.34
0.36
10
8
6
4
2
2.5. LIBRO SM
75
La funcion hnorm eval

ua el parametro de suavizacion, asintoticamente optimo,
para estimar una funcion de densidad cuando la distribucion subyacente es
normal.
Uso de la funci
on:
hnorm(x, weights)
de suavizacion asintoticamente optimo para el caso normal.
Ejemplo de uso:
> x <- rnorm(50)
> hnorm(x)
[1] 0.4273578
La funcion hsj utiliza el metodo plug-in de Sheather-Jones para la seleccion
de un parametro de suavizacion, que es adecuado para la construccion de la
densidad estimada unidimensional.
Uso de la funci
on:
hsj(x)
de suavizacion encontrado por el metodo de Sheather-Jones.
Ejemplo de uso:
> x <- rnorm(50)
> hsj(x)
[1] 0.3350144
76
CAP
Estimaci
on de la funci
on de regresi
on
Funcion sm.discontinuity: utiliza una comparacion de izquierda a derecha
de las curvas de regresion no parametricas para evaluar la evidencia de la
presencia de una o mas discontinuidades en una curva de regresion o en la
superficie.
Uso de la funci
on:
sm.discontinuity(x, y, h, hd, ...)
donde:
hd: es un parametro de suavizado que se utiliza para suavizar las diferencias del lado izquierdo y derecho de las estimaciones de la regresion no
parametrica.

on: esta funcion devuelve una lista que contiene
los siguientes elementos:
p: es el p-valor para la prueba de la hipotesis nula de que no se presentan
discontinuidades.
sigma: es la desviacion estandar estimada de los errores.
eval.points: son los puntos de evaluacion de las estimaciones de la regresion no parametrica.
st.diff: es un vector o matriz de las diferencias estandarizadas entre el
lado izquierdo y derecho de los estimadores en los puntos de evaluacion.
diffmat: cuando x es un vector, este contiene la ubicacion y las diferencias
estandarizadas, donde este u
ltimo es mayor que 2,5.
angle: cuando x es una matriz, este contiene los angulos estimados en las
diferencias estandarizadas que fueron construidas.
h: es el parametro de suavizacion principal.
hd: es el parametro de suavizacion utilizado para la doble suavizacion.
77
2.5. LIBRO SM
Ejemplo de uso:
> par(mfrow = c(3, 2))
> provide.data(nile)
> sm.discontinuity(Year, Volume, hd = 0)
Test of continuity:
significance =
0.006
location st.diff
1888.5 -2.85
1889.5 -3.65
1890.5 -3.12
1891.5 -2.82
1896.5 3.78
1897.5 4
1898.5 4.77
1899.5 3.35
1915.5 -3.02
1938.5 2.58
> sm.discontinuity(Year, Volume)
Test of continuity:
significance =
0.009
location st.diff
1887.5 -2.9
1888.5 -3.28
1889.5 -3.34
1890.5 -3.03
1896.5 3.47
1897.5 4.2
1898.5 4.54
1899.5 4.48
1900.5 4.06
1901.5 3.37
1902.5 2.51
>
>
>
>
ind <- (Year > 1898)

plot(Year, Volume)
h <- h.select(Year, Volume)
sm.regression(Year[ ind], Volume[ ind], h, add = TRUE)
78
CAP
>
>
>
+
+
+
+
hvec <- 1:15

p <- numeric(0)
for (h in hvec) {
result <- sm.discontinuity(Year, Volume, h,
display = "none", verbose = 0)
p <- c(p, result$p)
}
>
>
>
>
>
>
>
plot(hvec, p, type = "l", ylim = c(0, max(p)), xlab = "h")

lines(range(hvec), c(0.05, 0.05), lty = 2)
provide.data(trawl)
Position <- cbind(Longitude, Latitude)
ind <- (Longitude < 143.8)
sm.regression(Position[ind,], Score1[ind], theta = 35, phi = 30)
sm.discontinuity(Position[ind,], Score1[ind], col = "blue")
Test of continuity:
significance =
0.016
> # The following example takes longer to run.

> # Alternative values for nside are 32 and 64.
> # Alternative values of yjump are 1 and 0.5.
>
>
>
>
>
>
>
>
>
>
>
>
par(mfrow = c(1, 1))

nside <- 16
yjump <- 2
x1 <- seq(0, 1, length = nside)
x2 <- seq(0, 1, length = nside)
x <- expand.grid(x1, x2)
x <- cbind(x1 = x[, 1], x2 = x[, 2])
y <- rnorm(nside * nside)
ind <- (sqrt((x[, 1] - 0.5)^2 + (x[, 2] - 0.5)^2) <= 0.25)
y[ind] <- y[ind] + yjump
image(x1, x2, matrix(y, ncol = nside))
sm.discontinuity(x, y, df = 20, add = TRUE)
Test of continuity:
significance =
0.005
79
1900
1920
1940
1000
Volume
1880
600
1000
600
Volume
2.5. LIBRO SM
1960
1880
1900
1920
1940
1960
Year
0.00
0.15
1000
600
1880
1900
1920
1940
1960
10
Score1[ind]
Latitude
1.5
1.0
0.5
0.0
11.2
e
11.4
Lon
143.0
ud
143.2
11.6
gitu 11.8
143.4
tit
143.6
de La
14
2.5
143.0
143.2
143.4
Longitude
0.6
0.8
1.0
Figura 2.36: sm de sm-discontinuity1
0.4
2.
3.5
0.0
0.2
4.5
0.0
12
11.8 11.4
Year
x2
Volume
0.30
Year
0.2
0.4
0.6
0.8
1.0
x1
Figura 2.37: sm de sm-discontinuity2
143.6
143.8
80
CAP
Funcion sm.monotonicity utiliza la idea de un .ancho de banda crticopara

evaluar la evidencia de que una curva de regresion es no monotona.
Uso de la funci
on:
sm.monotonicity(x, y, N = rep(1, length(y)), h,
type = "continuous", ...)
on: esta funcion devuelve una lista que contiene
los siguientes elementos:
p: es el p-valor para la prueba de la hipotesis nula que la curva de verdad
es monotona.
hcrit: es el parametro de suavizacion crtico. Este es el valor mas peque
no
que cuando se aplica a los datos observados, hace la curva monotona.
h: es el parametro de suavizacion utilizado para la doble suavizacion.
Ejemplo de uso:
> # Radiocarbon dating data
> provide.data(radioc)
Data file being loaded
> ind <- (Cal.age>5000 & Cal.age<6000)
> cal.age <- Cal.age[ind]
> rc.age <- Rc.age[ind]
> sm.monotonicity(cal.age, rc.age, method = "aicc", nboot = 200)
The smallest h which produces an increasing curve is 56.436
Standard deviation of estimated errors is: 25.603
Smoothing parameter used for estimation of residuals is: 27.875
Test of monotonicity: p = 0.06
> # Hosmer & Lemeshow birth data
> provide.data(birth)
81
2.5. LIBRO SM
> sm.monotonicity(Lwt[Smoke == "N"], Low[Smoke == "N"],

+ type = "binomial")
rc.age
4600
4800
5000
5200
The smallest h which produces a decreasing curve is 36.495

Test of monotonicity: p = 0.35
5000
5200
5400
5600
5800
6000
cal.age
0.6
0.4
0.0
0.2
Low[Smoke == "N"]
0.8
1.0
Figura 2.38: sm de sm-monotonicity1
100
150
200
Lwt[Smoke == "N"]
Figura 2.39: sm de sm-monotonicity2
82
CAP
Funcion sm.regression crea una estimacion de la regresion no parametrica

de los datos, que constan de una sola variable respuesta y una o dos covariables. En dos dimensiones se produce una perspectiva, una imagen (image),
un contorno (slice) o una parcela rgl de la superficie de regresion estimada.
Uso de la funci
on:
sm.regression(x, y, h, design.mat = NA, model = "none",
weights = NA, group = NA, ...)
donde:
design.mat: es la matriz de dise
no utilizada para producir y, cuando estos
se suponen que son los residuos de un modelo lineal.
on: esta funcion devuelve una lista que contiene los valores de la estimacion de los puntos de evaluacion, el parametro de
suavizacion y los pesos del parametro de suavizacion. Si se ha especificado un
modelo de referencia y la prueba se establece en TRUE, entonces se devuelve
el p-valor de la prueba.
Cuando solo hay una covariable, tambien se devuelven los pesos asociados
con las diferentes observaciones, una estimacion de la desviacion del error
estandar y el error estandar de la estimacion. Si se ha especificado un modelo
de referencia, el error estandar se refiere a la comparacion entre la estimacion
y el modelo de referencia y la definicion de los valores del modelo de referencia
que seran devueltos tambien.
Ejemplo de uso:
> provide.data(trawl)
>
>
>
>
Zone92 <- (Year == 0 & Zone == 1)

Position <- cbind(Longitude - 143, Latitude)
dimnames(Position)[[2]][1] <- "Longitude - 143"
par(mfrow = c(2, 2))
> sm.regression(Longitude, Score1, method = "aicc", col = "red",

+ model = "linear")
83
2.5. LIBRO SM
Test of linear model:
significance =
> sm.regression(Position[Zone92, ], Score1[Zone92],

+ display = "image", theta = 120)
> sm.regression(Position[Zone92, ], Score1[Zone92], df = 12,
+ col = "se", theta = 120)
> sm.regression(Position[Zone92, ], Score1[Zone92], df = 12,
+ col = "se", model = "linear", theta = 120)
significance =
0.031
> par(mfrow = c(1, 1))

> sm.regression(Position[Zone92, 2:1], Score1[Zone92],
+ display = "rgl", df = 12)
> sm.regression(Position[Zone92, 2:1], Score1[Zone92],
+ display = "rgl", df = 12, alpha = c(0.9, 1), col = "se",
+ model = "linear")
significance =
0.031
> sm.regression(Position[Zone92, 1], Score1[Zone92],

+ panel = TRUE)
+ panel = TRUE)
+ panel = TRUE, display = "rgl")
84
11.4
11.6
Latitude
11.8
1.0
0.0
143.6
0.0
Longitude
0.2
0.4
0.6
0.8
Longitude 143
de
43
1.5
1.0
0.5
0.0
0.0
0.2
0.4
11.8
11.7
0.6
L11.6
ati11.5
tu11.4
0.8
de11.3
itu
43
de
itu
ng
one92]
Score1[Z
one92]
Score1[Z
1.5
1.0
0.5
0.0
0.0
0.2
0.4
11.8
11.7
0.6
L11.6
ati11.5
tu11.4
0.8
de11.3
ng
143.2
Lo
142.8
Lo
Score1
2.0
CAP
Figura 2.40: sm de sm-regression1
85
2.5. LIBRO SM
h = 0.0948
1.0
0.5
0.0
Score1[Zone92]
1.5
2.0
df = 6
0.0
0.2
0.4
0.6
0.8
Position[Zone92, 1]

df = 12
h = ( 0.167 , 0.0863 )
1.5
Score1[Zon
1.0
e92]
0.5
0.0
titu
La
11.3
11.4
11.5
11.6
de
0.8
11.7
11.8
0.0
0.6
0.4 143
0.2 itude
g
Lon
86
CAP
Captulo 3
Aplicaci
on pr
actica
3.1.
Estudio con datos reales
En el siguiente estudio vamos a considerar los datos que se encuentran en el

libro MASS, en concreto el objeto de tipo hoja de datos geyser. Estos datos
corresponden a un estudio sobre las erupciones de los geiseres del Old Faithful en el
Parque Nacional de Yellowstone, Wyoming.
La version que aqu consideramos se describe en el libro de Azzalini y Bowman (1990) y corresponde a la medicion continua desde el 1 de Agosto hasta el
15 de Agosto de 1985. Ademas, algunas mediciones de la duracion nocturna fueron
codificadas como 2, 3 o 4 minutos, tras haber sido inicialmente descrito como corto,
mediano o largo.
Geyser es una hoja de datos con 299 observaciones y dos variables, las cuales
son: duration que es el tiempo del n
umero de erupciones (en minutos) y waiting que
es el tiempo de espera para la proxima erupcion.
Con los datos de la hoja de datos geyser, procederemos a ilustrar los metodos de
regresion no parametrica univariantes descritos en el captulo 1 de este trabajo. Si
bien ilustraremos otros metodos notables (en este caso consideraremos un suavizamiento por Spline seg
un Heckman y Ramsay 1996), como ya se indico en dicho
captulo, nuestro interes se centrara fundamentalmente en el estimador polin
omico
local. Para ello utilizaremos las funciones disponibles en los libros KernSmooth y
locpol, que fueron ampliamente descritos en el captulo 2.
Para ajustar el estimador polinomico local de grado p utilizaremos la funcion
locpoly del libro KernSmooth. El uso de dicha funcion lo haremos considerando en
primer lugar una eleccion arbitraria del parametro de suavizado o ancho de banda.
En este caso hemos considerado h = 0,25, no obstante tal y como ilustraremos
despues, es posible utilizar elecciones automaticas, mas refinadas, como los criterio
de seleccion cross-validation y plug-in discutidos en el captulo 1.
En cuanto a la funcion n
ucleo considerada, la funcion locpoly por defecto usa
87
88
PRACTICA
CAP
ITULO 3. APLICACION
n
ucleos normales (argumento kernel =normal), y en este ejemplo hemos dejado
dicha eleccion por defecto. De este modo compararemos el resultado usando diferentes grados para lo cual hemos generado el siguiente codigo:
> data(geyser, package = "MASS")

> names(geyser)
[1] "waiting"
"duration"
> x <- geyser$duration

> y <- geyser$waiting
> plot(x, y)
> fit <- locpoly(x, y, bandwidth = 0.25,degree=1)
> lines(fit,col="red")
> lines(fit,col="blue")
> lines(fit,col="green")
> lines(fit,col="orange")
> legend(1,60,legend=c("p=1", "p=0", "p=2", "p=3"),lty=1,
+ col=c("red", "blue", "green", "orange"))
89
60
70
80
90
100
110
3.1. ESTUDIO CON DATOS REALES
50
p=1
p=0
p=2
p=3
Figura 3.1: Estimador polinomico local de grado p para los datos del geyser. La
variable x es duracion y la variable y es tiempo de espera y el tama
no de la hoja de
datos es n = 299
Como se puede observar en la Figura 3.1 conforme vamos aumentando el grado

del estimador polinomico local, las estimaciones son mas irregulares, intentando
capturar en mayor medida las observaciones consideradas. Notese que como ya se
discutio en el captulo 1, esto supone estimaciones que pagan la disminucion en el
sesgo con un incremento notable de la variabilidad. Tambien es notable ver como
el incremento de p = 1 hasta p = 2 no supone una mejora del modelo (ni en sesgo
ni en variabilidad), de hecho como tambien se discutio al inicio de este trabajo es
preferible usar grados impares frente a los inmediatamente consecutivos pares (ver el
libro de Fan y Gijbels 1996 para una discusion mas detallada y algunos resultados
numericos).
Finalmente podemos ver como las diferencias entre el estimador de NadarayaWatson que considera ajustes locales constantes (p = 0) presenta las mayores diferencias en la proximidades a las fronteras, esto es debido a que los ajustes lineales
locales (p = 1) permiten una correccion automatica de los efectos frontera (para
mas detalles se puede ver el libro de Wand y Jones 1994 o Fan y Gijbels 1996).
A continuacion vamos a ilustrar otro metodo de suavizamiento para dichos datos,
en concreto un estimador de tipo spline. Existen funciones para dicho proposito en
varios libros de R (SemiPar, ssplines, esplines, etc.), ademas de la funcion smooth.spline
dentro del libro base stats. En este caso hemos considerado esta u
ltima funcion y
la hemos comparado con el resultado ofrecido por la funcion sm.spline, que implementa el estimador descrito en Heckman y Ramsay (1996), y que esta contenida
90
PRACTICA
CAP
ITULO 3. APLICACION
en el libro pspline. Dicho estimador se define con un parametro de suavizado que

por defecto considera un criterio basado en validacion cruzada o validacion cruzada
generalizada.
Nosotros dejamos las definiciones por defecto que considera dicha funcion. A
efectos comparativos tambien hemos incluido el estimador lineal local con ancho
de banda plug-in (mas detalles sobre dicho procedimiento se daran en el ejemplo
siguiente, Figura 3.3). De este modo, el codigo generado para dicho ajuste y los
resultados ofrecidos son los siguientes:
110
y <- geyser$waiting
plot(x,y,title=Suavizado por Splines para los datos del geyser,
+ xlab=Duraci
on,ylab=Tiempo de espera)
geyser.spl <- sm.spline(x,y)
geyser.spl
lines(geyser.spl, lty=1,lwd=2, col = "green")
geyser.sts<-smooth.spline(x,y)
geyser.sts
lines(geyser.spl, lty=2,lwd=2, col = "blue")
h1 <- dpill(x, y)
fit1 <- locpoly(x,y,bandwidth = h1, degree=1)
lines(fit1,col="red",lty=3,lwd=2)
legend(topright,legend=c("sm.spline", "smooth.spline",
+ "locpoly (plug-in)"), lwd=2,lty=1:3,
+ col=c("green","blue","red"))
80
70
50
60
Tiempo de espera
90
100
sm.spline
smooth.spline
locpoly (plugin)
Duracin
Figura 3.2: Estimador de tipo spline.
91
La Figura 3.2 muestra los ajustes realizados. Podemos observar que los resultados
graficos de las funciones smooth.spline y sm.spline son identicos, si bien los algoritmos implementos difieren ligeramente y los resultados de los objetos generados
tambien. En concreto se obtiene:
> geyser.spl
Call:
smooth.Pspline(x = ux, y = tmp[, 1], w = tmp[, 2], method = method)
Smoothing Parameter (Spar): 5.388047e-05
Equivalent Degrees of Freedom (Df): 32.3831
GCV Criterion: 31.59988
CV Criterion: 60.01679
> geyser.sts
Call:
smooth.spline(x = x, y = y)
Smoothing Parameter spar= 1.044891 lambda= 0.05174406 (12 iterations)
Equivalent Degrees of Freedom (Df): 3.366694
Penalized Criterion: 9604.32
GCV: 110.0222
Con respecto a la bondad de los ajustes vemos que el estimador lineal local ofrece
una estimacion mas suavizada que los splines en este caso.
Nuestro siguiente objetivo sera comparar todos los procedimientos disponibles
para la seleccion del ancho de banda, asociado al estimador lineal local. Los procedimientos para seleccionar el ancho de banda considerado son los metodos plug-in,
validacion cruzada y la sencilla regla del pulgar.
La implementacion de dichos metodos se hace en diversas funciones disponibles
en los libros KernSmooth, Locpol y sm (si bien existen algunas versiones mas
disponibles en otros libros de R que contienen posibilidades para metodos no parametricos, como por ejemplo el libro Locfit que implementa versiones del estimador que
fijan el n
umero de observaciones en el entorno local en lugar del tama
no del mismo).
Todas ellas fueron ampliamente descritas en el captulo 2 de este trabajo, por
lo que aqu nos centraremos en su particular aplicacion a los datos con los que
estamos trabajando. Agrupando las funciones seg
un la metodologa de seleccion que
implementan, podemos enumerar las siguientes:
92
PRACTICA
CAP
ITULO 3. APLICACION
Selectores de tipo plug-in: la funcion dpill que forma parte del libro KernSmooth, implementando el metodo de Ruppert, Sheather y Wand (1995). Y
la funcion pluginBw dentro del libro Locpol, que implementa el metodo
descrito en el libro de Fan y Gijbels (1996) paginas 110-112.
Selectores basados en Validacion Cruzada: la funcion regCVBwSelC, que
se puede encontrar en el captulo 2 de este trabajo en el libro Locpol, y la
funcion h.select del libro sm.
Usando dichas funciones nuestro objetivo es el calcular el estimador lineal local,
usando distintos parametros ancho de banda. En concreto parametros seg
un metodos
plug-in y validacion cruzada. Ademas de ilustrar el resultado final, podremos discutir
las particularidades que conllevan la aplicacion de cada uno de ellos.
Dado que nuestro objetivo en este momento es el parametro de suavizado, volvemos a fijar la eleccion de la funcion n
ucleo de tipo gausiano y, como ya hemos
comentado antes, fijamos ajustes de grado p = 1. El codigo generado para tales
propositos es el siguiente:
y <- geyser$waiting
plot(x,y,title=Estimaci
on lineal local para los datos del geyser,
+ xlab=Duraci
on,ylab=Tiempo de espera)
h1 <- dpill(x, y)
h1
lines(fit1,col="red",lty=1,lwd=2)
h2<- pluginBw(x,y, deg=1, kernel=gaussK)
h2
fit2 <- locpoly(x, y , bandwidth = h2, degree=1)
lines(fit2,col="yellow",lty=2,lwd=2)
h3<- regCVBwSelC(x,y, deg=1, kernel=gaussK)
h3
lines(fit3,col="green",lty=3,lwd=2)
h4<-h.select(x, y,method = "cv")
h4
lines(fit4,col="blue",lty=4,lwd=2)
legend(topright,legend=c("h1-dpill", "h2-pluginBw","h3-regCVBwSelC",
+ "h4-h.select"),lwd=2,lty=1:4,col=c("red", "yellow", "green","blue"))
93
Los resultados obtenidos para los parametros de suavizado son los siguientes:
> h1
[1] 0.2342897
> h2
[1] 0.08805326
> h3
[1] 0.739138
> h4
[1] 0.6789217
110
En el calculo de h2, el parametro seg

un el metodo plug-in dentro del libro locpol
resulta notorio el valor tan peque
no que se obtiene en comparacion con los otros.
Habra que estudiar el procedimiento implementado puesto que si observamos la estimacion resultante (Figura 3.3) la curva estimada sufre de regularidades inadmisibles
en las proximidades de la frontera inicial, debido a la escasez de observaciones.
El resultado correspondiente a los criterios basados en validacion cruzada (h3
y h4) son analogos, observandose leves diferencias que tendra que ver con la implementacion concreta que se ha hecho del metodo (en concreto con la rejilla de
minimizacion definida para el criterio).
En este caso el mejor ajuste viene desde el metodo plug-in que implementa la
funcion dpill del libro KernSmooth.
80
70
50
60
Tiempo de espera
90
100
h1dpill
h2pluginBw
h3regCVBwSelC
h4h.select
Duracin
Figura 3.3: Estimador lineal local con distintos h.
94
PRACTICA
CAP
ITULO 3. APLICACION
3.2.
Estudio con datos simulados
Ilustraremos ahora los metodos de regresion no parametrica y en concreto el estimador polinomial local a partir de datos simulados. Nuestro objetivo ahora sera el
de cuantificar la bondad de las estimaciones (dado que conocemos los modelos exactos) y ademas ilustrar aspectos interesantes del problema de regresion con sera el
del efecto del tama
no muestral y la variabilidad de la muestra considerada.
De este modo estudiaremos el comportamiento de los estimadores con distintos
tama
nos de muestra (n = 25, 50, 100y500) y con distintas desviaciones tpicas para
los residuos del modelo (0,3, 0,4y0,1).
De esta forma lo que se pretende es observar la convergencia de la curva teorica
y asimismo ver como el problema de estimacion se hace mas difcil de resolver a
medida que vamos aumentando el valor de la desviacion tpica de los residuos del
modelo.
Para realizar lo anterior consideramos el siguiente modelo de regresion:
Y = m(x) +
donde
m(x) = sen(2x) + 2exp(16x2 )
y donde x se genera seg
un una distribucion uniforme continua en el intervalo (2, 2)
y los residuos se consideran normales con media 0 y desviacion tpica .
En primer lugar, empezaremos comparando el estimador polinomico local (EPL)
con con grados p = 0, 1, 3. El parametro ancho de banda lo fijamos en h = 0,15.
Y en segundo lugar, de forma similar al ejercicio que hicimos con los datos del
geyser, compararemos el EPL con los distintos metodos de seleccion para el ancho
de banda (plug-in, CV, regla del pulgar), fijando ajustes de grado p = 1.
Para cuantificar la precision de las estimaciones resultantes utilizaremos como criterio de error la suma residual de cuadrados sobre una rejilla de puntos de estimacion.
De este modo evaluaremos el estimador sobre una red de puntos xl , l = 1, ..., ngrid
equiespaciada en (2, 2) de tama
no ngrid = 500. Una vez calculadas las estimaciones sobre la rejilla calcularemos los errores con la formula:
500
X
1
ch (xl ))2
(m(xl ) m
500 i=1
(3.1)
y compararemos los resultados tomando la raz cuadrada.

Consideramos en primera lugar (caso 1) la estimacion considerando muestras
de tama
no n = 100 y = 0,4. El codigo generado, as como los resultados obtenidos
son los siguientes:
95
3.2. ESTUDIO CON DATOS SIMULADOS
n<-100
sigma<-0.4
nucleo<-gaussK
regFun<-function(x) sin(2*x)+2*exp(-16*x^2)
x<-runif(n,-2,2)
mx<-regFun(x)
y<-mx+rnorm(n,mean=0,sd=sigma)
plot(x,y,main="Datos simulados")
curve(sin(2*x)+2*exp(-16*x^2),col="black",lwd=2,add=T)
h<- 0.15
fit0 <- locpoly(x,y,bandwidth = h,degree=0)
lines(fit0,col="orange",lty=2,lwd=2)
fit1 <- locpoly(x,y,bandwidth = h, degree=1)
lines(fit1,col="blue",lty=3,lwd=2)
fit3 <- locpoly(x,y,bandwidth = h, degree=3)
lines(fit3,col="green",lty=4,lwd=2)
legend(topright,legend=c("Curva te
orica", "ajuste p=0",
"ajuste p=1","ajuste p=3"),lwd=2,lty=1:4,col=c("black",
"orange","blue","green"))
Datos simulados
Curva terica
ajuste p=0
ajuste p=1
ajuste p=3
Figura 3.4: Estimacion polinomial local a partir de datos simulados. El tama

no
muestral es 100 y la desviacion tpica residual 0,4
Como se puede observar en la Figura 3.4, el grafico cuando p = 0 y p = 1 son muy

parecidos los estimadores salvo en la frontera, debido a que p = 1 permite corregir
de forma automatica los efectos frontera. Los resultados para p = 3 muestran una
mayor irregularidad.
96
PRACTICA
CAP
ITULO 3. APLICACION
Consideramos ahora una disminucion en el tama

no muestral hasta n = 25 manteniendo el mismo = 0,4 (caso 2). El codigo generado a tal efecto es analogo,
cambiando n<-25, y los resultados obtenidos se muestran en la Figura 3.5.
2.0
Datos simulados
1.0
0.5
0.0
0.5
1.0
1.5
Curva terica
ajuste p=0
ajuste p=1
ajuste p=3
2.0
1.5
1.0
0.5
0.0
0.5
1.0
1.5

no
Como se puede observar en el grafico cuando hemos disminuido el tama
no muestral a 25 vemos que los estimadores presentan bastantes irregularidades, sobre todo
cuando intentamos ajustar un polinomio de grado alto (p = 3).
Considerando ahora un tama
no de muestra intermedio de n = 50 (caso 3), y
usando un codigo similar en R obtenemos los resultados que se muestran en la Figura
3.6.
Datos simulados
Curva terica
ajuste p=0
ajuste p=1
ajuste p=3

no
97
3.2. ESTUDIO CON DATOS SIMULADOS
La convergencia a la curva teorica la podemos observar considerando un tama

no
de muestra elevado como n = 500 (caso 4).
Datos simulados
Curva terica
ajuste p=0
ajuste p=1
ajuste p=3

no
A continuacion ilustraremos el comportamiento de los estimadores lineales locales

con diferentes metodos de seleccion del parametro de suavizado. Nos centramos
en el caso 1 en el que se simularon n = 100 datos con = 0,4. Los metodos
considerados son el selector basado en validacion cruzada calculado usando la funcion
regCVBwSelC, el de tipo plug-in calculado usando pluginBw y el calculado seg
un
la simple regla del pulgar, ofrecido por la funcion thumbBw, todas ellas contenidas
en el libro locpol. La comparacion de los metodos la haremos va la raz cuadrada
del error definido en 3.1. El codigo generado para ello se muestra a continuacion:
n<-100
sigma<-0.4
ngrid<-500
nucleo<-gaussK
xgrid<-runif(ngrid,-2,2)
regFun<-function(x) sin(2*x)+2*exp(-16*x^2)
x<-runif(n,-2,2)
mx<-regFun(x)
y<-mx+rnorm(n,mean=0,sd=sigma)
p<-1
est <- function(h, x,y, xgrid,p,nucleo) return(locPolSmootherC(x,y,
xgrid, h, deg=p,kernel=nucleo)$beta0)
error<-function(val,est)return(sqrt(mean((val-est)^2)))
98
PRACTICA
CAP
ITULO 3. APLICACION
cvBwSel <- regCVBwSelC(x,y, deg=p, kernel=nucleo)

teoricos <- regFun(xgrid)
estimados <- est(cvBwSel, x,y,xgrid,p,nucleo)
cvError <- error(teoricos, estimados)
thBwSel <- thumbBw(x, y, deg=p, kernel=nucleo)
estimados <- est(thBwSel, x,y,xgrid,p,nucleo)
thError <- error(teoricos, estimados)
piBwSel <- pluginBw(x, y, deg=p, kernel=nucleo)
estimados <- est(piBwSel, x,y,xgrid,p,nucleo)
piError <- error(teoricos, estimados)
resultado <- list(n=n,cv=cvError,th=thError,pi=piError)
resultado
Los resultado obtenidos comparando con los tama
nos muestrales n = 25, 50, 100, 500
se muestran de forma resumida en la siguiente tabla:
cv
n= 25
n= 50
n=100
n=500
th
0.3238807
0.2768001
0.2140991
0.07953835
0.2492679
0.2395432
0.2080944
0.09508265
pi
NA
0.2738346
0.2093302
0.08881742
De dichos resultados se puede observar que el comportamiento de los metodos

plug-in es ligeramente superior a validacion cruzada. No obstante la diferencia se
hace menos patente en tama
nos de muestra lmite. Tambien es de destacar que
cuando se consideran pocos datos n = 25 no es posible el calculo del selector de tipo
plug-in. Esto es debido a que dichos metodos requieren estimaciones de las derivadas
que no son posibles en este caso.
Finalmente, repetiremos el ejercicio de comparacion de los selectores variando la
dificultad del problema de estimacion. Esto lo haremos variando la desviacion tpica
de los residuos del modelo, considerando = 0,001, 0,1, 0,5. El tama
no de muestral
lo mantenemos en n = 100.
cv
sigma=
sigma=
sigma=
0.001
0.1
0.5
0.01265973
0.06306383
0.2278124
th
0.08732214
0.0912005
0.1755693
pi
0.03904589
0.06364838
0.1888333
Bibliografa
99
Desde la tabla de resultados obtenidos se observa un paralelismo entre aumentar

el tama
no muestral y reducir la variabilidad residual. Cuando hay alta variabilidad
muetral = 0,5 validacion cruzada ofrece resultados pobres. En este caso son los
selectores de tipo plug-in y en particular la sencilla regla del pulgar.
100
Bibliografa
Bibliografa
[1] Wand, M. P. and Jones, M. C. (1995). Kernel Smoothing. Chapman and
Hall, London.
[2] Wand, M. P. (1994). Fast Computation of Multivariate Kernel Estimators.
Journal of Computational and Graphical Statistics, 3, 433-445.
[3] Sheather, S. J. and Jones, M. C. (1991). A reliable data-based bandwidth
selection method for kernel density estimation. Journal of the Royal Statistical
Society, Series B, 53, 683690.
[4] Scott, D. W. (1979). On optimal and data-based histograms. Biometrika, 66,
605610.
[5] Wand, M. P. (1995). Data-based choice of histogram binwidth. University of
New South Wales, Australian Graduate School of Management Working Paper
Series No. 95011.
[6] Ruppert, D., Sheather, S. J. and Wand, M. P. (1995). An effective bandwidth selector for local least squares regression. Journal of the American Statistical Association, 90, 12571270.
[7] John Fox, 2002. Nonparametric Regression. Appendix to An R and S-PLUS
Companion to Applied Regression.
[8] Fan, J. and Gijbels, I. Local polynomial modelling and its applications. Chapman and Hall, London (1996).
[9] Wand, M. P. and Jones, M. C. Kernel smoothing. Chapman and Hall Ltd.,
London (1995).
[10] Crist
obal, J. A. and Alcal
a, J. T. (2000). Nonparametric regression estimators for length biased data. J. Statist. Plann. Inference, 89, pp. 145-168.
[11] Ahmad, Ibrahim A. (1995). On multivariate kernel estimation for samples
from weighted distributions. Statistics and Probability Letters, 22, num. 2, pp.
121-129
101
102
Bibliografa
[12] H
ardle W. (1990). Smoothing techniques. Springer Series in Statistics, New
York (1991).
[13] Loader, C. (1999). Local Regression and Likelihood. Springer, New York.
[14] Consult the Web page http://www.locfit.info/.
[15] Cleveland, W. and Grosse, E. (1991). Computational Methods for Local
Regression. Statistics and Computing 1.
[16] Sheather, S. J. and Jones, M. C. (1991). A reliable data-based bandwidth
selection method for kernel density estimation. JRSS-B 53, 683-690.
[17] Bowman, A.W. and Azzalini, A. (1997). Applied Smoothing Techniques for
Data Analysis: the Kernel Approach with S-Plus Illustrations. Oxford University
Press, Oxford.
[18] Hurvich, C.M., Simonoff, J.S. and Tsai, C.-L. (1998). Smoothing parameter selection in nonparametric regression using an improved Akaike information
criterion. J. R. Statistic. Soc., Series B, 60, 271-293.
[19] Bowman, A.W., Pope, A. and Ismail, B. (2006). Detecting discontinuities
in nonparametric regression curves and surfaces. Statistics and Computing, 16,
377390.
[20] Bowman, A.W., Jones, M.C. and Gijbels, I. (1998). Testing monotonicity
of regression. J.Comp.Graph.Stat. 7, 489-500.
[21] Bowman, A.W. (2006). Comparing nonparametric surfaces. Statistical Modelling, 6, 279-299.
[22] Hastie, T.J. and Tibshirani R.J. Generalized Additive Models. Chapman
and Hall. (2000).
[23] Fan, J., Gijbels, I. and Hu, T.-C. and Huang, L.-S. (1996). An asymptotic study of variable bandwidth selection for local polynomial regression with
application to density estimation. Statistica Sinica, Vol. 6, No. 1.
[24] Wand, M.P. and Jones, M.C. (1994). Multivariate Plug-in Bandwidth Selection. Computational Statistics, 9. pp. 97-116.
[25] Wand, M.P. and Jones, M.C. (1995). Kernel Smoothing. Monographs on
Statistics and Applied Probability 60. Ed. Chapman and Hall.
[26] Azzalini, A. and Bowman, A. W. (1990). A look at some data on the Old
Faithful geyser. Applied Statistics 39, 357-365.
Bibliografa
103
[27] Venables, W. N. and Ripley, B. D. (2002). Modern Applied Statistics with

S. Fourth edition. Springer.
[28] Heckman, N. and Ramsay, J. O. (1996). Spline smoothing with model based
penalties. McGill University, unpublished manuscript.
[29] Martnez Miranda, M.D., Raya Miranda, R., Gonz
alez Manteiga, W.
and Gonz
alez Carmona, A. (2008). A bootstrap local bandwidth selector for
additive models. Journal of Computational and Graphical Statistics, 17,38-55.
[30] Linton, O.B., and Nielsen, J.P. (1995). A Kernel Method of Estimating Structured Nonparametric Regression Based on Marginal Integration.
Biometrika, 82, 93100.
[31] Kim,W., Linton, O.B., and Hengartner, N.W. (1999). A Computationally
Efficient Oracle Estimator for Additive Nonparametric Regression with Bootstrap Confidence Intervals. Journal of Computational and Graphical Statistics,
8, 278297.
[32] Kauermann, G., and Opsomer, J.D. (2003). Local Likelihood Estimation in
Generalized Additive Models. Scandinavian Journal of Statistics, 30, 317337.
[33] Nielsen, J.P., and Sperlich, S. (2005). Smooth Backfitting in Practise. Journal of the Royal Statistical Society, Ser. B, 67, 4361.
[34] Mammen, E., and Park, C. (2005). Bandwidth Selection for Smooth Backfitting in Additive Models. The Annals of Statistics, 33, 12601294.
[35] Severance-Lossin, E., and Sperlich, S. (1997). Estimation of Derivatives for Additive Separable Models. Discussion paper, SBF 373. HumboldtUniversity, Berlin.
[36] Opsomer, J.D., and Ruppert, D. (1997). Fitting a Bivariate Additive Model
by Local Polynomial Regression. The Annals of Statistics, 25, 186211.
[37] Nadaraya, E.A. (1964). On estimating regression. Theory Probab. Appl, No.9,
pp. 141-142.
[38] Watson, G.S. (1964). Smooth regression analysis. Sankhya Ser. A, No. 26, pp.
101-116.
[39] Eubank, R.L. (1988). Spline Smoothing and Nonparametric Regression. Marcel Dekker, New York.
104
BIBLIOGRAF
IA
[40] Cleveland, W.S. (1979). Robust Locally Weighted Regression and Smoothing
Scatterplots. Journal of the American Statistical Association, Vol. 74, No. 368.
Theory and Methods, pp. 829-836.
[41] Ruppert,D. Wand, M.P. and Carroll, R.J. (2003). Semiparametric Regression. Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge
University Press.
[42] Ruppert, D. and Wand, M.P. (1994). Multivariate Locally Weighted Least
Squares Regression. The Annals of Statistics, Vol. 22, No. 3, pp. 1346-1370.
[43] Bellman, R.E. (1961). Adaptive control processes. Princeton University Press.
[44] Friedman, J.H. and Stuetzle, W. (1981). Projection pursuit regression.
Journal of the American Statistical Association, Vol. 76, No. 376, pp. 817-823.
[45] Breiman, L. and Friedman, J.H. (1985). Estimating optimal transformations for multiple regression and correlation (with discussion). Journal of the
American Statistical Association, Vol. 80, pp. 580-619.
[46] Hastie, T.J. and Tibshirani R. (1990). Generalized additive models. Washington, D.C.;Chapman and Hall.
[47] Buja, A., Hastie, T.J. and Tibshirani, R. (1989). Linear smoothers and
additive models (with discussion). The Annals of Statistics, Vol. 17, No. 2, pp.
453555.
[48] Kim,W., Linton, O.B., and Hengartner, N.W. (1997). A nimble
method of estimating additive nonparametric regression. Electronic article,
http://www.stats.yale.edu.
[49] Hengartner, N.W. (1996). Rate optimal estimation of additive regression via
the integration method in the presence of many covariates. Preprint, Department of Statistics, Yale University.
[50] Cook and Weisberg (1994). An Introduction to Regression Graphics. Wiley,
New York.

Regresion No Parametrica en R

Cargado por

Copyright:

Formatos disponibles

Regresion No Parametrica en R

Cargado por

Información del documento

Título original

Derechos de autor

Formatos disponibles

Compartir este documento

Compartir o incrustar documentos

Opciones para compartir

¿Le pareció útil este documento?

¿Este contenido es inapropiado?

Copyright:

Formatos disponibles

Regresion No Parametrica en R

Cargado por

Copyright:

Formatos disponibles

REGRESION

Autora: Nisa Boukichou Abdelkader

sus deseables caractersticas de interpretabilidad y sencillez en la practica.

realizar conclusiones acerca de dichos procedimientos.

Consideremos en primer lugar el planteamiento de un problema de regresion,

3. Estimar las derivadas de m, as como algunas propiedades relacionadas.

La teora y los metodos de suavizamiento o regresion no parametrica han cobrado

DEL MODELO DE REGRESION

metodologas ha tenido dos razones principales: la primera, que los estadsticos se

Figura 1.1: Ancho de Banda

(N adaraya W atson, 1964)

3. Suavizamiento por Splines: se define como la solucion a un problema de

(m(r) (x))2 dx ( > 0)

4. Estimadores basados en Desarrollos en Serie Ortogonal o los que estan

DEL MODELO DE REGRESION

Ruppert, Wand y Carroll (2003) establecen los siguientes factores a tener en

Si suponemos que la funcion de regresion,m tiene p derivadas en un punto x0 ,

As, se obtienen estimaciones de los coeficientes bj con j = 0, . . . , p y entonces,

h es un parametro denominado ancho de banda o par

utilizados por teoricos y analistas de datos, es el ajuste lineal el que ha mostrado

4. Gaussiano: k(u) = (2)1 exp( x2 )1|u|1

- El Error Cuadratico Medio Condicional en su version global (MISE):

Por tanto, como sabemos el error cuadratico medio se puede descomponer en

En cuanto, a la eleccion entre un criterio u otro dependera de que busquemos

La validacion cruzada por mnimos cuadrados.

Consideramos ahora la extension del modelo (1.1) al caso en que se considera

(Xi xj )lj KH (Xi x)

El problema en forma matricial y su solucion, se definen mediante:

cp (x) = eT1 XTx Wx,H Xx

donde e1 es un vector (pD + 1), con un 1 en la primera posicion y 0 en el resto,

1 (X1 x)T ((X1 x)p )T

p = 1, es el estimador lineal local multivariante (Ruppert y Wand, 1994):

c1 (x) = eT1 XTx,H Wx,H Xx

Aunque la generalizacion al caso multidimensional de la mayora de tecnicas

Modelos aditivos no param

Propuestos por Hastie y Tibshirani (1990), son una tecnica de regresion no

parametrica multivariante, los modelos aditivos pueden interpretar facilmente los

m(x) = E Y |X = (x1 , . . . , xD )T = + m1 (x1 ) + . . . + mD (xD )

La seleccion de una matriz de anchos de banda en el caso general multivariante,

El metodo de validacion cruzada de este modelo se puede encontrar en el libro

El programa R es un entorno de analisis y programacion estadstico que forma

En este libro llamado KernSmooth se trabajara con funciones para el suavizado

2.2. LIBRO KERNSMOOTH

> est <- bkde(x, bandwidth=0.25)

Figura 2.1: KernSmooth de bkde

> est <- bkde2D(x, bandwidth=c(0.7, 7))

Figura 2.2: KernSmooth de bkde2D1

Figura 2.3: KernSmooth de bkde2D2

2.2. LIBRO KERNSMOOTH

Valor devuelto por la funci

> h <- dpih(x)

1.3598478 1.5686859 1.7775241 1.9863623

Figura 2.4: KernSmooth de dpih

2.2. LIBRO KERNSMOOTH

Figura 2.5: KernSmooth de dpik

2.2. LIBRO KERNSMOOTH

Valor devuelto por la funci

data(geyser, package = "MASS")

> fit <- locpoly(x, y, bandwidth = h)