Resumo

Dados envolvendo medidas angulares estão presentes nas mais diversas áreas do
conhecimento. Para analisá-los é necessário utilizar uma teoria estatística específica e
apropriada, diferente da que utilizamos para dados lineares. Particularmente, quando
o interesse for formular, ajustar e fazer diagnósticos em modelos de regressão, uma
vez que, neste contexto, a natureza da variável deve ser considerada. Neste traba-
lho, utilizamos os modelos de regressão von Mises para investigar a associação tipo
circular-linear e apresentamos dois resíduos padronizados que foram obtidos a partir
da componente da função desvio e cujas distribuições de probabilidades podem ser
aproximadas pela distribuição normal padrão, definida para dados lineares.
Abstract

Data involving angular are present in the most diverse areas of science. To analyze
them is necessary to introduce an appropriate theory and to study specific and appro-
priate statistics as well, different from that we use for linear data. When the interest
is to formulate, to adjust and to make diagnostics on regression models, the nature of
the variables must be considered. In this work, we use the von Mises regression models
to investigate the circular-linear association and discuss two standardized residuals de-
fined from the component of the deviance function whose probability distributions can
be approximated by the normal standard distribution defined for linear data.
 Universidade Federal de Campina Grande
Centro de Ciências e Teconologia
Programa de Pós-Graduação em Matemática
Curso de Mestrado em Matemática

Análise de Resíduos em
Modelos de Regressão von Mises

por

Grayci-Mary Gonçalves Leal

sob orientação do

Prof. Dr. Francisco Antônio Morais de Souza

Dissertação apresentada ao Corpo Docente do Programa

de Pós-Graduação em Matemática - CCT - UFCG, como
requisito parcial para obtenção do título de Mestre em
Matemática.

Campina Grande - PB
Abril/2006
 Análise de Resíduos em Modelos de
Regressão von Mises
por

Grayci-Mary Gonçalves Leal

Dissertação apresentada ao Corpo Docente do Programa de Pós-Graduação em

Matemática - CCT - UFCG, como requisito parcial para obtenção do título de Mestre
em Matemática.

Área de Concentração: Matemática

Linha de Pesquisa: Probabilidade e Estatística

Aprovada por:

Prof. Dr. André Gustavo Campos Pereira

Prof. Dr. Antonio José da Silva

Prof. Dr. Francisco Antônio Morais de Souza

Orientador

Universidade Federal de Campina Grande

Centro de Ciências e Tecnologia
Programa de Pós-Graduação em Matemática
Curso de Mestrado em Matemática

Abril/2006

ii
 Agradecimentos

Primeiramente, ao meu Senhor e Salvador Jesus Cristo que me concede mais esta
vitória, pois sem Ele nada posso fazer;

Ao meu esposo Janilson por seu amor, incentivo e compreensão sempre dispensado
e ao meu primeiro filho Carlos Alberto, por está presente neste momento tornando-o
mais especial;

Aos meus pais, Afonso e Fátima, por todo amor, confiança e apoio em todas as
áreas, não só agora, mas sempre. Aos meus irmãos Kelly e Affonso por sempre acredi-
tarem e confiarem em mim, enfim agradeço a toda minha família;

Ao professor Francisco Morais, pela dedicação, humildade e paciência ao me ori-

entar em mais este trabalho;

Aos professores André Gustavo Campos Pereira e Antonio José da Silva, por gen-
tilmente aceitarem ao convite de avaliar este trabalho;

A todos os meus professores de graduação e pós-graduação da UAME/UFCG,

entre eles Aparecido Jesuíno, Antonio José, Daniel Pelegrino e Marco aurélio. E aos
funcionários Salete, D. Argentina, Valdir e Marcelino;

Aos meu colegas Lya, Jaqueline, Areli, Tatiana, Ana Cristina, Rosângela, Marta,
Jesualdo, Lino, Aluízio, Cícero, Luciano, Moisés, Orlando e Lauriclécio, pela amizade
e companherismo de todos;

À CAPES, pelo apoio financeiro.

iii
Dedicatória

Ao meu esposo Janilson e aos

meus pais Afonso e Fátima.

iv
Conteúdo

Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

1 Dados Direcionais 6
1.1 Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.2 Representação Gráfica de Dados Direcionais . . . . . . . . . . . . . . . 7
1.3 Medidas Descritivas no Círculo . . . . . . . . . . . . . . . . . . . . . . 9
1.3.1 Medidas de Locação . . . . . . . . . . . . . . . . . . . . . . . . 10
1.3.2 Medidas de Concentração . . . . . . . . . . . . . . . . . . . . . 11
1.3.3 Medidas de Dispersão . . . . . . . . . . . . . . . . . . . . . . . . 12
1.3.4 Outras Medidas de Dispersão . . . . . . . . . . . . . . . . . . . 12
1.4 Momentos Trigonométricos . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.5 Conceitos Teóricos Básicos para Dados Direcionais . . . . . . . . . . . . 15

2 Modelos Probabilísticos para Dados Circulares 21

2.1 Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.2 Modelo Exponencial . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.3 Modelos de Transformação . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.4 A Distribuição Circular Uniforme - U . . . . . . . . . . . . . . . . . . . 23
2.5 As Distribuições Arqueadas . . . . . . . . . . . . . . . . . . . . . . . . 24
2.5.1 A Distribuição Normal Arqueada - WN (µ, σ 2 ) . . . . . . . . . 25
2.5.2 A Distribuição Poisson Arqueada - WP(λ) . . . . . . . . . . . 27
2.5.3 A Distribuição Cauchy Arqueada - WC(µ, ρ) . . . . . . . . . . 28
2.6 A Distribuição von Mises - VM(µ, λ) . . . . . . . . . . . . . . . . . . 29
2.6.1 Estimação dos parâmetros da Distribuição von Mises . . . . . . 31
 2

2.7 Outras Distribuições Circulares . . . . . . . . . . . . . . . . . . . . . . 33

3 Modelos de Regressão para uma Resposta Angular 34

3.1 Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.2 Alguns Modelos para Regressão Angular . . . . . . . . . . . . . . . . . 35
3.3 Os Modelos de Regressão von Mises . . . . . . . . . . . . . . . . . . . . 36
3.3.1 Modelo de Médias von Mises . . . . . . . . . . . . . . . . . . . 37
3.3.2 Modelo de Dispersão von Mises . . . . . . . . . . . . . . . . . . 41
3.3.3 Modelo Misto von Mises . . . . . . . . . . . . . . . . . . . . . . 44

4 Resíduos para uma Resposta Angular 46

4.1 Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
4.2 A componente da função desvio . . . . . . . . . . . . . . . . . . . . . . 47
4.3 Resíduos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
4.3.1 Resíduo d∗i . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
4.3.2 Resíduo ri . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

5 Aplicação 55
5.1 Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
5.2 Descrição e Análise dos Dados . . . . . . . . . . . . . . . . . . . . . . 56
Conclusões . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

A Demonstração das Equações (1.5) e (1.6) 61

A.1 Demonstração . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

B Alguns Resultados sobre Dados Direcionais 63

B.1 Demonstração das Propriedades da Função Característica φY (p) . . . . 63
B.2 Demonstração para a Desigualdade tipo Tchebychev . . . . . . . . . . . 65

C Demonstração da Equação (4.10) 66

C.1 Enunciado do Teorema da Transformação Integral . . . . . . . . . . . 66
C.2 Demonstração da Equação (4.10) . . . . . . . . . . . . . . . . . . . . . 66

D Alguns resultados sobre as funções de Bessel 68

D.1 Derivadas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
D.2 Fórmulas de recorrência . . . . . . . . . . . . . . . . . . . . . . . . . . 69
 1

E Programas Computacionais 70
E.1 Comandos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70

Bibliografia 78
Introdução

Os modelos de regressão fornecem meios capazes de realizar previsões dos valores

da variável resposta, através de informações sobre as variáveis explicativas, ou seja,
eles nos permite investigar possíveis associações entre estas variáveis. No entanto, na
prática, nem sempre encontramos fenômenos que nos forneçam apenas dados lineares,
pois existem muitos exemplos, nas mais diversas áreas de conhecimento, em que os
dados a serem modelados são também direcionais, como no Meio Ambiente, na Mete-
orologia, na Geologia, entre outras. Por isso, faz-se necessário o desenvolvimento de
uma teoria apropriada para dados como esses. Portanto, quando se investiga associação
envolvendo variáveis direcionais e variáveis lineares, é necessário levar em consideração
os seguintes tipos de associações:

• circular-circular que mede o grau de associação entre duas variáveis direcionais;

• circular-linear que pode ser usada para predizer valores de uma variável dire-
cional a partir de valores de variáveis lineares;

• linear-circular que pode ser usada para predizer valores de uma variável linear
dado o valor de uma variável direcional.

Após a formulação do modelo, uma etapa de grande importância na análise de um

ajuste de regressão é a aplicação das técnicas de diagnósticos, dentre elas a análise de
resíduos que tem como objetivo detectar a presença de pontos extremos e avaliar a
adequação da distribuição proposta para a variável resposta. Em se tratando de dados
direcionais, faz-se necessário a utilização de resíduos convenientes, que destaquem as
características desses dados.
 5

Neste trabalho, consideramos a associação do tipo circular-linear, por isso foi conve-
niente considerar os modelos de regressão von Mises, por eles gozarem de boas pro-
priedades. Apresentamos dois resíduos padronizados que foram obtidos a partir da
componente da função desvio e cujas distribuições de probabilidades podem ser com-
paradas com a distribuição normal padrão.

A dissertação encontra-se estruturada da seguinte forma.

No Capítulo 1, apresentamos os principais conceitos relacionados com os dados
direcionais e suas estatísticas resumo.
No Capítulo 2, mostramos as principais famílias dos modelos probabilísticos para os
dados direcionais com as principais distribuições que ocorrem sobre o círculo.
No Capítulo 3, discutimos os modelos de regressão von Mises.
No Capítulo 4, apresentamos duas caracterizações para os resíduos de dados direcionais.
No Capítulo 5, fazemos uma aplicação a respeito da associação entre direções de ventos
e concentração de ozônio na atmosfera.
Capítulo 1

Dados Direcionais

Neste capítulo introduzimos os conceitos relacionados com os dados direcionais,

bem como apresentar as principais estatísticas resumo a esse tipo de observação.

1.1 Introdução
Em diversas áreas de conhecimento há o interesse por medidas angulares. Po-
demos citar como exemplo a Geologia, Biologia, Meteorologia, Medicina, Física, En-
genharia, Meio Ambiente, entre outras. Alguns exemplos típicos de dados direcionais
nestas áreas são:

(i) Direção de ventos e correntes marinhas;

(ii) Direção de migrações de animais;

(iii) Orientação da fratura de uma rocha;

(iv) Estudos de mapas mentais;

(v) Transporte de materiais geológicos.

Observações considerando direções podem ser interpretadas como dados direcio-

nais. Cada observação direcional pode ser representada, geometricamente, como um
ponto no círculo unitário, no caso bidimensional, ou como um ponto sobre a superfície
de uma esfera, no caso tridimensional, desde que a orientação e a direção inicial sejam
 7

definidas. Na verdade, os dados direcionais podem ser tratados em um espaço de qual-

quer dimensão. No entanto, neste trabalho, nosso interesse será trabalhar no espaço
de dimensão 2.
Os dados direcionais podem ser expressos em qualquer medida de ângulo, como
por exemplo: graus, grado, radiano, etc. Mas, por convenção, trabalha-se com os dados
em radianos, com intervalos de variação [0, 2π) ou (−π, π], dependendo do interesse da
análise.
Uma observação de uma variável circular é representada por um ponto P sobre
a circunferência de um círculo unitário centrado na origem do sistema de coordenadas
−→
cartesianas e, portanto, tem-se que OP é um vetor unitário em R2 . Podemos ver esta
−→
observação ainda como um ângulo y, formado pelo vetor OP e o semi-eixo positivo Ox ,
medido no sentido anti-horário. Sendo assim, uma observação de uma variável circular
pode ser representada pelo par ordenado (cos(y), sen(y)), em coordenadas cartesianas.
Existem casos de observações em eixos, em que não é possível distinguir qual é o início
ou o fim da observação, como por exemplo: o estudo de fraturas em rochas. Esse tipo
de dado é denominado axial, diferentemente do anterior, denominado vetorial (Fisher,
1995). Para se trabalhar com dados axiais deve-se, inicialmente, convertê-los a dados
vetoriais, transformando a observação y em 2y, ou seja, duplicando o ângulo (Mardia
e Jupp, 1999).

1.2 Representação Gráfica de Dados Direcionais

Os dados direcionais, assim como os dados lineares, podem ser representados por
diagramas e gráficos, com a finalidade de facilitar sua interpretação e análise. Quando
os dados não são agrupados, sua representação consiste em corresponder a cada obser-
vação bruta um ponto no círculo unitário. No entanto, quando os dados são agrupados
eles podem ser representados por uma das três maneiras seguintes: histograma circular,
histograma linear ou diagrama de rosas.
Os histogramas circulares são formados por um conjunto de barras que represen-
tam as freqüências de cada classe, cujas bases estão sobre a circunferência do círculo
unitário dividido de acordo com o tamanho de cada classe, com e áreas proporcionais
às freqüências. Se fizermos um corte no eixo deste histograma e o transportá-lo ao
 8

eixo de coordenadas X iremos transformá-lo num histograma linear, semelhante ao de

dados lineares.
O diagrama de rosas é um tipo de histograma circular formado por setores circu-
lares cujas áreas são proporcionais às freqüências de cada classe e os raios são propor-
cionais à raiz quadrada de cada freqüência.
Exemplo: Consideremos o conjunto de dados representando 19 observações de
direções de ventos, expresso em graus, da tabela abaixo. Representaremos os dados
não-agrupados, e em seguida os agrupado em histograma circular e diagrama de rosas.

Tabela 1.1: Medidas de Direção de Ventos

327 8 91 204 88 86 305 333 344 18
Direção
270 57 67 6 21 11 281 27 84

Figura 1.1: Representação Gráfica dos Dados Não-Agrupados

 9

Figura 1.2: Representação Gráfica em Histograma Circular

Figura 1.3: Representação Gráfica em Diagrama de Rosas

1.3 Medidas Descritivas no Círculo

Depois que os dados direcionais são coletados e representados graficamente, preci-
saremos estabelecer algumas estatísticas descritivas adequadas. A princípio, é tentador
querermos utilizar as mesmas medidas empregadas para os dados lineares. Entretanto,
pelo exemplo fictício, a seguir, podemos verificar a inadequação destas medidas no
círculo.
Exemplo Fictício: Imaginemos uma amostra de tamanho 2 em que os ângulos
 10

observados tenham sido y1 = 10 e y2 = 3590 que são “próximos entre si” e eqüidistastes
da direção 00 . Neste caso, a média aritmética amostral é y = (10 + 3590 )/2 = 1800
que corresponde ao ponto no círculo diametralmente oposto ao resultado intuitivo. Em
contra partida, o desvio padrão amostral fica dado por: Sn−1 = 2530 , ou seja, um valor
que não descreve o fato de y1 e y2 estarem “próximos entre si”.
Como as medidas resumo usadas para variáveis lineares não fazem sentido para
as variáveis circulares, pelo que vimos no exemplo anterior, precisamos estabelecer
algumas medidas descritivas relacionadas com o conceito de posição, dispersão e con-
centração que façam sentido para as variáveis circulares.

1.3.1 Medidas de Locação

1.3.1.1 Média Direcional

Considere uma amostra aleatória de n observações circulares y1 , y2 , . . . , yn , às quais

−→ −→ −→
podemos associar n vetores unitários OP 1 , OP 2 , . . . , OP n , conforme mencionamos an-
teriormente. Em termos algébricos, calculemos as quantidades
n
X n
X
S= sen(yi ), C = cos(yi ) e R2 = C 2 + S 2 .
i=1 i=1

A média direcional de y1 , y2 , . . . , yn é definida como sendo o ângulo µ, correspondente

−→ −→ −→
ao vetor resultante da soma OP 1 + OP 2 + · · · + OP n . Esse ângulo deve satisfazer as
condições cos(µ) = C/R e sen(µ) = S/R, ou seja, µ = arctan(S/C) (Fisher, 1995).
Assim,


 arctan(S/C), se S ≥ 0 e C > 0


µ= arctan(S/C) + π, se C < 0 (1.1)



 arctan(S/C) + 2π, se S < 0 e C > 0·
√
A quantidade R = C 2 + S 2 representa o comprimento do vetor resultante e fornece
informação sobre o grau de concentração dos ângulos observados. Em vez de R, é mais
comum usar o comprimento médio do vetor resultante, definido por R = R/n, o qual
tem a vantagem de variar no intervalo [0,1]. Portanto, µ pode ser determinado também
como soluções das equações
C = R cos µ (1.2)

S = Rsenµ (1.3)
 11

onde
S C
S= e C= . (1.4)
n n
Das expressões (1.2), (1.3) e (1.4), podemos mostrar (vide Apêndice A) que
n
1X
cos(yi − µ) = R e (1.5)
n i=1
n
1X
sen(yi − µ) = 0. (1.6)
n i=1

Para dados axiais, a média axial é definida da seguinte forma: Sejam y1 , y2 , . . . , yn

ângulos que representam observações de dados axiais e seja µ a média direcional dos
µ µ
ângulos vetoriais correspondentes. Então a média axial é dada por 2
e 2
+ π (Mardia
e Jupp, 1999).

1.3.1.1 Mediana Direcional

Definição: A mediana direcional Ye dos ângulos y1 , y2 , . . . , yn é um ângulo φ que

satisfaça as seguintes condições:

(i) A metade dos pontos situam-se no arco [φ, φ + π);

(ii) A maioria dos pontos são mais próximos de φ do que φ + π.

Para n ímpar, a mediana é um único ângulo e, quando n for par a mediana é a média
aritmética entre os dois ângulos centrais.

1.3.2 Medidas de Concentração

1.3.2.1 O Comprimento Médio do Vetor Resultante

Como definido anteriormente, R = R/n onde R é o comprimento do vetor resultante

−→ −→ −→ −→
OP 1 + OP 2 + · · · + OP n . Como cada OP i é unitário, temos

−→ −→ −→ −→ −→
0 ≤k OP 1 + OP 2 + · · · + OP n k≤k OP 1 k + · · · + k OP n k= 1 + · · · + 1 = n.

Portanto,
0 ≤ R ≤ n ⇒ 0 ≤ R/n ≤ 1 ⇒ 0 ≤ R̄ ≤ 1.

A quantidade R̄ fornece informações sobre o grau de concentração dos ângulos ob-

servados. No caso em que todas as observações forem iguais a um mesmo ângulo,
 12

digamos y, tem-se que R̄ = 1, representando uma concentração máxima dos dados, e

corresponde a uma variabilidade mínima entre eles. Outro caso extremo é quando os
dados encontram-se uniformemente distribuídos no círculo, nesta situação tem-se uma
concentração mínima (R̄ = 0) e uma variabilidade máxima. É importante ressaltarmos
que R̄ = 0 não significa uma distribuição uniforme no círculo, pode acontecer algum
tipo de estrutura de agrupamento, como o exemplo citado em Fisher (1995, p.32).

1.3.3 Medidas de Dispersão

1.3.3.1 A Variância Circular

A variância circular amostral é definida por

V = 1 − R. (1.7)

O conceito é análogo à variância de dados lineares, no sentido de que quanto menor o

valor da variância circular, maior a concentração nos dados. Observe que

0 ≤ V ≤ 1, uma vez que 0 ≤ R̄ ≤ 1.

1.3.3.2 O Desvio Padrão Circular

Para dados direcionais, a definição de desvio padrão é diferente da definição para

dados lineares, que corresponde à raiz quadrada da variância. O desvio padrão circular
amostral é definido por

σ = {−2 log(1 − V )}1/2 = {−2 log R̄}1/2 .

Uma boa aproximação para o desvio padrão circular amostral é


 (2V )1/2 , para V pequeno (V < 0, 18)
σ=
 {2(1 − R)}1/2 , para R grande (R > 0, 82)

1.3.4 Outras Medidas de Dispersão

1.3.4.1 A Dispersão Circular Amostral

A dispersão circular amostral é definida por

1 − R2
δb = 2 , (1.8)
2R
onde R2 é o comprimento médio resultante dos ângulos duplicados 2y1 , . . . , 2yn . Com
q
2 2 P P
R2 = C 2 + S 2 onde C 2 = n1 ni=1 cos(2yi ) e S 2 = n1 ni=1 sen(2yi ).
 13

1.3.4.2 O Desvio Médio Circular

Uma alternativa para medir a distância entre dois ângulos α e β é


 |α − β|, se |α − β| ≤ π
min((α − β), 2π − (α − β)) =
 2π − |α − β|, se |α − β| > π

Desta forma, podemos escrever


 π − π + |α − β|, se π − |α − β| ≥ 0
min((α − β), 2π − (α − β)) =
 π + π − |α − β|, se π − |α − β| < 0

Então,

 π − (π − |α − β|), se π − |α − β| ≥ 0
min((α − β), 2π − (α − β)) =
 π − [−(π − |α − β|)], se π − |α − β| < 0

Portanto,
min((α − β), 2π − (α − β)) = π − |π − |α − β||.

Com isto, uma medida de dispersão dos ângulos y1 , . . . , yn em torno de um dado ângulo
αé
n n
1X 1X
d0 (α) = {π − |π − |yi − α||} = π − |π − |yi − α||
n i=1 n i=1
A mediana amostral ye é o ponto que minimiza d0 (Mardia e Jupp, 1999), logo, o desvio
médio circular é definido por
n
1X
d0 (e
y) = {π − |π − |yi − ye||}. (1.9)
n i=1

1.3.4.3 A Diferença Circular Média

A diferença circular média é a distância média entre observações circulares é dada por
n
1 1X
D0 = [d0 (y1 ) + d0 (y2 ) + · · · + d0 (yn )] = d0 (yj ).
n n j=1

Portanto,
n n
1 XX
D0 = 2 {π − |π − |yi − yj ||}, (1.10)
n i=1 j=1

que representa a distância média entre pares de observações.

 14

1.3.4.4 A Amplitude Circular

A amplitude circular é o comprimento do menor arco que contém todas as observa-

ções. Para determinarmos a distância circular inicialmente consideraremos os ângulos
y1 , . . . , yn de forma ordenada y(1) ≤ · · · ≤ y(n) . Os comprimentos de arco entre obser-
vações adjacentes são

 y
(i+1) − y(i) , i = 1, . . . , n − 1
Ti =
 2π − y + y , i = n
(n) (1)

Assim, definimos a amplitude circular w por:

w = 2π − max(T1 , . . . , Tn ) (1.11)

1.4 Momentos Trigonométricos

As expressões S e C definidos em (1.4), são importantes para determinação da
média direcional e variância amostral, como visto anteriormente. Definimos o primeiro
momento trigonométrico em torno de 00 por
′
m1 = C + iS. (1.12)

Observe que, de (1.2) e (1.3)

′
m1 = R cos µ + iRsenµ = R(cos µ + isenµ) = Reiµ .

Por extensão, definimos o p-ésimo momento trigonométrico em torno de 00 , como sendo

′
mp = ap + ibp , para p = 1, 2, . . . (1.13)

em que
n n
1X 1X
ap = cos(pyi ) e bp = sen(pyi ).
n i=1 n i=1
Portanto,
′
mp = Rp eiµp , (1.14)

onde µp e Rp denotam a média direcional e o comprimento médio do vetor resultante

amostral de py1 , . . . , pyn , respectivamente. o p-ésimo momento trigonométrico em torno
da média direcional µ é
mp = ap + ibp , (1.15)
 15

onde
n n
1X 1X
ap = cos[p(yj − µ)] e bp = sen[p(yj − µ)].
n j=1 n j=1
Em particular, temos
n n
1X 1X
m1 = a1 + ib1 = cos[p(yj − µ)] + i sen[p(yj − µ)].
n j=1 n j=1

Dos resultados, (1.5) e (1.6) temos,

m1 = R + i.0 = R.

Posteriormente, apresentaremos as versões populacionais dos momentos trigonométri-

cos, que desempenham importante papel na teoria de distribuições de probabilidade
sobre o círculo.

1.5 Conceitos Teóricos Básicos para Dados Direcio-

nais
1.5.1 A Função de Distribuição

Entendemos por uma variável aleatória circular como sendo uma variável aleatória já
convertida em ângulos e assumindo valores no intervalo [0, 2π).
Definição: Definimos a função de distribuição de uma variável aleatória circular como
sendo uma função F tal que

F (y) = P (0 < Y ≤ y), ∀ y ∈ [0, 2π)

(1.16)
F (y + 2π) − F (y) = 1, ∀ y ∈ R.

Sejam y1 e y2 satisfazendo −∞ < y1 < y2 < ∞, com y2 − y1 < 2π. A probabilidade de

Y pertencer ao arco (y1 , y2 ) é dada por
Z y2
P (y1 < Y ≤ y2 ) = F (y2 ) − F (y1 ) = dF (y), (1.17)
y1

onde a integral é a integral de Lebesgue-Stieltjes.

De fato, (1.17) é verdadeira, pois para y1 < y2 temos que,

P (Y ≤ y2 ) = P (Y ≤ y1 ) + P (y1 < Y ≤ y2 )

⇒ P (y1 < Y ≤ y2 ) = P (Y ≤ y2 ) − P (Y ≤ y1 ).
 16

Logo,
P (y1 < Y ≤ y2 ) = F (y2 ) − F (y1 ).

Observação: Embora F (0) = 0 e F (2π) = 1, para as variáveis aleatórias circulares, a

função de distribuição no círculo não satisfaz as condições lim F (x) = 0 e lim F (x) =
x→−∞ x→∞
1, que são satisfeitas por funções de distribuições de variáveis aleatórias lineares.
Definição: Seja F a função de distribuição de uma variável aleatória circular Y . Se F
for absolutamente contínua, definimos a função densidade de probabilidade de Y como
sendo a função f tal que
Z y2
f (y)dy = F (y2 ) − F (y1 ), −∞ < y1 ≤ y2 < ∞. (1.18)
y1

Observação: Uma função f é função de densidade de probabilidade de uma distribui-

ção absolutamente contínua se, e somente se:

(i) f (y) ≥ 0, quase sempre em (−∞, ∞);

(ii) f (y + 2π) = f (y), quase sempre em (−∞, ∞);

R 2π
(iii) 0
f (y)dy = 1.

Definição: Sejam Y uma variável aleatória circular. Os momentos circulares de Y são

dados por
R 2π
αp = E[cos(pY )] = 0
cos(py)dF (y) (1.19)
R 2π
βp = E[sen(pY )] = 0
sen(py)dF (y) (1.20)

onde p ∈ Z. Como conseqüência das propriedades das funções trigonométricas tem-se

que

αp = α−p , βp = −βp , |αp | ≤ 1 e |βp | ≤ 1.

1.5.2 A Função Característica

Definição: Definimos a função característica de uma variável aleatória circular Y como

sendo o conjunto de números complexos {φY (p) : p ∈ R} dada por
Z 2π
ipY
φY (p) = E(e ) = eipY dF (y), com p ∈ R. (1.21)
0
 17

Podemos escrever (1.21) da forma,

φY (p) = αp + iβp ,

onde αp e βp são os momentos circulares de Y definidos em (1.19) e (1.20). Como

está mostrado no Apêndice B, a função característica goza das seguintes propriedades:
Propriedades:

(i) φY (0) = 1

(ii) φY (p) = φ−Y (p) onde φY (p) é o conjugado complexo de φY (p)

(iii) |φY (p)| ≤ 1,

(iv) É suficiente considerarmos a função φY (p) apenas para valores inteiros e positivos
de p.

1.5.3 Independência e Convolução

Definição: Sejam Y1 e Y2 duas variáveis aleatórias circulares. A função característica

conjunta (Y1 , Y2 ) é definida pela seqüencia dupla {φY1 Y2 (p, q) com p, q ∈ Z} dada por

φY1 Y2 (p, q) = E(eipY1 +iqY2 ). (1.22)

Definição: As variáveis Y1 e Y2 são independentes se, e somente se, φY1 Y2 (p, q) =

φY1 (p)φY2 (q) para todo p, q ∈ Z.
Sejam Y1 , . . . , Yn variáveis aleatórias circulares independentes. A função característica
P
da soma Sn = ni=1 Yi é análoga ao caso linear. Se as variáveis forem identicamente
distribuídas temos que
n
Y
φSn (p) = φYi (p) = φY1 (p). . . . .φYn (p) = [φY (p)]n .
i=1

Definição: Sejam Y1 e Y2 variáveis aleatórias circulares independentes com respectivas

funções de distribuição F1 e F2 e considere S = Y1 + Y2 . A função de distribuição F
de S é dada por Z 2π
dFS (s) = dF2 (s − y1 )dF1 (y1 ).
0
 18

1.5.4 Momentos e Medidas de Locação e Dispersão

Vamos considerar as versões populacionais das medidas resumo e dos momentos trigo-
nométricos que foram apresentadas nas seções (1.3) e (1.4).

1.5.4.1 Momentos Trigonométricos

Sejam αp = E[cos(pY )] e βp = E[sen(pY )], com p ∈ Z, os momentos trigonométricos de

uma variável aleatória circular Y . A seqüencia {(αp , βp ), p ∈ Z} caracteriza a função
característica de Y . Portanto, no círculo, qualquer distribuição é determinada pelos
seus momentos trigonométricos.
Para p ≥ 0 temos que a versão populacional do p-ésimo momento em torno de 00 é
dada por
φY (p) = ρp eiµp , ρp ≥ 0 (1.23)

onde ρp é o comprimento médio do vetor resultante. O p-ésimo momento trigonométrico

populacional em torno da média direcional é definido por

φY (p) = αp + iβ p , (1.24)

onde, αp = E{cos[p(Y − µ)]} e β p = E{sen[p(Y − µ)]}. Logo, a equação (1.24) pode

ser escrito como

φY (p) = E[eip(Y −µ) ].

1.5.4.2 Medidas de Locação e Dispersão

1.5.4.2.1 A Média Direcional e o Comprimento Médio Resultante

Para p = 1, em (1.23), temos

φY (1) = ρ1 eiµ1 = ρeiµ ,

onde µ é a versão populacional da média direcional µ e ρ é a versão populacional do

comprimento médio resultante R e pode ser calculado pela expressão ρ = (α12 + β12 )1/2 .
Observe que, ao considerarmos o efeito de rotação por −ψ na variável circular Y, temos
Y ∗ = Y − ψ. A média direcional de Y ∗ é µ∗ = µ − ψ, que é equivariante sob rotação.
Então
φ∗Y (p) = e−ipψ φY (p). (1.25)
 19

De fato,

∗
φ∗Y (p) = E(eipY ) = E(eip(Y −ψ) ) = e−ipψ E(eipY ) = e−ipψ φY (p).

Considere p = 1 em (1.25), então

φ∗Y (1) = e−iψ φY (1) = e−iµ ρeiµ = ρ + i0.

Mas,
φ∗Y (1) = E[cos(Y − µ)] + iE[sen(Y − µ)].

Logo,
E[cos(Y − µ)] + iE[sen(Y − µ)] = ρ + i0.

Portanto,

E[cos(Y − µ)] = ρ, e (1.26)

E[sen(Y − µ)] = 0. (1.27)

1.5.4.2.2 A Variância Circular

Definição: A variância circular ν, de uma variável aleatória circular Y é definida como

ν = 1 − ρ. (1.28)

De (1.26) temos,
ν = 1 − E[cos(Y − µ)].

E, analogamente à situação amostral, 0 ≤ ν ≤ 1.

1.5.4.2.3 O Desvio Padrão Circular

A versão populacional do desvio padrão circular é dada por

σ = {−2 log(1 − ν)}1/2 = {−2 log ρ}1/2 . (1.29)

Esta definição é motivada pelo fato de que a distribuição normal ao redor do círculo
gera a distribuição normal arqueada, como veremos posteriormente.
 20

1.5.4.2.4 A Dispersão Circular

A versão populacional da dispersão circular é

1 − ρ2
δ= , (1.30)
2ρ22
p
onde ρ2 = α22 + β22 , com α2 = E[cos(2Y )] e β2 = E[sen(2Y )].

1.5.4.2.5 A Mediana Circular

A versão populacional da mediana circular Ye é a mediana direcional µ

e definida como
a direção φ que minimiza E[π − |π − |Y − φ|] e satisfaz

1 1
P (Y ∈ (e e + π]) ≥
µ, µ e P (Y ∈ (e
µ − π, µ
e)) ≥ .
2 2

1.5.5 A Desigualdade Tipo Tchebyshev

Várias desigualdades para distribuições são observadas no círculo, muitas delas foram
estudadas por Marshal & Olkin (1961). Usando a desigualdade de Tchebyshev na
variável aleatória sen( Y −µ
2
), temos:

Y −µ ν
P (|sen( )| ≥ ε) ≤ 2 , 0 < ε ≤ 1. (1.31)
2 2ε

A demonstração encontra-se no Apêndice B.

Capítulo 2

Modelos Probabilísticos para Dados

Circulares

Neste capítulo, nosso objetivo é fazer um estudo das duas principais famílias de
modelos probabilísticos para dados direcionais. Os chamados modelos exponenciais e
os modelos de transformação. Além disso, apresentar algumas das mais importantes
distribuições sobre o círculo, dando um destaque à distribuição von Mises.

2.1 Introdução
Em inferência estatística, a análise de modelos probabilísticos adequados é de
grande relevância, tanto para dados lineares quanto para dados direcionais. A maio-
ria dos modelos para dados direcionais pertence a uma das duas classes de modelos
paramétricos, ou os modelos exponenciais ou os modelos de transformação (Mardia e
Jupp, 1999). Assim sendo, faremos um estudo das principais propriedades dessas duas
classes de modelos. Dentre as famílias de distribuições no círculo, temos a distribuição
uniforme circular como a mais simples e a distribuição von Mises como a mais impor-
tante, pois desempenha um papel em inferência no círculo semelhante à distribuição
normal na reta. Apresentaremos também, outras distribuições circulares importantes
que são as distribuições arqueada.
 22

2.2 Modelo Exponencial

Definição: Um modelo é do tipo exponencial se sua função densidade de proba-
bilidade for da forma

′
f (x; ω) = b(x) exp{φ(ω) t(x) − ψ(ω)}, x ∈ χ e ω ∈ Ω, (2.1)

com relação a algum domínio de medida λ, onde χ é o espaço amostral e Ω é o espaço

paramétrico. A estatística t : χ → Rm é a estatística canônica e φ(ω) é o parâmetro
canônico (φ : Ω → Rm ).

Observação: Segundo Mardia e Jupp (1999), considere d = dim[φ(Ω)].

(i) Se (2.1) for minimal, então o modelo é chamado um modelo exponencial (m, d);

(ii) Se d < m, o modelo é chamado modelo exponencial curvo;

(iii) Se d = m podemos fazer uma reparametrização para o modelo exponencial (2.1).

Considere θi = φi (ω) para i = 1, . . . , m, então (2.1) pode ser escrita na forma

′
f (x; θ) = b(x)exp{θ t(x) − ψ(θ)}, (2.2)

onde θ varia em um subconjunto Θ ⊂ Rm . Quando Θ for um aberto em Rm , e ψ

pertencer ao conjunto das funções integráveis, então o modelo é chamado de Modelo
Exponencial Regular.

Principais Propriedades dos Modelos Exponenciais Regulares

(1) Os primeiros momentos da estatística canônica são dados por:

∂ψ ∂ 2ψ
Eθ (t) = e V arθ (t) = ′
∂θ ∂θ∂θ
(2) A matriz inversa de Informação de Fisher é V arθ (t).
(3) Se X1 , . . . , Xn forem independentes e identicamente distribuídas com função de
densidade de probabilidade dada por (2.2), então a distribuição amostral de t (média
amostral da estatística t ) tem função densidade de probabilidade proporcional a

′
exp{n[θ t(x) − ψ(θ)]},
 23

e t é suficiente para θ.
(4) O estimador de máxima verossimilhança θb de θ é único e é dado pela solução da
equação
Eθb (t) = t.

2.3 Modelos de Transformação

Definição: O modelo de transformação é um modelo em que um grupo G atua
tanto sobre o espaço amostral χ quanto sobre o espaço paramétrico Ω, de modo que, se
a variável aleatória X tem função densidade de probabilidade f (x, ω), então a variável
aleatória gX ∈ G tem função densidade de probabilidade f (gx, gω) em que

f (gx, gω) = f (x, ω)χ(g, x), (2.3)

para algumas funções g sobre G × χ.

Em outras palavras, um modelo de transformação é um modelo composto de transfor-
mação em que o grupo G age transitivamente sobre o espaço paramétrico, isto é, para
′ ′
quaisquer ω e ω ∈ Ω, existe um g ∈ G tal que gω = ω . Em particular, quando o
grupo G é o grupo aditivo R e a medida é a medida de Lebesgue, então (2.3) fica da
forma
f (x + g, ω + g) = f (x, ω). (2.4)

Para nossos propósitos, a principal propriedade dos modelos compostos de transforma-

ção é que o estimador de máxima verossimilhança ω
b de ω é equivariante, isto é,

ω
b (gx1 , . . . , gxn ) = gb
ω (x1 , . . . , xn ).

2.4 A Distribuição Circular Uniforme - U

Como mencionado anteriormente, a distribuição uniforme é a distribuição mais
simples sobre o círculo. Ela é a única que tem a propriedade de ser equivariante sob
rotação e reflexão (Mardia e Jupp, 1999).
Definição: Uma variável aleatória circular Y tem distribuição uniforme sobre o círculo
se sua função densidade de probabilidade for da forma
1
f (y) = , α ≤ y < α + 2π, com α ∈ [0, 2π). (2.5)
2π
 24

Observe que, se α ≤ β ≤ α + 2π, então

Z β
Pr (α < Y ≤ β) = F (β) − F (α) = f (y)dy.
α

Logo,
1 β β−α
Pr (α < Y ≤ β) = y|α = ,
2π 2π
ou seja, essa probabilidade é proporcional ao comprimento de arco.
A função característica da distribuição uniforme é dada por

 0, se p = ±1, ±2, . . .
φY (p) =
 1, se p = 0.

De fato, pela definição de função característica temos que

Z 2π Z 2π Z 2π
ipY ipY 1 1
φY (p) = e dF (y) = e dy = eipY dy.
0 0 2π 2π 0

Então, para p = 0, temos Z 2π

1
φY (p) = d(y) = 1.
2π 0

Para p = ±1, ±2, . . . , temos

1 1 ipy 2π e2ipπ − 1
φY (p) = e |0 = = 0.
2π ip 2ipπ

Portanto, 
 0, se p = ±1, ±2, . . .
φY (p) =
 1, se p=0

2.5 As Distribuições Arqueadas

Muitos modelos probabilísticos para dados direcionais são derivados a partir de
transformações de modelos probabilísticos lineares usuais. Em outras palavras, dada
uma distribuição sobre a reta podemos arqueá-la em torno da circunferência do círculo
de raio unitário e assim obter uma distribuição de probabilidade angular.
Definição: Seja X uma variável aleatória sobre a reta, definimos a variável aleatória
Xw da distribuição arqueada por

Xw = X{mod(2π)}. (2.6)
 25

Se consideramos o círculo como o conjunto de números complexos de módulo unitário,

a função de arqueamento X 7→ Xw pode ser escrita da seguinte forma

X 7→ e2πiX .

Se X tem função de distribuição FX , então a função de distribuição FY de Y é dada

por
∞
X
FY (y) = {F (y + 2kπ) − F (2kπ)}, 0 ≤ y < 2π. (2.7)
k=−∞

Se X tem uma função densidade de probabilidade f então a função densidade de

probabilidade correspondente fY de Y é
∞
X
fY (y) = f (y + 2kπ), 0 ≤ y ≤ 2π. (2.8)
k=−∞

Propriedades:

(i) (X + Z)ω = Xω + Zω

(ii) Se X tem função característica φ então a função característica φY (p) de Xω é:

φY (p) = φ(p)
P∞
(iii) Se φ é integrável então Y tem densidade com 2
p=1 (αp + βp2 ) < ∞, e
∞
X X∞
1
fw (y) = f (y + 2kπ) = [1 + 2 (αp cos py + βp senpy)]
k=−∞
2π p=1

2.5.1 A Distribuição Normal Arqueada - WN (µ, σ 2 )

Um caso especial dentre as distribuições arqueadas é aquela que se obtém a partir de

uma distribuição Gaussiana, por meio do seu arqueamento sobre o círculo.
Se uma variável aleatória X tem distribuição normal N (µ, σ 2 ), então a variável aleatória
Y = X{ mod (2π)} tem distribuição normal arqueada com parâmetros µ e σ 2 . Sua
função densidade de probabilidade é dada por
∞
X
2 1 (y + 2kπ − µ)2
fY (y; µ, σ ) = √ exp{− }· (2.9)
k=−∞
2πσ 2 2σ 2

De fato, se X ∼ N (µ, σ 2 ) sua função densidade de probabilidade é da forma

1 (x − µ)2
f (x) = √ exp{− }.
σ 2π 2σ 2
 26

Então, de (2.8) a função densidade de probabilidade de Y é dada por

∞
X ∞
X
2 1 (y + 2kπ − µ)2
fY (y; µ, σ ) = f (y + 2kπ) = √ exp{− }
k=−∞ k=−∞
2πσ 2 2σ 2

A função característica de uma variável aleatória X ∼ N (µ, σ 2 ) é φX (t) = exp{iµt −

t2 σ 2 /2}. Da propriedade (ii) da seção 2.5, segue-se que a função característica de Y é

φ(p) = exp{iµp − p2 σ 2 /2}

2 σ 2 /2
= eiµp e−p
2 /2)p2
= eiµp e(−σ
2
= eiµp ρp
2
= [cos(µp) + isen(µp)]ρp
2 2
= ρp cos(µp) + iρp sen(µp).

Portanto,
φ(p) = αp + iβp , (2.10)

onde

 α = ρp2 cos(µp)
p
(2.11)
 β = ρp2 sen(µp)·
p

Com este resultado podemos representar fω de outra maneira, a saber,

X∞
1 2
fw (y; µ, ρ) = {1 + 2 ρp cos p(y − µ)}. (2.12)
2π p=1

Para verificarmos isso, basta ver que, da propriedade (iii) da seção 2.5, temos que
∞
X X∞
1
fω (y) = f (y + 2kπ) = [1 + 2 (αp cos(py) + βp sen(py))].
k=−∞
2π p=1

Logo, de (2.11), temos que

2 2
αp cos(py) + βp sen(py) = ρp cos(µp) cos(py) + ρp sen(µp)sen(py)
2 σ 2 /2 2 σ 2 /2
= e−p cos(µp) cos(py) + e−p sen(µp)sen(py)
2 σ 2 /2
= e−p cos p(y − µ).

Logo,
X∞
1 2
fw (y; µ, ρ) = {1 + 2 ρp cos[p(y − µ)]}.
2π p=1

Conclusões:
 27

• Podemos determinar a função de distribuição Fω integrando, termo a termo, a

expressão (2.12).

• Quando σ 2 ≥ 2π, fω pode ser aproximada pelos três primeiros termos de (2.12).

• Quando σ 2 ≤ 2π, fω pode ser aproximada pelo termo k = 0 de (2.9).

Observações:

(i) A distribuição normal arqueada é unimodal e simétrica em torno da moda µ.

(ii) Quando ρ → 0 a WN (µ, σ 2 ) tende à distribuição uniforme.

(iii) Quando ρ → 1 a WN (µ, σ 2 ) tende à distribuição degenerada, concentrada em

µ.

2.5.2 A Distribuição Poisson Arqueada - WP(λ)

A redução módulo 2π “transforma” a reta num círculo. Da mesma forma, temos que se
m ∈ Z∗+ a redução módulo 2πm “transforma” os inteiros no grupo das m-ésimas raízes
de 1, visto como um subgrupo do círculo (Mardia e Jupp, 1999). Em particular, se X
é uma variável aleatória sobre os inteiros, então Xw , definida por

Xw = 2πX{mod(2π)}, (2.13)

é uma variável aleatória distribuída sobre o círculo. A função de probabilidade de Xw

é:
X∞
2πr
Pr (Xw = )= PX (r + km), r = 0, 1, . . . , m − 1, (2.14)
m k=−∞

Desta forma, se X ∼ Poisson (λ), Xω tem distribuição Poisson Arqueada com função
de probabilidade dada por
X∞
2πr −λ λr+km
Pr (Xω = )=e , r = 0, 1, . . . , m − 1, (2.15)
m k=0
(r + km)!

e denotamos por Xw ∼ WP(λ).

A função característica de Xw ∼ WP(λ) é dada por

φ(p) = exp{−λ(1 − e2πip/m )}, onde p = 0, 1, 2, . . . (2.16)

 28

Uma propriedade importante da distribuição Poisson Arqueada é que o conjunto das

suas funções densidades de probabilidades é fechado com relação à soma.
De fato, isto ocorre, pois se considerarmos Y1 ∼ WP (λ1 ) e Y2 ∼ WP (λ2 ) indepen-
dentes, então, a função característica de Y1 + Y2 é dada por

φY1 +Y2 (p) = φY1 (p)φY2 (p)

= exp{−λ1 (1 − e2πip/m )} exp{−λ2 (1 − e2πip/m )}

= exp{−(λ1 + λ2 )(1 − e2πip/m )}.

Portanto, Y1 + Y2 ∼ WP (λ1 + λ2 ).

2.5.3 A Distribuição Cauchy Arqueada - WC(µ, ρ)

Na reta, a distribuição de Cauchy tem função densidade de probabilidade dada por

1 a
f (x; µ, a) = , a > 0, −∞ < µ < ∞, x ∈ R.
π a2 + (x − µ)2

onde µ é a mediana (parâmetro de locação) e a representa a distância entre a mediana

e o primeiro quartil (parâmetro de escala).
A função característica da distribuição de Cauchy é

φ(t) = e−a|t|+itµ , onde t ∈ R

Daí, a função característica da distribuição de Cauchy Arqueada para p ∈ Z é da forma

φ(p) = e−a|p|+ipµ

= e−a|p| eipµ

= e−a|p| (cos pµ + isenpµ).

Logo,
φ(p) = e−a|p| (cos pµ) + e−a|p| (senpµ) = αp + iβp .

onde,

αp = e−a|p| cos(pµ) e
(2.17)
βp = e−a|p| sen(pµ).
 29

Assim, a função densidade de probabilidade da distribuição Cauchy Arqueada fica dada

por
X∞
1
c(y; µ, ρ) = {1 + 2 ρp cos p(y − µ)}. (2.18)
2π k=−∞

Da propriedade (iii) da seção (2.5) temos que,

∞
X X∞
1
c(y; µ, ρ) = f (y + 2kπ; µ, a) = [1 + 2 (αp cos py + βp senpy)].
k=−∞
2π p=1

Substituindo (2.17) na expressão acima, e considerando ρ = e−a , segue-se que

X∞
1
c(y; µ, ρ) = {1 + 2 e−a|p| cos p(y − µ)}
2π k=−∞
X∞
1
= {1 + 2 ρ|p| cos p(y − µ)}.
2π k=−∞

Propriedades:

(i) A média direcional de WC(µ, ρ) é µ{ mod (2π)} e o comprimento médio do vetor

resultante é ρ;

(ii) A distribuição WC(µ, ρ) é unimodal e simétrica em torno de µ;

(iii) Quando ρ → 0 a WC(µ, ρ) converge para a distribuição circular uniforme;

(iv) Quando ρ → 1 a WC(µ, ρ) fica concentrada no ponto µ;

(v) A função de distribuição da WC(µ, ρ) é dada por

½ ¾
1 −1 (1 + ρ2 ) cos(y − µ) − 2ρ
F (y) − F (µ) = cos ;
2π 1 + ρ2 − 2ρ cos(y − µ)

(vi) Tal como a Distribuição Poisson Arqueada, o conjunto das funções densidade de
probabilidades das distribuições de Cauchy Arqueadas é fechado com relação à
soma.

2.6 A Distribuição von Mises - VM(µ, λ)

Definição: Dizemos que uma variável aleatória circular Y tem distribuição von
Mises com parâmetros µ e λ se sua função densidade de probabilidade for dada por
1
fY (y; µ, λ) = eλ cos(y−µ) , 0 ≤ y < 2π, 0 ≤ µ < 2π, e λ > 0, (2.19)
2πI0 (λ)
 30

onde, I0 (λ) é a função de Bessel modificada de primeiro tipo e ordem zero, que é dada
pela série
X∞
1 1 2r
I0 (λ) = ( λ) ·
r=0
r!2 2

A função de distribuição de Y, é dada por

Z y
1
FY (y) = eλ cos(t−µ) dt, (2.20)
2πI0 (λ) 0

que não admite uma forma fechada. Algumas características da distribuição von Mises
Ip (λ)
envolvem a função Ap (λ) = I0 (λ)
, onde Ip (λ) é a função de Bessel modificada do primeiro
tipo e ordem p, que é definida pela série
∞
X 1 1
Ip (λ) = ( λ)2r+p , p = 0, 1, . . .
r=0
(p + r)!r! 2

As propriedades desta função podem ser encontradas em Abramowitz e Stegun (1970).

Na distribuição von Mises VM(µ, λ), o parâmetro µ representa a média circular de Y
enquanto que λ é denominado parâmetro de concentração, e é tal que o comprimento re-
sultante médio satisfaz ρ = A1 (λ). A dispersão circular de Y é dada por δ = [λA1 (λ)]−1
e as componentes de seus momentos trigonométricos são αp = Ap (λ) e βp = 0, para
p ≥ 1. A distribuição von Mises goza de propriedades bastante interessantes, dentre as
quais, destacamos

(i) Ela é unimodal e simétrica em torno da média direcional µ;

(ii) Quando λ → 0, a distribuição de von Mises converge para a distribuição circular

uniforme;

(iii) Quando λ → ∞, a distribuição converge para a distribuição degenerada no ponto

y = µ;

(iv) Segundo Stephens (1963), para valores não muito pequenos de λ, a distribui-
ção von Mises VM(µ, λ) é próxima de uma distribuição normal arqueada WN
(µ, 1/λ).

(v) A razão entre a densidade na moda e a densidade na antimoda é e2λ . Isto implica
que, quanto maior for o valor de λ maior será o agrupamento em torno da moda.
 31

2.6.1 Estimação dos parâmetros da Distribuição von Mises

Seja Y1 , . . . , Yn uma amostra aleatória independente de distribuição VM(µ, λ), com

valores observados y1 , . . . , yn . A função de log-verossimilhança fica dada por
n
X
l(µ, λ; y1 , . . . , yn ) = −n log 2π − n log I0 (λ) + λ cos(yi − µ).
i=1

Utilizando a equação (1.5) na expressão acima temos

l(µ, λ; y1 , . . . , yn ) = −n{log 2π − log I0 (λ) + λR cos(µ − µ)}. (2.21)

O estimador de máxima verossimilhança para µ é obtido quando a função de log-

verossimilhança atinge seu valor máximo. Como cos x tem valor máximo para x = 0,
então
cos(µ − µ) = 0 ⇒ µ − µ = 0 ⇒ µ
b=µ

ou seja, o estimador de máxima verossimilhança da média direcional µ é exatamente a

média circular amostral, definida pela equação (1.1). A propósito, Bingham e Mardia
(1975) mostram que a distribuição von Mises é a única distribuição circular onde a
direção média amostral é o estimador de máxima verossimilhança para a direção média
populacional.
Por sua vez, obtemos o estimador de máxima verossimilhança para λ da seguinte forma
∂l 1 ′
= n{R cos(µ − µ) − I (λ)}.
∂λ I0 (λ) 0
Mas, I0′ (λ) = I1 (λ), então
∂l I1 (λ)
= n{R cos(µ − µ
b) − }.
∂λ I0 (λ)
I1 (λ)
Como A1 (λ) = I0 (λ)
temos,

n{R cos(µ − µ
b) − A1 (λ)} = 0 ⇒ R cos(µ − µ
b) − A1 (λ) = 0 ⇒ A1 (λ) = R.

b para λ é tal que

Portanto, o estimador de máxima verossimilhança λ

b = R ou seja, λ
A1 (λ) b = A−1 (R)
1

b é viciado e uma aproximação para o vício foi obtida por Best e Fisher
O estimador λ
(1981). Eles concluem que, para amostras pequenas e valores pequenos de λ, a distribui-
b tem uma grande assimetria, com uma longa cauda à direita, conseqüentemente,
ção de λ
 32

as inferências baseadas em valores pequenos de λ podem ser suspeitas, a menos que

n seja suficientemente grande. Best e Fisher (1981) sugerem a seguinte aproximação
para A−1
1 .


 2r + r3 + (5r5 )/6, r < 0, 53


A−1
1 (r) ≃ −0, 4 + 1, 39r + 0, 43/(1 − r), 0, 53 ≤ r ≤ 0, 85



 (r3 − 4r2 + 3r)−1 , r > 0, 85.

Para pequenas amostras, n ≤ 15, Fisher(1995, p.88) sugere o seguinte estimador


 max{λ b − 2(nλ)
b −1 , 0}, se λ
b<2
λb∗ =
 (n−1)3 λb , b ≥ 2,
se λ
n3 +n

b é o estimador de máxima verossimilhança, λ

onde λ b = A−1 (R).
1

Distribuições Assintóticas

A partir da teoria padrão dos estimadores de máxima verossimilhança, segue-se que a

distribuição assintótica para grandes amostras de (b b é
µ, λ)

√
b − λ) ∼ N2 (0, I −1 ),
µ − µ, λ
n (b (2.22)

onde I é a matriz de informação de Fisher, dada por

  
∂2l ∂2l
∂µ2 ∂µ∂λ
I = E −   .
∂2l ∂2l
∂λ∂µ ∂λ2

Calculando as derivadas parciais e considerando as propriedades do valor esperado,

verificamos que a matriz de informação de Fisher é representada da seguinte forma:
 
λA(λ) 0
I= . (2.23)
2 A(λ)
0 1 − A (λ) − λ

Logo, para n grande, temos



 nV ar(b
µ) ≈ 1

 λA(λ)
b
nV ar(λ) ≈ 1
A(λ) (2.24)

 1−A2 (λ)− λ

 nCov(b b ≈
µ, λ) 0
 33

Portanto, µ b são aproximadamente distribuídos de acordo com duas normais inde-

beλ
pendentes, ou seja

√ ³ 1 ´
µ − µ) ∼ N 0,
n(b e
λA(λ)
à !
√ 1
b − λ) ∼ N 0,
n(λ .
2 A(λ)
1 − A (λ) − λ

2.7 Outras Distribuições Circulares

Apresentamos aqui as distribuições mais importantes para variáveis aleatórias
que assumem valores no círculo. Além dessas, existem outras distribuições, tais como
Cardióide, Lattice, entre outras. Mais detalhes podem ser encontrados em Batschelet
(1981), Fisher (1995) e Mardia (1972).
Capítulo 3

Modelos de Regressão para uma

Resposta Angular

Neste capítulo, consideraremos modelos de regressão onde a variável resposta é

uma variável angular sobre um conjunto de variáveis explicativas lineares. Apresenta-
remos modelos em que tanto a média direcional quanto a dispersão, de uma variável
com distribuição von Mises, estão relacionadas com as variáveis explicativas por meio
de funções de ligação apropriadas. Finalizaremos fazendo inferências para os modelos
de média, de dispersão e misto.

3.1 Introdução
Os modelos de regressão têm como objetivo fornecer um mecanismo capaz de
fazer previsões dos valores de uma variável, denominada resposta, utilizando informa-
ções sobre outras variáveis, que chamamos explicativas. Quando a variável resposta
é angular, a teoria dos modelos de regressão apropriada tem sido pouco considerada.
Em geral, a variável é tratada como sendo linear, embora problemas desse tipo não são
raros e estão presentes nas mais diversas áreas. Por exemplo, a partir de informações
sobre velocidades de ventos, fazer previsões com respeito a direções, ou ainda, a partir
de dados de distâncias percorridas por determinados animais, prever direções a serem
seguidas, entre outras aplicações. Seguindo esta abordagem, podemos destacar alguns
estudos como pioneiros.
 35

• Gould (1969), Laycock (1975) e Mardia (1972) consideraram várias propostas de

modelos em que a média da variável resposta Y como função de x é uma curva
em volta de infinitos espirais na superfície de um cilindro infinito.

• Jonhson e Wehrly (1978) propuseram uma classe de modelos onde a variável

resposta preenche apenas uma única espiral ao longo de x.

Devido a relevância destes modelos discutiremos os resultados dos mesmos e apresen-

taremos algumas extensões dos modelos de Jonhson e Wehrly.

3.2 Alguns Modelos para Regressão Angular

A modelagem de uma variável aleatória circular Y em termos de variáveis expli-
cativas pode ser feita sob as três abordagens seguintes:

• Modelo de Médias: Modela-se a média direcional de Y em termos do vetor das

covariáveis X = (X1 , . . . , Xp )T .

• Modelo de Dispersão: Modela-se a dispersão de Y em termos do vetor das cova-

riáveis X = (X1 , . . . , Xp )T .

• Modelo Misto: Modelam-se ambas, a média direcional e a dispersão de Y, em

termos do vetor das covariáveis X = (X1 , . . . , Xp )T .

Segundo Fisher e Lee (1992), quando a resposta é uma variável circular, os modelos de
regressão apresentam algumas características peculiares, dentre elas podemos destacar
a possibilidade de se modelar a dispersão circular tão bem quanto se modela a média,
uma vez que conjuntos de dados encontrados na prática exibem características de au-
mento de variabilidade de Y para pequenos valores de X, e outra razão é o fato de
existir distribuições, tais como a uniforme, em que a média direcional não está definida.
Vale salientar que modelar a dispersão não é uma tarefa simples, apresenta certas di-
ficuldades, porque para as distribuições circulares não existe uma medida natural de
escala. Por isso, é conveniente trabalhar com a família de distribuição von Mises,
que apresentam uma medida de dispersão intrínseca (λ) e, como foi dito no capítulo
anterior, esta distribuição compartilha de propriedades inferênciais semelhantes à dis-
tribuição normal nos dados lineares.
 36

Fazendo uma analogia à teoria de regressão normal linear, Gould (1969) propõe um
modelo para média direcional que apresenta a estrutura
p
X
µ = µ0 + βj xj , (3.1)
j=1

onde µ0 , β1 , . . . , βp são os parâmetros a serem estimados. Gould (1969) obteve um

sistema de equações de estimação referente aos estimadores de máxima verossimilhança
para os parâmetros desse modelo e desenvolveu um método iterativo para determinar
b que depende do ponto de iniciação do processo. Laycok (1975), estudando este
β,
modelo, verificou que o método da máxima verossimilhança era equivalente ao método
dos mínimos quadrados.
Posteriormente, Jonhson e Wehrly (1978) verificaram que a função de verossimilhança
para este modelo não atendia ao critério de unicidade, podendo gerar muitos picos
igualmente importantes. Sendo assim, sugeriram uma abordagem para uma variável
explicativa, por meio de um modelo específico para a distribuição conjunta de Y e
X, com uma função de distribuição marginal F (x) completamente especificada. A
distribuição condicional de Y dado X = x é VM (µ + 2πF (x), λ), um modelo que
permite estimação direta de µ e λ, o qual chamamos de modelo de médias de Jonhson
e Wehrly. Para o modelo de dispersão, com apenas uma variável explicativa, eles
sugeriram modelar a distribuição condicional de Y dado X = x como VM (µ, λx), que
também permite a estimação direta dos parâmetros.

3.3 Os Modelos de Regressão von Mises

Sejam Y1 , . . . , Yn observações angulares, que seguem distribuição von Mises com
médias direcionais µ1 , . . . , µn e parâmetros de concentração λ1 , . . . , λn , respectiva-
mente. Segundo Fisher e Lee (1992), esses parâmetros direcionais podem ser modelados
pelos modelos de regressão von Mises. Basta considerarmos que as n observações in-
dependentes estão relacionadas com x1 , . . . , xn , em que xi é um vetor de p variáveis
explicativas, com p ≤ n, através de funções de ligação adequadas.
 37

3.3.1 Modelo de Médias von Mises

No modelo de médias von Mises assumimos que todos os parâmetros de concentração

sejam iguais λ1 = · · · = λn = λ e que as médias direcionais µ1 , . . . , µn são relacionados
com p-covariáveis explicativas pela forma

µi = µ + g(β T xi ), (3.2)

onde g é uma função de ligação monótona, duas vezes continuamente diferenciável,

cujo objetivo é mapear a linha real para o círculo e β = (β1 , . . . , βp )T é um vetor de
parâmetros a serem estimados. Outra importante consideração que precisa ser feita
sobre a função de ligação g é que g(0) = 0, fazendo com que o parâmetro µ possa
ser visto como uma origem, uma vez que este conceito não está bem definido quando
a variável resposta assume valores no círculo. Fisher e Lee (1992) discutem possíveis
candidatas para a função g, no entanto uma função de ligação bastante utilizada é
g(t) = 2arctan(t).

Inferências para o Modelo de Média von Mises

Vamos considerar o modelo em que as variáveis aleatórias Y1 , . . . , Yn são independentes

e Yi ∼ VM (µi , λ), para i = 1, . . . , n, onde µi = µ + g(β T xi ) e g é conhecida. Neste
caso, a log-verossimilhança é
n
X
l ≡ l(µ, β, λ; y) = −n log I0 (λ) + λ cos[yi − µ − g(β T xi )] (3.3)
i=1

e o vetor paramétrico é dado por θ = (µ, β T , λ)T .

 38

Estimadores de Máxima Verossimilhança para o Modelo de Médias von

Mises

Inicialmente definimos:

ui = sen(yi − µ − g(β T xi ))

u = (u1 , . . . , un )T

X = (x1 , . . . , xn )T

G = diag(g ′ (β T x1 ), . . . , g ′ (β T xn ))
Xn
S = sen(yi − g(β T xi ))/n (3.4)
i=1
Xn
C = cos(yi − g(β T xi ))/n (3.5)
i=1

R = (S + C 2 )1/2 .
2
(3.6)

Calculando as derivadas parciais de (3.3), com relação aos parâmetros µ, β e λ, fazendo

as devidas substituições temos que
X n
∂l
= λ sen(yi − µ − g(β T xi ))
∂µ i=1
n
X
= λ sen[(yi − g(β T xi )) − µ]
i=1
n
X
= λ {[sen(yi − g(β T xi )) cos µ] − [senµ cos(yi − g(β T xi ))]}
i=1
= λ(nS cos µ − nCsenµ)

= λn(S cos µ − Csenµ).

Então,
∂l
= 0 ⇒ S cos µ
b − Csenb
µ = 0 ⇒ S cos µ
b = Csenb
µ.
∂µ
Desta forma, o estimador de máxima verossimilhança µ
b deve satisfazer

S 2 cos2 µ
b = C 2 sen2 µ
b.

Sabemos que cos2 θ + sen2 θ = 1, substituindo este resultado na expressão acima temos,

S 2 cos2 µ
b = C 2 (1 − cos2 µ
b) ⇒ S 2 cos2 µ
b + C 2 cos2 µ
b = C 2.

Assim,
R2 cos2 µ
b = C 2 ⇒ R| cos µ
b| = |C| ⇒ R cos µ
b = C.
 39

De fato, temos que se

π 3π
0≤µ
b< ou ≤µ b > 0 e C > 0 ⇒ | cos µ
b < 2π ⇒ cos µ b e |C| = C,
b| = cos µ
2 2
e se
π 3π
b < π ou π ≤ µ
≤µ b< b < 0 e C < 0 ⇒ | cos µ
⇒ cos µ b e |C| = −C.
b| = − cos µ
2 2
Portanto,
b| = |C| ⇒ R cos µ
R| cos µ b = C.

De forma análoga,

S 2 cos2 µ
b = C 2 sen2 µ
b ⇒ S 2 (1 − sen2 µ
b) = C 2 sen2 µ
b

⇒ R2 sen2 µ
b = S2

⇒ Rsenb
µ = S.

Podemos verificar que se

π π
0≤µ
b< ou ≤ µ
b < π ⇒ senb
µ > 0 e S > 0 ⇒ |senb
µ| = senb
µ e |S| = S,
2 2
e se
3π 3π
π≤µ
b< ou b < 2π ⇒ senb
≤µ µ < 0 e S < 0 ⇒ |senb µ e |S| = −S.
µ| = −senb
2 2
Portanto,
µ| = |S| ⇒ Rsenb
R|senb µ = S.

Logo, o estimador µ
b deve satisfazer o sistema de equações

 R cos µ b = C
(3.7)
 Rsenb µ = S.

Por outro lado,

Xn
∂l
= sen(yi − µ − g(β T xi ))g ′ (β T xi )xi
∂β i=1
n
X
= xi g ′ (β T xi )sen(yi − µ − g(β T xi )).
i=1

Então,
X n
∂l
=0⇒ xi g ′ (β T xi )sen(yi − µ − g(β T xi )) = 0.
∂β i=1
 40

Logo,
n
X
xi g ′ (β T xi )ui = 0 ⇒ xT Gu = 0.
i=1
Agora, com relação ao parâmetro de concentração λ, temos
X n
∂l
= −nA1 (λ) + cos[(yi − g(β T xi )) − µ]
∂λ i=1
n
X n
X
T
= −nA1 (λ) + cos µ cos(yi − g(β xi )) + senµ sen(yi − g(β T xi ))
i=1 i=1
= −nA1 (λ) + nC cos µ + nSsenµ

= n[−A1 (λ) + C cos µ + Ssenµ].

Portanto,
∂l b
= 0 ⇒ C cos µ
b + Ssenb
µ = A1 (λ).
∂λ
Utilizando a expressão (3.7) temos,

b ⇒ A1 (λ)
C 2 /R + S 2 /R = A1 (λ) b = R.

Portanto, os estimadores de máxima verossimilhança de µ, β e λ, devem satisfazer o

seguinte sistema de equações

X T Gu = 0 (3.8)

Rsenb
µ = S (3.9)

R cos µ
b = C (3.10)
b = R.
A1 (λ) (3.11)

Este sistema não possui uma solução explícita e deve ser resolvido por um processo
iterativo. Podemos obter uma solução através do algoritmo dos mínimos quadrados
iterativamente reponderados (IRLS), Green (1984). O vetor βb é atualizado pela equa-
ção
(m+1) (m)
X T G2 X(βb − βb ) = X T G2 y*, (3.12)
µ
bi
onde y* = (y1∗ , . . . , yn∗ )T e yi∗ = b ′ (β
[A1 (λ)g bT xi )] ·

Processo iterativo

b, βb e λ,
A obtenção das estimativas de máxima verossimilhança µ b é feita através do

seguinte processo iterativo:

 41

(i) Obter uma estimativa inicial βb0 ;

(ii) Calcular S, C e R usando as expressões (3.4), (3.5) e (3.6);

(iii) Determinar µ b por (3.9)-(3.11);

beλ

(iv) Atualizar βb pela equação (3.12);

(v) Repetir os passos (ii) a (iv) até obter a convergência.

3.3.2 Modelo de Dispersão von Mises

Como dito anteriormente, existem situações em que o interesse consiste em modelar o

parâmetro de concentração λ da distribuição von Mises. Para o modelo de dispersão,
assumimos que as médias direcionais são todas iguais (µ1 = · · · = µn = µ) e os
parâmetros de concentração λ1 , . . . , λn são modelados por

λi = h(γ T z i ), (3.13)

onde γ = (γ0 , γ1 , . . . , γq )T é um vetor de parâmetros a serem estimados, z i = (1, zi1 , . . . , ziq )T

é o vetor das covariáveis e h : R → [0, ∞) duas vezes continuamente diferenciável e
atua como uma função de ligação.

Inferências para o Modelo de Dispersão von Mises

Vamos considerar o modelo em que as variáveis aleatórias Y1 , . . . , Yn são independentes

e Yi ∼ VM (µ, λi ), para i = 1, . . . , n, onde λi = h(γ T z i ) e h é conhecida. Neste caso,
a log-verossimilhança é expressa como
n
X n
X
T
l ≡ l(µ, γ; y) = − log I0 (h(γ z i )) + h(γ T z i ) cos(yi − µ). (3.14)
i=1 i=1

e o vetor paramétrico é dado por θ = (µ, γ T )T .

 42

Estimadores de Máxima Verossimilhança para o Modelo de Dispersão von

Mises

Inicialmente definimos:
n
X
S = λi sen(yi ), (3.15)
i=1
n
X
C = λi cos(yi ) e (3.16)
i=1

R = (S + C 2 )1/2 .
2
(3.17)

Calculando as derivadas parciais de (3.14), com relação aos parâmetros µ, γ0 , γ1 , . . . , γq

e fazendo as devidas substituições temos que
Xn
∂l
= h(γ T z i )sen(yi − µ)
∂µ i=1
n
X
= λi sen(yi − µ)
i=1
Xn
= λi (senyi cos µ − cos yi senµ)
i=1
n
X n
X
= cos µ λi senyi − senµ λi cos yi
i=1 i=1
= S cos µ − Csenµ.

Então,
∂l
= 0 ⇒ S cos µ µ ⇒ S 2 cos2 µ
b = Csenb b = C 2 sen2 µ
b.
∂µ
Sabemos que cos2 θ + sen2 θ = 1. Substituindo este resultado na expressão acima temos,

S 2 cos2 µ
b = C 2 (1 − cos2 µ
b) ⇒ S 2 cos2 µ
b + C 2 cos2 µ
b = C 2.

Assim,
R2 cos2 µ
b = C 2 ⇒ R| cos µ
b| = |C| ⇒ R cos µ
b = C.

De forma análoga,

S 2 cos2 µ
b = C 2 sen2 µ
b ⇒ S 2 (1 − sen2 µ
b) = C 2 sen2 µ
b ⇒ R2 sen2 µ
b = S 2 ⇒ Rsenb
µ = S.

Portanto, o estimador µ
b deve satisfazer

 R cos µ
b = C
 Rsenbµ = S.
 43
d Ip (t)
Por outro lado, considerando que I (t)
dt 0
= I1 (t) e Ap (t) = I0 (t)
, conforme Apêndice D,
tem-se que
n
X I1 (λi ) ∂λi X n
∂l ∂λi
= − + cos(yi − µ)
∂γ0 I (λ ) ∂γ0 i=1
i=1 0 i
∂γ0
n
X ∂
= − A1 (λi )h′ (γ T z i ) (γ0 + γ1 z i1 + · · · + γq z iq ) +
i=1
∂γ0
n
X ∂
cos(yi − µ)h′ (γ T z i ) (γ0 + γ1 z i1 + · · · + γq z iq )
i=1
∂γ0
n
X
= h′ (γ T z i ){−A1 (λi ) + cos(yi − µ)}.
i=1

Então o estimador de γ0 deve satisfazer

n
X
γ T z i ){−A1 (λbi ) + cos(yi − µ
h′ (b b) = 0
i=1

Agora, para j = 1, . . . , q,
n
X n
X
∂l I1 (λi )
= − h′ (γ T zi )zij + cos(yi − µ)h′ (γ T z i )zij
∂γj i=1
I0 (λi ) i=1
Xn Xn
= − A1 (λi )h′ (γ T z i )zij + cos(yi − µ)h′ (γ T z i )zij
i=1 i=1
n
X
= h′ (γ T z i ){−A1 (λi ) + cos(yi − µ)}zij .
i=1

Desta forma, o estimador de máxima verossimilhança de γj deve satisfazer

n
X
γ T z i ){−A1 (λbi ) + cos(yi − µ
h′ (b b)}z ij = 0 para j = 1, . . . , q.
i=1

Portanto, o estimador de máxima verossimilhança θb = (b b T )T deve satisfazer o

µ, γ
sistema
n
X
γ T z i ){−A1 (λbi ) + cos(yi − µ
h′ (b b)} = 0, (3.18)
i=1
n
X
γ T z i ){−A1 (λbi ) + cos(yi − µ
h′ (b b)}zij = 0 para j = 1, . . . , q, (3.19)
i=1
R cos µ
b = C, (3.20)

Rsenb
µ = S. (3.21)

Como no modelo de médias, este sistema não admite solução explícita, logo a mesma
deve ser obtida através de algum método iterativo. A solução através do algoritmo do
 44

mínimos quadrados reponderados (IRLS), Green (1984), usa a equação de atualização

γ (m+1) − γ
Z T WZ(b b (m) ) = Z T Wy*, (3.22)

bi )
γ T z i )}2 A1 (λ
onde Z = (z 1 , . . . , z n )T , W é uma matriz diagonal com elementos wi = {h′ (b
cos(yi −b
µ)−A1 (λ bi )
e y* = (y1∗ , . . . , yn∗ )T , com yi∗ = bi ) .
γT zi )A1 (λ
h′ (b

Processo iterativo

A obtenção das estimativas de máxima verossimilhança µ

beγ
b , é feita através do seguinte
processo iterativo:

(i) Obter uma estimativa inicial γ

b0;

(ii) Calcular S, C e R usando (3.15) - (3.17) ;

(iii) Determinar µ
b por (3.20)-(3.21);

(iv) Atualizar γ
b pela equação (3.22);

(v) Repetir os passos (ii) a (iv) até obter a convergência.

3.3.3 Modelo Misto von Mises

Segundo Fisher (1995), podemos encontrar situações, envolvendo variáveis direcionais,

onde tanto a média direcional quanto o parâmetro de concentração dependem das
covariáveis. Nesses casos, nos referimos a mais uma classe de modelos de regressão,
denominada modelo misto von Mises, que combina o modelo de médias com o modelo
de dispersão. Para este modelo, a função de log-verossimilhança é dada por
n
X n
X
l ≡ l(µ, β, γ; y) = − T
log I0 (h(γ z i )) + h(γ T z i ) cos(yi − µ − g(β T xi )), (3.23)
i=1 i=1

onde o vetor paramétrico é θ = (µ, β T , γ T )T . Essa função é maximizada pela utilização

simultânea dos dois processos iterativos citados nas Secções 3.3.1 e 3.3.2, embora sejam
necessárias fazer as seguintes modificações:

(i) A matriz G2 da equação de atualização para β, (3.12), deve ser substituída por
bi A1 (λ
GΛG, onde Λ é uma matriz diagonal com elementos λ bi );

b por A1 (λ
(ii) Na definição do vetor y∗ em (3.12), deve-se substituir A1 (λ) bi );
 45

b , a i-ésima coordenada de y∗
(iii) Na equação (3.22), equação de atualização para γ
fica dada por
T
b − g(βb z i )) − A1 (λ
cos(yi − µ bi )
yi∗ = .
bi )
γ T z i )A1 (λ
h′ (b

Processo iterativo

b, βb e γ
Para obtermos as estimativas de máxima verossimilhança µ b , devemos seguir o
seguinte processo iterativo:

(i) Obter estimativas iniciais βb0 e γ

b0;

b por (3.12), considerando a devida modifica-

(ii) Obter uma estimativa atualizada β,
ção;

(iii) Obter uma estimativa atualizada γ

b , por (3.22), considerando a devida modifica-
ção;

(iv) Determinar R, S e C, por (3.15)-(3.17);

(v) Determinar µ
b por (3.20)-(3.21);

(vi) Repetir os passos (ii) a (v) até obter a convergência.

Capítulo 4

Resíduos para uma Resposta Angular

Neste capítulo estudaremos os resíduos para respostas angulares. Inicialmente,

faremos um estudo da componente da função desvio para a distribuição von Mises e
em seguida, apresentaremos duas formas padronizadas de resíduos para o modelo de
regressão von Mises e investigaremos suas respectivas distribuições de probabilidade.

4.1 Introdução
Podemos destacar a análise de diagnóstico como uma etapa de grande relevância
no ajuste de um modelo de regressão, uma vez que nela é possível verificar possíveis
afastamentos das suposições feitas para o modelo, bem como detectar a presença de
observações extremas, nos resultados do ajuste. A análise de diagnóstico foi iniciada
com a análise de resíduos, com a finalidade de detectar a presença de pontos extremos e
avaliar a adequação da distribuição proposta para a variável resposta, vide Paula (2004,
p.29). Quando a variável de interesse é direcional, como enfatizado anteriormente, já
podemos esperar que a metodologia utilizada para variáveis lineares não faça sentido.
Logo, faz-se necessário o desenvolvimento de uma teoria apropriada para dados dire-
cionais.
 47

4.2 A componente da função desvio

Após o ajuste de um modelo de regressão a um conjunto de dados, uma avaliação
da qualidade desse ajuste sempre é realizada por meio de uma estatística especial. Nos
modelos lineares generalizados, uma medida que avalia a qualidade do ajuste bastante
P
b = ni=1 d2i , onde
utilizada é a função desvio, definida por D(y, µ)
√
di ≡ di (yi ; µbi ) = ± 2[li (yi ; µei ) − li (yi ; µbi )]1/2 , (4.1)

em que li (yi ; . . . ) é a contribuição de yi para a log-verossimilhança, µei é o estimador de

máxima verossimilhança baseado apenas em yi e µbi é o estimador de máxima verossi-
milhança baseado na amostra completa, segundo (Souza, 1999) e (McCullagh e Nelder,
1989). O sinal de di representa a direção de discrepância. Em particular, se Y1 , . . . , Yn
representam variáveis aleatórias independentes com Yi ∼VM(µi , λ), então a função de
log-verossimilhança é dada por
n
X
l(µi ; λ, yi ) = −n log I0 (λ) + λ cos(yi − µi ),
i=1

de onde podemos escrever

li (yi ; µi , λ) = − log I0 (λ) + λ cos(yi − µi ).

Logo,

li (yi ; µei , λ) = − log I0 (λ) + λ cos(yi − µei ) e

li (yi ; µbi , λ) = − log I0 (λ) + λ cos(yi − µbi ).

Portanto, substituindo as expressões acima em (4.1), temos

√
di = ± 2[− log I0 (λ) + λ cos(yi − µei ) + log I0 (λ) − λ cos(yi − µbi )]1/2

Como µei é o estimador de máxima verossimilhança baseado apenas em yi , temos que

cos(yi − µei ) = 1. Logo,
√
di = ± 2[λ − λ cos(yi − µbi )]1/2
√
= ± 2[λ(1 − cos(yi − µbi ))]1/2 .

Então, a componente da função desvio é

√
di ≡ di (yi ; µbi , λ) = ± 2λ[1 − cos(yi − µbi )]1/2 . (4.2)
 48

Utilizando algumas propriedades das funções trigonométricas podemos escrever (4.2)

de uma outra maneira mais apropriada,
· ¸
√ 1
di = di (yi ; µbi , λ) = ±2 λsen (yi − µbi ) . (4.3)
2

De fato, observe inicialmente que

1 1
cos(yi − µbi ) = cos[ (yi − µbi ) + (yi − µbi )],
2 2

de onde se obtém

1 1
cos(yi − µbi ) = cos2 [ (yi − µbi )] − sen2 [ (yi − µbi )].
2 2

Assim,

1 1 1
1 − cos(yi − µbi ) = 1 − {cos2 [ (yi − µbi )] − sen2 [ (yi − µbi )]} = 2sen2 (yi − µbi ).
2 2 2

Portanto, substituindo esta expressão em (4.2) temos,

· ¸1/2
√ 21
di = di (yi ; µbi , λ) = ± 2λ 2sen (yi − µbi )
2
· ¸
√ 1
= ±2 λsen (yi − µbi ) .
2

4.3 Resíduos
De uma maneira geral, o resíduo para a i-ésima observação pode ser definido
como uma função ri = r(yi , µ
bi ), que tem como objetivo medir a discrepância entre
o valor observado e o valor ajustado para a i-ésima observação (Cox e Snell, 1968).
No modelo de regressão normal linear, em que Yi ∼ N (xTi β, σ 2 ), o resíduo ordinário
b Todavia há outras formas de definir resíduos. Nesse
é definido por ri = yi − xTi β.
modelo, a estimativa do vetor de médias é µ b onde βb = (XT X)−1 XT y. Logo
b = Xβ,
b = Hy, com H = X(XT X)−1 XT . A matriz H é a matriz de
podemos escrever µ
projeção ortogonal dos vetores de Rn no subespaço gerado pelas colunas da matriz
modelo X, conhecida como matriz hat. Portanto, considerando o vetor de resíduos
ordinários definido por r = (r1 , . . . , rn )T , segue que

r=y−µ
b = y − Hy = (I − H)y.
 49

Para este modelo, segue que E(r) = 0 e V ar(r) = σ 2 (I−H), ou seja, ri tem distribuição
′
normal com média zero e variância V ar(ri ) = σ 2 (1 − hii ). Como os ri s possuem
variâncias distintas , é conveniente que eles sejam expressos de uma forma padronizada
para que seja possível fazer comparações entre os mesmos. Uma forma natural é obter
o resíduo studentizado, dividindo ri pelo respectivo desvio padrão, logo
ri
ti = , i = 1, . . . , n, (4.4)
s(1 − hii )1/2
Pn ri2
onde, s2 = i=1 (n−p) e p é o número das covariáveis. No entanto, ti não segue
uma distribuição t-student, já que ri não é independente de s2 . Para contornar este
problema, substituímos s2 por s2(i) , que é o erro quadrático médio correspondente ao
modelo sem a i-ésima observação. Assim o novo resíduo studentizado é dado por
ri
t∗i = , (4.5)
s(i) (1 − hii )1/2
que segue uma distribuição tn−p−1 central.
Para os modelos lineares generalizados, a definição de um resíduo studentizado pode
ser feita de forma análoga à regressão normal linear. Todavia algumas propriedades
não continuam valendo, por isso faz-se necessária a definição de outros resíduos cujas
propriedades sejam conhecidas ou se aproximam das propriedades de t∗i . Nos MLGs os
resíduos mais utilizados são definidos a partir das componentes da função desvio, cuja
versão padronizada (vide McCullagh, 1987; Davison e Gigli, 1989) é dada por
d(yi ; µi )
tDi = q . (4.6)
1−b hii
Williams(1984) verificou, por meio de simulações, que a distribuição de tDi tende a
estar mais próxima da normalidade do que as distribuições de outros resíduos citados
na literatura.
Podemos encontrar muitas informações sobre resíduos em Cox e Snell (1968), McCul-
lagh e Nelder (1989), entre outros.

4.3.1 Resíduo d∗i

McCullagh (1987) mostra que, para os modelos lineares generalizados, a distribuição

de probabilidade de
d(yi ; µi ) + ρ3i /6
p
1 + (14ρ23i − 9ρ4i )/36
 50

é aproximadamente N (0, 1), em que ρ3i e ρ4i são os coeficientes de assimetria e cur-
∂L(ηi )
tose de ∂ηi
, respectivamente, e d∗ (yi ; µi ) é a i-ésima componente da função des-
vio D∗ (y; µ).
b Usando os resultados de Cox e Snell (1986) é possível mostrar que
E(d∗ (yi ; µi )) = 0 e V ar(d∗ (yi ; µi )) = 1 − hii , em que os termos negligenciados são
o(n−1 ). Para este caso, a matriz hat é dada por H = W1/2 X(XT WX)−1 XT W1/2 onde
W é uma matriz diagonal de pesos wi = E(−∂ 2 li /∂ηi2 ), calculado na convergência do
processo iterativo do MLG. Baseado nesta característica da matriz hat, Souza e Paula
(2002) consideraram uma correção para a componente da função desvio, equação (4.3),
obtendo o resíduo padronizado d∗i dado por
∗
√ sen 1 (yi − µbi )
2
di = ±2 λ , (4.7)
(1 − b ∗
hii )1/2
onde a matriz hat é H∗ = ĠX(X G˙ 2 X)−1 X Ġ com Ġ = diag(g ′ (βxi )) e b
′ ′ ′
h∗ii é o ele-
mento da diagonal principal de H∗ , avaliada na estimativa de máxima verossimilhança.
Utilizando procedimentos análogos a Williams (1984), para os modelos lineares genera-
lizados, Souza e Paula (2002) investigaram a distribuição de probabilidade do resíduo
(4.7), por meio de estudo de simulação, e verificaram que a distribuição deste resíduo
apresenta concordância com a distribuição normal padrão.

4.3.2 Resíduo ri

Souza e Paula (2002) fizeram outra transformação na componente da função desvio,

baseando-se na abordagem empregada por Davison e Gigli (1989). Tal abordagem
tem como objetivo expressar a componente do desvio como função de uma variável
aleatória com distribuição conhecida, e em seguida fazer uma expansão em série de
Taylor até segunda ordem. Considerando a transformação y − µ = √z , na função
λ

densidade de probabilidade da distribuição von Mises, dada por (2.19), e usando o fato
de que cos t = 1 − 21 t2 + o(t2 ), quando t → 0, verificamos que, como a função densidade
de probabilidade da von Mises é da forma
1
fY (y; µ, λ) = eλ cos(y−µ) ,
2πI0 (λ)
tem-se
z
exp{λ cos(y − µ)} = exp{λ cos( √ )}
λ
2
z z2
= exp{λ − + o( )}.
2 λ
 51

Mas, pela continuidade da função exponencial exp{o(z 2 )} = 1 + o(z 2 ). Então a função

densidade de probabilidade VM(µ, λ) fica da forma

z 1 z2
f ( √ ; µ, λ) = exp{λ − + 1 + o(z 2 )}
λ 2πI0 (λ) 2
−z 2 /2 λ
e e
= √ √ {1 + o(z 2 )} (4.8)
2π 2πI0 (λ)

Assim verificamos que a função densidade de probabilidade da VM(µ, λ) é aproxima-

damente proporcional à função densidade de probabilidade da normal padrão N (0, 1)
para valores próximos de µ. É válido lembrar que, na propriedade (iv) da seção (2.6),
Stephens (1963) mostrou que existe relação entre a distribuição von Mises e a distri-
buição normal.
A abordagem de Davison e Gigli (1989) requer que seja conhecida a mediana da distri-
buição von Mises VM(µ, λ), mas, devido à simetria desta distribuição, vamos conside-
rar a seguinte transformação para podermos definir esta medida. Seja Y ∼ VM(µ, λ)
no intervalo [−π, π), considere W0 = (Y − µ)( mod 2π) e a transformação


 W0 + 2π, se −2π < W0 < −π


W = W0 , se −π ≤ W0 < π (4.9)



 W − 2π, se
0 π ≤ W < 2π.0

Então, W ∼ VM(0, λ) no intervalo [−π, π), logo a mediana de W é wm = 0, similar

à média direcional. De acordo com a abordagem usada por Davison & Gigli (1989),
considere que W tem função de distribuição acumulada FW (w; ν, λ) onde ν representa
a média direcional.
Pelo teorema da Transformação Integral, podemos escrever o desvio residual da seguinte
forma:
−1
d(W ; ν, λ) = d(FW (Φ(Z)); ν, λ) = T (Z; ν, λ). (4.10)

onde, Z ∼ N (0, 1) e T é uma função desconhecida. A equação (4.10) está demonstrada

no Apêndice C. De (4.10) podemos concluir que o desvio residual d pode ser obtido
pelo estudo da transformação T . Mas, como na VM(µ, λ) a função de distribuição
acumulada não possui uma forma explícita, então não podemos determinar uma forma
analítica para a função T . No entanto, podemos recorrer aos dois primeiros termos da
expansão em série de Taylor de T (z; ν, λ) em torno de z = 0. Como wn é a mediana de
 52

W ∼ VM(ν, λ) e f como sendo sua respectiva função de densidade de probabilidade,

pela demonstração de Davison & Gigli temos,
1
T (z; ν, λ) = T (0; ν, λ) + T ′ (z; ν, λ)z + T ′′ (z; ν, λ)z 2 + R2 (z, 0), (4.11)
2
onde,

T (0; ν, λ) = d(wn ; ν, λ) (4.12)

d(wn ; ν, λ)
T ′ (0; ν, λ) = (4.13)
(2π)1/2 f (wn )
1 ³ f ′ (wn )d′ (wn ; ν, λ) ´
T ′′ (0; ν, λ) = d ′′
(wn ; ν, λ) − (4.14)
2πf 2 (wn ) f (wn )
Demonstração das expressões (4.12) - (4.14):
(4.12) De (4.10) temos que,

−1
T (0; ν, λ) = d(FW (Φ(0)); ν, λ) = d(wn ; ν, λ)

(4.13) Calculando a derivada da função T em relação à variável z, temos,

−1
∂T ∂[d(FW (Φ(z)))] −1 1 1 2
= = d′ [FW (Φ(z))] 1/2
e−z /2 .
∂z ∂z f (wn ) (2π)
Quando z → 0,
d′ (wn ; ν, λ)
T ′ (0; ν, λ) = .
(2π)1/2 f (wn )
(4.14) Encontramos a expressão (4.14), derivando a expressão (4.13) e analisando-a
quando z → 0.
h d′ (w ; ν, λ)e−z2 /2 i′
n
T ′′ (z; ν, λ) =
(2π)1/2 f (wn )
2 /2 2 /2
[d′ (wn ; ν, λ)e−z ]′ (2π)1/2 f (wn ) − d′ (wn ; ν, λ)e−z [(2π)1/2 f (wn )]′
=
2πf 2 (wn )
Agora
−1
2 /2 d −1 d[FW (Φ(z))] d[Φ(z)] −z2 /2
[d′ (wn ; ν, λ)e−z ]′ = −1 d′
([FW (Φ(z))]) e
d[FW (Φ(z))] d[Φ(z)] dz
−1 2 /2
= +d′ (d[FW (Φ(z))])[−ze−z ]
−z 2 /2
d′′ (wn ; ν, λ)e 2 /2
= − zd′ (wn ; ν, λ)e−z
(2π)1/2 f (wn )

£ ¤′ −1
−1 d[FW (Φ(z))] d[Φ(z)]
(2π)1/2 f (wn ) = f ′ (FW (Φ(z)))
d[Φ(z)] dz
2 /2
f ′ (wn )e−z
= .
f (wn )
 53

Substituindo estas expressões verificamos que,

h 2 i
′ )e−z /2
e−z /2 d′′ (wn ; ν, λ) − zd′ (wn ; ν, λ)f (wn )(2π)1/2 − d′ (wn ; ν, λ) f (wfn(w
2

n)
T ′′ (z; ν, λ) = 2
.
2πf (wn )
Portanto, quando z → 0 temos,
µ ¶
′′ 1 ′′ f ′ (wn )d′ (wn ; ν, λ)
T (0; ν, λ) = 2
d (wn ; ν, λ) − .
2πf (wn ) f (wn )
Para W ∼ VM(ν, λ), vimos que a componente da função desvio pode ser determinado
pela equação (4.3) ficando, neste caso, da forma
· ¸
√ 1
d(w; ν, λ) = 2 λsen (w − ν) . (4.15)
2
Derivando a expressão (4.15) duas vezes temos que,
· ¸
′
√ 1
d (w; ν, λ) = λ cos (w − ν)
2
· ¸
′′ −1 √ 1
d (w; ν, λ) = λsen (w − ν) .
2 2
Como W ∼ VM(ν, λ), sua função de densidade de probabilidade é dada por
1
f (w) = eλ cos(w)
2πI0 (λ)
e a primeira derivada por
−λsen(w) λ cos(w)
f ′ (w) = e .
2πI0 (λ)
Substituindo estas expressões em (4.11) e usando (4.8) segue-se que

d(w; 0, λ) = T (z; 0, λ) = (2λπ)1/2 I0 (λ)e−λ z + o(w2 ), quando w → 0 (4.16)

Uma vez que, de (4.11), temos,

1
T (z; 0, λ) = T (0; 0, λ) + T ′ (0; 0, λ)z + T ′′ (0; 0, λ)z 2 + R2 (z, 0),
2
onde,
µ ¶
√ 1
T (0; 0, λ) = d(wn ; 0, λ) = 2 λsen w
2
√ ¡ ¢
d(wn ; 0, λ) 2 λsen 21 w
T ′ (0; 0, λ) = =
(2π)1/2 f (wn ) (2π)1/2 f (w)
1 ³ f ′ (wn )d′ (wn ; 0, λ) ´
T ′′ (0; 0, λ) = ′′
d (wn ; 0, λ) −
2πf 2 (wn ) f (wn )
µ −λsen(w)eλ cos(w)
√ ¶
1 −1 2πI0 (λ)
λ cos 12 (w)
= e2(λ cos(w))
sen(w) − eλ cos(w)
.
2
2
2πI0 (λ) 2πI0 (λ)
 54

Portanto, quando w → 0 tem-se

d(w; 0, λ) = T (z; 0, λ) = (2λπ)1/2 I0 (λ)e−λ z + o(w2 ).

Com este resultado podemos aproximar a componente da função desvio de VM(0, λ)

por uma função linear da distribuição normal padrão N (0, 1). Mas, temos também
que, para µ 6= 0, segue de (4.9), que
√
d(wi ; νi , λ) = 2λ(1 − cos(yi − µi ))1/2 = d(yi ; νi , λ), (4.17)

ou seja, a componente da função desvio da VM(µ, λ) segue a mesma distribuição do

componente da função desvio para VM(0, λ). Pelos termos acima, Souza e Paula (2002)
apresentaram o seguinte resíduo padronizado
µ ¶1/2
2 sen( 12 (yi − µbi ))
ri = ± , (4.18)
π b −λb
I0 (λ)e

que foi tratado como o resíduo d∗i e comprovado que sua distribuição de probabilidade
é equivalente à distribuição normal padrão N (0, 1).
Capítulo 5

Aplicação

Neste capítulo apresentaremos a aplicação que realizamos, considerando o con-

junto de dados descrito em Fisher (1995, p.251), representando 19 observações de
direção de ventos e concentração de ozônio.

5.1 Introdução
Desde os anos 90, o mundo tem se voltado para as questões ambientais, dentre
elas, a poluição, os problemas na camada de ozônio, o efeito estufa entre outros. Re-
centemente, o aumento da temperatura global da Terra, tem preocupado, não apenas
os ambientalistas, como também a comunidade científica em todo o mundo. O efeito
estufa, bem como o buraco na camada de ozônio, tem aumentado a temperatura no
planeta e este aumento tem causado uma série de fenômenos catastróficos marcados
pelas mudanças de tempo em vários lugares, intensificação de ventos, chuvas e tem-
pestades. Por isso, nosso objetivo, neste trabalho, é investigar a associação entre a
direção dos ventos e a concentração de ozônio, de forma a verificar se as direções dos
ventos são influenciadas pela concentração de ozônio da atmosfera. É bem verdade que
vale a recíproca, a concentração de ozônio é influenciada pela direção dos ventos, tal
abordagem pode ser vista em Fisher (1995).
 56

5.2 Descrição e Análise dos Dados

Nesta aplicação, consideramos o conjunto de dados descrito no Apêndice B.18 de
Fisher (1995, p. 251), representando 19 observações de direção de ventos, expressa em
graus e concentração de ozônio, expressa em g/m3 . Essas observações foram coletadas
sempre às 6:00 h da manhã, num intervalo de 4 dias, no período de 18 de abril a 29
de junho de 1975, por uma estação de tempo em Milwaukee. Os valores desses dados
podem ser vistos pela Tabela 5.1 a seguir.

Tabela 5.1: Medidas de Direção de Ventos e Concentração de Ozônio

Direção Concentração Direção Concentração
327 28.0 8 112.0
91 85.2 204 20.0
88 80.5 86 72.5
305 4.7 333 16.0
344 45.9 18 45.9
270 12.7 57 32.6
67 72.5 6 56.6
21 56.6 11 112.0
281 31.5 27 91.8
84 55.2
Fonte: Fisher (1995)

Como estamos considerando a variável resposta como sendo a direção dos ven-
tos (circular) e a concentração de ozônio como variável explicativa (linear), então a
associação entre as variáveis é do tipo circular-linear com o modelo de regressão von
Mises.
Para termos melhor visão das componentes sistemáticas nas observações circula-
res, na Figura 5.1 apresentamos um gráfico sugerido por Fisher (1995), onde o eixo das
abscissas corresponde à concentração de ozônio e no eixo das ordenadas consideramos
yi como sendo a i-ésima observação de direção de ventos e yi + 2π.
 57

Figura 5.1: Gráfico da Concentração de Ozônio X Direção dos Ventos

Pelo gráfico, percebemos a possível existência de uma tendência na direção média.

Desta forma vamos ajustar o seguinte modelo de médias:

(i) yi segue uma distribuição von Mises VM(µi , λ) para i = 1, . . . , 19;

(ii) µi = µ + 2 arctan(βxi ), onde xi representa a concentração de ozônio de i-ésima

observação.

Aplicando o procedimento iterativo descrito em (3.12), através do aplicativo R (dispo-

nível em http://www.r-project.org), obtivemos as seguintes estimativas de máxima
verossimilhança, com respectivos desvios padrão assintóticos, que estão apresentados
na Tabela 5.2, indicando que há indícios de que ajuste é aceitável.
 58

Tabela 5.2: Estimativas dos parâmetros do modelo de médias von Mises

Parâmetro Estimativa D. Padrão Estim/ D. Padrão Nível descritivo
µ -1.492 0.1903 -7.8402 −
β 0.025817 0.005707 4.5237 <0.0001
λ 2.135 0.5988 3.5655 −

Análise dos Resíduos

Nas Figuras 5.2 e 5.3, apresentamos os gráficos normal de probabilidade para os re-
síduos padronizados ri e d∗i , respectivamente. Em cada gráfico apresentamos bandas
de confiança construídas através da técnica de envelopes, cujos limites foram obtidos
a partir de 19 simulações, segundo Atkison (1985). Desta forma, a probabilidade do
resíduo observado ri ou (d∗i ) exceder o limite superior do envelope é aproximadamente
1/20. Ao analisarmos esses gráficos não verificamos afastamentos das suposições feitas
inicialmente para o ajuste, indicando assim que o modelo de médias foi ajustado de
forma satisfatória.

Figura 5.2: Gráfico dos envelopes para o resíduo-r

 59

Figura 5.3: Gráfico dos envelopes para o resíduo-d∗

Conclusões

Reforçando o que já foi dito neste trabalho, a análise de dados direcionais não
precisa, nem deve, ser feita utilizando a teoria desenvolvida para dados lineares, pois
dispomos de ferramentas eficientes e mais apropriadas para análises de dados direcio-
nais, quer sejam análises descritivas exploratórias ou mesmo quando se deseja ajustar
modelos de regressão, dentre outras. Os resíduos considerados neste trabalho, são ca-
pazes de indicar a presença de informações discrepantes, como também apresentam
propriedades que nos permitem fazer inferências baseadas na aproximação pela distri-
buição normal padrão, tão importante e conhecida da estatística linear.
Apêndice A

Demonstração das Equações (1.5) e

(1.6)

Neste apêndice apresentaremos a demonstração das equações (1.5) e (1.6) citadas

no Capítulo 1.

A.1 Demonstração
Aplicando a propriedade do cosseno de uma diferença na expressão (1.5), obser-
vamos que
n n
1X 1X
cos(yi − µ) = (cos yi cos µ + senyi senµ)
n i=1 n i=1
" n n
#
1 X X
= cos µ cos yi + senµ senyi
n i=1 i=1

De (1.4) podemos escrever

n
1X
cos(yi − µ) = C cos µ + Ssenµ
n i=1

Substituindo (1.2) e (1.3) na expressão acima temos

n
1X
cos(yi − µ) = R cos2 µ + Rsen2 µ
n i=1
= R(cos2 µ + sen2 µ)

= R.
 62

Portanto,
n
1X
cos(yi − µ) = R.
n i=1
Por outro lado, de maneira análoga temos
n n
1X 1X
sen(yi − µ) = (senyi cos µ − senµ cos yi )
n i=1 n i=1
" n n
#
1 X X
= cos µ senyi − senµ cos yi .
n i=1 i=1

Fazendo as devidas substituições temos que,

n
1X
sen(yi − µ) = S cos µ − Csenµ
n i=1
SC CS
= −
R R
= 0.

Logo,
n
1X
sen(yi − µ) = 0.
n i=1
Apêndice B

Alguns Resultados sobre Dados

Direcionais

Neste apêndice demonstraremos as propriedades da função característica φY (p)

e apresentaremos uma demonstração para desigualdade tipo Tchebychev, que foram
citadas no Capítulo 1.

B.1 Demonstração das Propriedades da Função Ca-

racterística φY (p)
Como visto anteriormente, a função característica φY (p) goza das seguintes pro-
priedades:

(i) φY (0) = 1;

(ii) φY (p) = φ−Y (p), onde φY (p) é o conjugado complexo de φY (p);

(iii) |φY (p)| ≤ 1;

(iv) É suficiente considerarmos a função φY (p) apenas para valores inteiros e positivos
de p.

Demonstração: A propriedade (i) é consequência da definição. Ou seja,

φY (0) = E(ei0Y ) = E(1) = 1.

 64

(ii) Seja φY (p) = αp − iβp . De (1.21) temos que

φY (−p) = E(e−ipY ) = E[cos(−pY )] + iE[sen(−pY )]

= E[cos(pY )] − iE[sen(pY )]

= αp − iβp

= φY (p).

(iii) A propriedade (iii) é consequência das propriedades dos momentos trigonométricos

e do módulo. Ou seja,

|φY (p)| = |E(eipY )| ≤ E|eipY | = 1.

(iv) Usando a teoria de séries de Fourier, verifica-se que é suficiente considerar a função
φY (p) apenas para valores inteiros positivos de p. Sejam os números complexos φY (p)
são os coeficientes de Fourier de F (Feller, 1996, p.595 ou Zygmund, 1959, p.11).
Quando os φY (p) são dados por (1.21), temos que a expansão da série de Fourier da F
é
∞
1 X
dF (y) ∼ φY (p)e−ipY .
2π p=−∞

A variável aleatória circular Y tem função de densidade dada por

∞
1 X
f (y) = φY (p)e−ipY , (B.1)
2π p=−∞
P∞ 2
uma vez que, p=−∞ (αp + βp2 ) é convergente em L2 .
Pelo Teorema da Inversão, considerado em variáveis aleatórias contínuas na reta,podemos
escrever (1.22) na forma
−1 ∞
1 X X
f (y) = { (αp cos py + βp senpy) + 1 + (αp cos py + βp senpy)}.
2π p=−∞ p=1

Pelas propriedades das funções trigonométricas e dos momentos circulares, temos que
X ∞
1
f (y) = {1 + (αp cos py + βp senpy)}. (B.2)
2π p=1

Portanto, é suficiente considerarmos a função φY (p) apenas para valores inteiros e

positivos de p.
 65

B.2 Demonstração para a Desigualdade tipo Tcheby-

chev
£ ¤
Temos que ρ = E[cos(Y − µ)] = E{cos ( Y −µ
2
) + ( Y −µ
2
) }. Aplicando a propri-
edade do cosseno de uma soma, bem como algumas propriedades do valor esperado
temos
· ¸ · ¸
2 Y −µ 2 Y −µ
ρ = E cos ( ) − E sen ( )
2 2
· ¸ · ¸
2 Y −µ 2 Y −µ
= E 1 − sen ( ) − E sen ( )
2 2
· ¸
2 Y −µ
= 1 − 2E sen ( ) .
2

Logo, · ¸
Y −µ
2 1−ρ ν
E sen ( ) = = .
2 2 2
Assim, · ¸ Z π
Y −µ2 y−µ
ν = 2E sen ( ) =2 sen2 ( )dF (y)·
2 −π 2
Agora, dado 0 < ε ≤ 1, considere Aε = {y : |sen( y−µ
2
)| ≥ ε}. Assim,
·Z Z ¸ Z
2 y−µ 2 y−µ y−µ
ν=2 sen ( )dF (y) + sen ( )dF (y) ≥ 2 sen2 ( )dF (y).
Aε 2 Acε 2 Aε 2

Então, Z Z
2 2
ν≥2 ε dF (y) = 2ε dF (y) = 2ε2 P (Y ∈ Aε ).
Aε Aε

Portanto, µ ¶
Y −µ ν
P |sen( )| ≥ ε ≤ 2 .
2 2ε
Apêndice C

Demonstração da Equação (4.10)

Neste apêndice enunciaremos o Teorema da Transformação Integral e uma de-

monstração para a equação (4.10) citada no Capítulo 4.

C.1 Enunciado do Teorema da Transformação Inte-

gral
Teorema: Se X é uma variável aleatória cuja função de distribuição F é uma
função contínua na reta, então a distribuição de Y = F (X) é U [0, 1].

C.2 Demonstração da Equação (4.10)

−1
De fato, as igualdades da equação (4.10) são observadas. Seja W = FW (FW (W ))
pelo teorema da Transformação Integral, podemos considerar que,

FW (W ) = U1 e Φ(Z) = U2

onde Ui ∼ U [0, 1] e U1 = U2 quase sempre.

Assim,
−1 −1
FW (U1 ) = W = FW (Φ(Z))

como U1 = U2 ⇒ FW (W ) = Φ(Z) quase sempre. Logo,

−1
W = FW (Φ(Z)).
 67

Portanto,

−1
∀ W tal que FW (W ; ν, λ) > 0 ∃! z ∈ R tal que W = FW (Φ(Z))

Logo,
−1
d(W ; ν, λ) = d(FW (Φ(Z)); ν, λ) = T (Z; ν, λ).
Apêndice D

Alguns resultados sobre as funções de

Bessel

Neste apêndice apresentaremos alguns resultados relacionados com a função de

Ip (t)
Bessel Ip (t), e também com a função Ap (t) = I0 (t)
. Como visto anteriormente, essas
funções são de grande importância para a distribuição von Mises. Destacaremos aqui,
alguns resultados que foram extraídos de Fisher (1995) e Abramowits & Stegun (1970).

D.1 Derivadas
A função de Bessel modificada do primeiro tipo e ordem p, é definida pela série
∞
X 1 1
Ip (λ) = ( λ)2r+p , p = 0, 1, . . .
r=0
(p + r)!r! 2

As derivadas de maior interesse no nosso caso, de funções envolvendo I0 (t) e I1 (t) são:

d
(i) I (t)
dt 0
= I1 (t);

d I1 (t)
(ii) I (t)
dt 1
= I0 (t) − t
= 12 {I0 (t) + I2 (t)};

d
(iii) dt
{tI1 (t)} = tI0 (t);

d A1 (t)
(iv) A (t)
dt 1
=1− t
− A21 (t);

d2 A1 (t)
© ª
(v) A (t)
dt2 1
= t2
− 2A1 (t) + 1t dtd A1 (t).
 69

D.2 Fórmulas de recorrência

Para p = 2, 3, . . . e t > 0 dispomos das seguintes fórmulas de recorrência:

2(p−1)
(i) Ip (t) = Ip−2 (t) − t
Ip−1 (t)

2(p−1)
(ii) Ap (t) = Ap−2 (t) − t
Ap−1 (t).
Apêndice E

Programas Computacionais

Neste apêndice apresentaremos os programas computacionais, escritos para o apli-

cativo R, que foram utilizados para realizar o ajuste do modelo de média von Mises,
para construção dos gráficos dos envelopes e dos demais gráficos deste trabalho.

E.1 Comandos
###Este programa contém os comandos para ajustar o modelo de médias,
usando a função "lm.circular" # do pacote CIRCULAR ###

#LEITURA DOS DADOS#

ozonio.dat<-read.table("ozonio.dat",header=TRUE)
ozonio<-ozonio.dat$concentracao direcao<-ozonio.dat$direcao

direcao.real <- circular(ozonio.dat$direcao, type = c("angles"),

units = c("degrees"), modulo = c("2pi"), zero = 90 , rotation =
c("clock"))
ajuste <- lm.circular(y=direcao.real, x=ozonio, init=0,type=c("c-l"))
print(ajuste)
#--------------------------------------------------------------------#

### # Este programa faz o gráfico de ’ENVELOPES’ # para componente

do desvio padronizado no # modelo de MEDIAS von Mises.

##LEITURA DOS DADOS ##

 71

ozonio.dat <-read.table("ozonio.dat", header=TRUE)

attach(ozonio.dat)

namos <- length(ozonio.dat$direcao) direcao <-

matrix(c(ozonio.dat$direcao),nrow=namos,ncol=1)

##TRASFORMANDO A DIRECAO DE ’graus’ PARA ’radianos’

require(circular)
require(CircStats)
direcao <- rad(direcao)# <--- angulos em [0 , 2*pi)

#TRANSFORMANDO A DIRECAO PARA O INTERVALO [-pi , +pi] ## for(k in

1:namos) {if (direcao[k] > pi) direcao[k] <- direcao[k]-2*pi}

#DEFININDDO AS FUNCOES DE LIGACAO E SUAS DERIVADAS

g <-function(a) { 2*atan(a) }
dg <-function(a) { 2/(1+a^2) }
d2g<-function(a) { -(4*a)/((1+a^2)^2) }
h <-function(a) { exp(a) }
dh <-function(a) { exp(a) }
d2h<-function(a) { exp(a) }

#source("bessel1.fun") #source("bestfish.fun") source("modulo.fun")

source("arctg.fun") source("lambda.fun") source("davgigli.fun")

#AJUSTE: estimativas dos parâmetros para os dados de OZÔNIO

#Modelo de MEDIAS
mu <- -1.492 # d.p. =0,1903) # mu ~ ?? graus
beta <- matrix(c(0.025817),nrow=1,ncol=1) # d.p.=0,005707
lambda <- 2.135 # d.p.=0,5988
# #Definindo matrizes e calculando as predicoes ’y_hat’ #
x <- matrix(c(ozonio.dat$concentracao),nrow=namos,ncol=1)
X <- x - mean(x)
y <- direcao
Xbeta <- X%*%beta
#A <- sapply(lambda,Bessel.FUN)
#dA <- 1 - A/lambda - A*A

#OBTENDO VALORES PREDITOS EM (0, 2*pi)

 72

yhat <- array(data=mu + 2*atan(Xbeta) + 2*pi, dim=namos, dimnames=NULL)

for(i in 1:namos)
{
yhat[i] <- mod.ab(yhat[i], 2*pi)
# if (yhat[i] > pi) yhat[i] <- yhat[i] - 2*pi
}
yhat<-matrix(c(yhat),nrow=namos,ncol=1)

#CALCULO DOS RESIDUOS: ’r’ e ’d.estrela’ # Formula Original:

d_i<-2*sign(y-yhat)*sqrt(lambda)*abs(sin((y-yhat)/2)) # comp.desvio
<- 2*sqrt(lambda)*sin((y-yhat)/2) # A expressao a seguir esta’ com
a correcao # proposta de Davison & Gigli (1989) para # obter uma
aproximacao para a normal N(0,1). #
r <- matrix(NA,nrow=namos,ncol=1)
for(i in 1:namos)
{
r[i] <- DAVGIGLI.FUN( lambda ) * sin( (y[i]-yhat[i])/2 )
}

gvet <- dg(Xbeta)

G <- diag(c(gvet),nrow=namos)
# K <- lambda * aux1 # K <- diag(c(K), nrow=namos)
G2 <- G %*% G
#H representa a matriz ’hat’
H <- G %*% X %*% solve(t(X) %*% G2 %*% X)%*% t(X) %*% G
Hdiag <- diag(H)
d.estrela <- 2*sqrt(lambda)*sin((y-yhat)/2) / sqrt(1-Hdiag)
#FIM DO CALCULO DOS RESIDUOS DO MODELO AJUSTADO

# SIMULACAO DAS ’BANDAS DE CONFIANCA’ #

residuos.r <- matrix(NA,namos,19)

residuos.d <- matrix(NA,namos,19)
y <- array (data=NA,dim=namos,dimnames=NULL)
mu.fit <- yhat
lambda.fit <- lambda

#INICIO DO CALCULO DAS ’BANDAS’

 73

for(j in 1:19)
{
#USANDO O ALGORITMO DE BEST & FISHER PARA GERAR ’y’
for (i in 1:namos)
{
y[i] <- rvm(1,mu.fit[i],lambda) #BestFisher.FUN(mu.fit[i], lambda.fit)
# if (y[i] < -pi) {y[i] <- y[i] + 2*pi}
# if (y[i] > pi) {y[i] <- y[i] - 2*pi}
}
y<-matrix(c(y),nrow=namos)

### FAZENDO O AJUSTE DA AMOSTRA SIMULADA # # #

* # * Ajuste do Modelo de MEDIAS * #

Inicio do algoritmo GERAL

### Estimativa inicial de ’mu’

Fisher (1993, eq. 2.9)
’mu’= media circular (origem) ###

Sen0 <- sum(sin(y))/namos

Cos0 <- sum(cos(y))/namos
mu <- arctg(Sen0, Cos0)
if (mu > pi) mu <- mu-2*pi #Transf. ’mu’ para o intervalo (-pi , pi)

### Inicializacao do vetor beta ###

X <- x
beta.ini <- solve(t(X)%*%X) %*% t(X) %*% tan((y-mu)/2)
ajuste.amostra.simulada <-
lm.circular.cl(y , X, init=beta.ini, verbose=FALSE)#, tol=1e-10)
beta<-ajuste.amostra.simulada$coeff

#CALCULO DAS ESTIMATIVAS ’y_hat’ DA ’j-ESIMA’ AMOSTRA SIMULADA

#OBTENDO OS VALORES PREDITOS EM (0, 2*pi) E TRANSFORMANDO PARA (-pi,
pi)
Xbeta <- X%*%beta
yhat <- array(data=mu + 2*atan(Xbeta) + 2*pi, dim=namos,dimnames=NULL)
for(i in 1:namos)
 74

{
yhat[i] <- mod.ab(yhat[i], 2*pi)
# if (yhat[i] > pi) yhat[i] <- yhat[i] - 2*pi
}
yhat<-matrix(c(yhat),nrow=namos,ncol=1)

#CALCULO DOS RESIDUOS: ’r’ e ’d.estrela’ DA ’j-ESIMA’ AMOSTRA

SIMULADA

r.aux <- matrix(NA,nrow=namos,ncol=1)

for(i in 1:namos)
{
r.aux[i] <- DAVGIGLI.FUN( lambda ) * sin( (y[i]-yhat[i])/2 )
}

gvet <- dg(Xbeta)

G <- diag(c(gvet),nrow=namos)
G2 <- G %*% G
# K <- lambda * aux1 # K <- diag(c(K), nrow=namos)

#H representa a matriz ’hat’

H <- G %*% X %*% solve(t(X) %*% G2 %*% X)%*% t(X) %*% G
Hdiag <- diag(H)
d.aux <- 2*sqrt(lambda)*sin((y-yhat)/2) / sqrt(1-Hdiag)
#FIM DO CALCULO DOS RESIDUOS DA ’j-ESIMA’ AMOSTRA SIMULADA

residuos.r[,j] <- sort( r.aux )

residuos.d[,j] <- sort( d.aux )
} #FIM DO CALCULO DOS RESIDUOS DA ’j-ESIMA’ AMOSTRA SIMULADA

#FIM DOS CALCULOS PARA OBTENCAO DAS ’BANDAS’ (item 7)

### INICIO DO ’LOOP-2’ (item 8) . # .

k<-19 alfa<-0.05
linf.r<-numeric(namos) lsup.r<-numeric(namos)
linf.d<-numeric(namos)
lsup.d<-numeric(namos)
 75

alfa1<-ceiling(k*alfa) alfa2<-ceiling(k*(1-alfa)) for(i in 1:namos)

{
eo<-sort(residuos.r[i,])
linf.r[i]<-eo[alfa1]
lsup.r[i]<-eo[alfa2]
eo<-sort(residuos.d[i,])
linf.d[i]<-eo[alfa1]
lsup.d[i]<-eo[alfa2]
}

xb.r<-apply(residuos.r,1,mean)
xb.d<-apply(residuos.d,1,mean)

#GRAFICO DO RESIDUO ’r’ ISOLADO win.graph()

faixa<-range(r,linf.r,lsup.r) par(pty="s")
qqnorm(linf.r,axes=F,xlab="",ylab="",ylim=faixa,type="l",lty=1,lwd=1)
par(new=T)
qqnorm(lsup.r,axes=F,xlab="",ylab="",ylim=faixa,type="l",lty=1,lwd=1)
par(new=T) qqnorm(xb.r,xlab="Quantiles of Standard
Normal",ylab="Deviance Residual",ylim=faixa,type="l",lty=2,lwd=1)
text(-1.2,1.2,"Residual r") par(new=T) #oldpar
<-par(pch=15,cex=0.1,csi=0.1) # set new values, remember old oldpar
<-par(pch=16,mkh=0.1) # set new values, remember old
qqnorm(r,axes=F,xlab="",ylab="",ylim=faixa)

par(oldpar) # set parameters back to remembered values #title("(2a)") #

#GRAFICO DO RESIDUO ’d.estrela’ ISOLADO win.graph()

faixa<-range(d.estrela,linf.d,lsup.d) par(pty="s")
qqnorm(linf.d,axes=F,xlab="",ylab="",ylim=faixa,type="l",lty=1,lwd=1)
par(new=T)
qqnorm(lsup.d,axes=F,xlab="",ylab="",ylim=faixa,type="l",lty=1,lwd=1)
par(new=T) qqnorm(xb.d,xlab="Quantiles of Standard
Normal",ylab="Deviance Residual",ylim=faixa,type="l",lty=2,lwd=1)
text(-1.2,1.2,"Residual d*") par(new=T) #oldpar
<-par(pch=15,cex=0.1,csi=0.1) # set new values, remember old oldpar
<-par(pch=16,mkh=0.1) # set new values, remember old
qqnorm(d.estrela,axes=F,xlab="",ylab="",ylim=faixa) par(oldpar) #
set parameters back to remembered values #title("(b)")
 76

#qqnorm(d.estrela,xlab="Percentis da N(0,1)",ylab="Residuo
d*",ylim=faixa) #title("Grafico dos envelopes para o modelo von
Mises") #

#GRAFICO DOS RESIDUOS ’r’ e ’d.estrela’ CONJUNTAMENTE #

win.graph() par(mfrow=c(1,2))
#PRIMEIRO: ’r’
faixa<-range(r,linf.r,lsup.r) par(pty="s")
qqnorm(linf.r,axes=F,xlab="",ylab="",ylim=faixa,type="l",lty=1,lwd=1)
par(new=T)
qqnorm(lsup.r,axes=F,xlab="",ylab="",ylim=faixa,type="l",lty=1,lwd=1)
par(new=T) qqnorm(xb.r,xlab="Quantiles of Standard
Normal",ylab="Deviance Residual",ylim=faixa,type="l",lty=2,lwd=1)
text(-1.2,1.2,"Residual r")

par(new=T)#
oldpar<-par(pch=15,cex=0.1,csi=0.1) # set new values, remember old
oldpar<-par(pch=16,mkh=0.1) # set new values, remember old
qqnorm(r,axes=F,xlab="",ylab="",ylim=faixa)
par(oldpar) # set parameters back to remembered values #
title("(2a)")#
qqnorm(r,xlab="Percentis da N(0,1)",ylab="Residuo r",ylim=faixa)
title("Grafico dos envelopes para o modelo von Mises") #
#SEGUNDO:’d.estrela’
faixa<-range(d.estrela,linf.d,lsup.d) par(pty="s")
qqnorm(linf.d,axes=F,xlab="",ylab="",ylim=faixa,type="l",lty=1,lwd=1)
par(new=T)
qqnorm(lsup.d,axes=F,xlab="",ylab="",ylim=faixa,type="l",lty=1,lwd=1)
par(new=T)
qqnorm(xb.d,xlab="Quantiles of Standard
Normal",ylab="Deviance Residual",ylim=faixa,type="l",lty=2,lwd=1)
text(-1.2,1.2,"Residual d*") par(new=T)
oldpar<-par(pch=15,cex=0.1,csi=0.1) # set new values, remember old
oldpar<-par(pch=16,mkh=0.1) # set new values, remember old
qqnorm(d.estrela,axes=F,xlab="",ylab="",ylim=faixa)
par(oldpar) #set parameters back to remembered values #
title("(b)")
qqnorm(d.estrela,xlab="Percentis da N(0,1)",ylab="Residuo
d*",ylim=faixa)
 77

title("Grafico dos envelopes para o modelo von Mises")

#-----------------------------------------------------------------------#

### Este programa faz o gráfico da Concentração do Ozônio versus

Direção do Vento ###

#LEITURA DOS DADOS

ozonio.dat <-read.table("ozonio.dat", header=TRUE)

attach(ozonio.dat)

##1◦ Gráfico##
plot(ozonio.dat$concentracao, ozonio.dat$direcao,
xlim=c(0,120), ylim=c(0,720), xlab="Concentração de Ozônio",
ylab="Direção do Vento",axes=F, pch=1) par(new=TRUE)
plot(ozonio.dat$concentracao, ozonio.dat$direcao+360, xlim=c(0,120),
ylim=c(-10,720), xlab="", ylab="",axes=F, pch=16)
axis(side=1,c(0,20,40,60,80,100,120))
axis(side=2,c(0,90,180,270,360,450,540,630,720)) box()
title(main = "Gráfico . . .")
# Cálculos auxiliares
namos <-length(ozonio.dat$direcao)
direcao <-matrix(c(ozonio.dat$direcao),nrow=namos,ncol=1)

require(circular) teta<-rad(ozonio.dat$direcao) teta<- -teta + pi/2

abcissa<-ozonio.dat$concentracao*cos(teta)
ordenada<-ozonio.dat$concentracao*sin(teta)

win.graph()
##2◦ Gráfico##
plot(abcissa,ordenada,xlim=c(-50,90),ylim=c(-50,110),xlab="",
ylab="", axes=F) axis(side=1,c(0)) axis(side=2,c(0)) par(new=T)
plot.circular(as.circular(teta), stack=TRUE, bins=150, shrink=0.025,
type="n", axes=F,xlim=c(-50,90),ylim=c(-50,110))

box() points(0, 0, pch="+") #title(main = "Gráfico . . .")

Bibliografia

[1] Abramowitz, M. e Stegun, I. A. (1970). Handbook of Mathematical Functions. New

York: Dover Pub.

[2] Batschelet, E. (1981). Circular Statistics in Biology. London: Academic Press.

[3] Best, D. J. e Fisher, N. I. (1981). The bias of the maximum likelood estimators of
the von Mises-Fisher concentration parameters. Communications in Statistics
- Simulation and Computations 10, 493-502.

[4] Bingham, M. S. e Mardia, K. V. (1975). Maximum likelihood characterization of

the von Mises distribution. In: Statistical Distributions in Scientific Work (C.
P. Patial et al., eds), Vol. 3, pp.387-398. Reidel, Dordrecht.

[5] Cox, D. R. e Snell, E. J. (1968). A general definition of residuals (with discussion).

Journal of the Royal Statistical Society B 30, 248-275.

[6] Davison, A. C. e Gigli, A. (1989). Deviance residuals and normal scores plots.
Biometrika 76, 211-221.

[7] Feller, W. (1966). An Introduction to Probability Theory and its Applications,

Vol.2 Wiley, New York. (27, 54,57, 64).

[8] Fisher, N. I. (1993). Statistical Analysis of Circular Data. Cambridge: University

Press.

[9] Fisher, N. I. e Lee, A. J. (1992). Regression models for an angular response.

Biometrics 48, 665-667.

[10] Gould, A. L. (1969). A regresion technique for angular variates. Biometrics 25,
683-700.
 79

[11] Green, P. J. (1984). Iteratively reweighted least squares for maximum likelihood
estimation, and some robust and resistant alternatives. Journal of the Royal
Statistical Society B 46, 2, 149-192.

[12] Johnson, R. A. e Wehrly, T. E. (1978). Some angular-linear distributions and

related regression models. Journal of the American Statistical Association 73,
602-606.

[13] Laycock, P. J. (1975). Optimal desing: Regression models for directions. Biome-
trika 62, 305-311.

[14] Mardia, K. V. (1972). Statistics of Directional Data. London: Academic Press.

[15] Mardia, K. V. e Jupp, P. E. (1999). Directional Statistics. Wiley series in proba-

bility and statistics.

[16] Marshal,A. W. e Olkin, I. (1961). Game theoretic proof that Chebyshev inequali-
ties are sharp. Pacific J. Math., 11, 1421-1429. (31).

[17] McCullagh, P. (1987). Tensor Methods in Statistics. New York: Chapman and
Hall.

[18] McCullagh,P. e Nelder, J. A. (1989). Generalized Linear Models. 2nd edition. Lon-
don: Chapman and Hall.

[19] Paula, G. A. (2004). Modelos de Regressão com Apoio Computacional. Instituto

de Matemática e Estatística - USP.

[20] Souza, F. A. M. (1999). Influência Local e Análise de Resíduos em Modelos de

Regressão von Mises. Tese de Doutorado. Instituto de Matemática e Estatística
- USP.

[21] Souza, F. A. M. e Paula, G. A. (2002). Deviance Residuals for an Angular Res-

ponse. Australian Statistical Publishing Association Inc. 2002 44, 345-356.

[22] Stephens, M. A. (1963). Random walk on a circle. Biometrika 50, 385-390.

[23] Williams, D. A. (1984). Residuals in generalized linear models. In: Proceedings of

the 12th. International Biometrics Conference, Tokyo, pp. 59-68.

[24] Zygmund, A. (1959). Trigonometrics Series. Cambridge University Press, Vol.1,

Cambridge. (27)