Estatística: Aplicação ao Sensoriamento Remoto
SER 202 - ANO 2015
Simulação Estocástica
Camilo Daleles Rennó
[email protected]
http://www.dpi.inpe.br/~camilo/estatistica/
Simulação
O que é Simulação?
é um experimento realizado a partir de modelos (reais ou virtuais)
Pode ser:
determinística: as entradas do modelo são fixas e para uma determinada
combinação de valores de entrada o resultado final é sempre o mesmo
estocástica (ou probabilística): o modelo e/ou as entradas incorporam variações
aleatórias de modo que os resultados são diferentes a cada simulação
X1
X2
..
.
Xn
MODELO
Y
(µ, , , , etc)
entrada fixa + modelo determinístico
X1
X2
..
.
Xn
X1
Y
Y
MODELO
X2
MODELO
(µ, , , , etc)
Xn
(µ, , , , etc)
entrada fixa + modelo estocástico
entrada estocástica + modelo determinístico
Para que fazer Simulação?
avaliar propagação de incertezas (quando a solução analítica é inviável)
avaliar cenários futuros (resultados possíveis)
testar a sensibilidade de parâmetros de um modelo
estimar pontualmente ou por intervalo um determinado resultado de um modelo
testar a significância de um resultado num teste de hipótese
2
Simulação de Monte Carlo
É um método que avalia um modelo determinístico através da aleatorização das entradas
deste modelo.
É particularmente útil quando o modelo é complexo, não-linear, ou quando envolve muitos
parâmetros de entrada (com diferentes graus de incerteza), o que dificultaria uma
solução analítica.
Através de um grande número de repetições (acima de 1000), garante-se que
praticamente todas as combinações de entradas sejam avaliadas.
O termo Monte Carlo foi dado em homenagem a roleta, jogo muito popular de Monte
Carlo, Mônaco.
5
4
X ~ Binomial (n = 5; p = 0,4)
x
0
1
2
3
4
5
P(X = x)
7,78%
25,92%
34,56%
23,04%
7,68%
1,02%
3
0
1
2
área de cada fatia é proporcional a
probabilidade do valor correspondente
3
Geração de Números Aleatórios
Originalmente os números aleatórios eram gerados usando dados, roletas, tabelas, etc.
Atualmente os computadores são usados para gerar números chamados pseudo-aleatórios,
que constituem uma sequencia de valores que, embora sejam gerados de forma
determinística, simulam variáveis aleatórias uniformes [0,1] independentes.
Qualquer variável aleatória pode ser simulada a partir de uma variável aleatória uniforme
[0,1] desde que se conheça a função de distribuição acumulada F(x) = P(X  x).
X ~ Binomial (n = 5; p = 0,4)
P(X = x)
7,78%
25,92%
34,56%
23,04%
7,68%
1,02%
P(X  x)
7,78%
33,70%
68,26%
91,30%
98,98%
100,00%
1,0
U = 0,4367
0,8
P(X  x)
x
0
1
2
3
4
5
X=2
0,6
0,4
0,2
0,0
0
1
2
3
X
4
5
4
Geração de Números Aleatórios
Originalmente os números aleatórios eram gerados usando dados, roletas, tabelas, etc.
Atualmente os computadores são usados para gerar números chamados pseudo-aleatórios,
que constituem uma sequencia de valores que, embora sejam gerados de forma
determinística, simulam variáveis aleatórias uniformes [0,1] independentes.
Qualquer variável aleatória pode ser simulada a partir de uma variável aleatória uniforme
[0,1] desde que se conheça a função de distribuição acumulada F(x) = P(X  x).
Sorteio de 8 valores
0,4138
0,9155
0,6218
0,3848
0,1058
0,2763
0,3855
0,8036
Z ~N(0,1) X ~N(10,4)
-0,2177
1,3753
0,3101
-0,2929
-1,2492
-0,5938
-0,2910
0,8546
9,5645
12,7505
10,6202
9,4142
7,5017
8,8124
9,4180
11,7092
1,0
Z 
X  10
~ N 0,1
2
X  2 Z  10
P(Z  z)
u(0,1)
X ~ Normal ( = 10; 2 = 4)
0,75
0,5
0,25
0,0
-4
-2
0
2
4
Y
5
Avaliação das Simulações
• Estimação da Função de Probabilidade através das frequências relativas
observadas (variáveis discretas)
• Métricas de tendência central e de dispersão:
média, desvio padrão, mediana, quantis, amplitude, mínimo/máximo, etc
• Intervalos de Credibilidade
os limites são definidos, desprezando-se os valores extremos (mesma proporção
para ambos os lados)
• Box-plot
mediana, 1o e 3o quartis e valores extremos (outliers)
6
Boxplot
É uma ótima alternativa para mostrar graficamente a dispersão de observações
de uma amostra e são muito úteis para comparar conjuntos de dados pois
causam grande impacto visual e são fáceis de entender.
Há muitas variações de boxplot mas em geral representam:
a) mediana
b)
1o
e
3o
quartis
c) mínimos e máximos
d) “outliers”
Ex: amostra com 20 valores
120
outliers
extremos
100
outliers
80
último ponto
superior
1,5*DIQ
60
40
20
0
1,5*DIQ
3o quartil
mediana
1o quartil
último ponto
inferior
DIQ (distância interquartil)
1,5*DIQ
7
Boxplot
120
a) qual é a distribuição mais simétrica?
D
100
b) qual é a distribuição mais assimétrica?
80
A
60
c) quais as 2 distribuições que mais se
confundem entre si?
40
AeB
20
0
d) quais as 2 distribuições que mais se
A
B
C
D
distinguem entre si?
BeC
8
Exemplos de Aplicações
Exemplo 1:
estimar função de probabilidade de um experimento complexo (urnas)
(ver Simulacao_exemplo1.xls)
Exemplo 2:
simular dados correlacionados (componentes principais)
(ver Simulacao_exemplo2.xls)
Exemplo 3:
determinar o valor crítico de um teste estatístico (KS para duas amostras)
(ver Simulacao_exemplo3.xls)
9
Exemplo 1
I
A
Etapas:
I) Das urnas A e B, sorteia-se uma
bola de cada. As duas bolas são
colocadas na urna C
B
C
10
Exemplo 1
Etapas:
I) Das urnas A e B, sorteia-se uma
bola de cada. As duas bolas são
colocadas na urna C
II) Da urna C, sorteiam-se duas bolas
(sem reposição)
A
B
II
C
11
Exemplo 1
A
Etapas:
I) Das urnas A e B, sorteia-se uma
bola de cada. As duas bolas são
colocadas na urna C
II) Da urna C, sorteiam-se duas bolas
(sem reposição)
III)Se as bolas forem da mesma cor,
ambas são colocadas na urna A.
Caso contrário, ambas são
colocadas na urna B
B
C
Sim
bolas de
mesma
cor?
III
Não
12
Exemplo 1
IV
A
B
50%
FR
7,85%
47,00%
40,19%
4,31%
0,65%
0,00%
40%
P(X = x)
X
0
1
2
3
4
5
30%
20%
10%
0%
Etapas:
I) Das urnas A e B, sorteia-se uma
bola de cada. As duas bolas são
colocadas na urna C
II) Da urna C, sorteiam-se duas bolas
(sem reposição)
III)Se as bolas forem da mesma cor,
ambas são colocadas na urna A.
Caso contrário, ambas são
colocadas na urna B
IV) Escolhe-se aleatoriamente a urna A
ou B e dela retiram-se 5 bolas (sem
reposição)
C
Definindo-se X como o número de bolas
azuis nas 5 observações, qual a
distribuição dos valores de X?
13
Exemplo 1 em R
>A<-c("R","R","R","G","B","B")
>B<-c("R","G","G","G","G","B")
>C<-c("R","G","B","B","B")
>n<-10000
>p<-rep(0,6)
>for (i in 1:n) {
>Af<-A
>Bf<-B
>sorteio1<-c(sample(A,size=1),sample(B,size=1))
>Cf<-c(C,sorteio1)
>sorteio2<-sample(Cf,size=2)
>if (sorteio2[1] == sorteio2[2]) Af<-c(A,sorteio2) else Bf<-c(B,sorteio2)
>if (runif(1,0,1) < 0.5) sorteio3<-sample(Af,5) else sorteio3<-sample(Bf,5)
>nB<-length(which(sorteio3 == "B"))
>p[nB+1]<-p[nB+1]+1
>}
>p<-p/n
>p
[1] 0.0817 0.4702 0.3939 0.0474 0.0068 0.0000
14
Exemplo 2
Simulação de 3 variáveis correlacionadas sendo que
X 1 ~ N (  X 1  10,3;  X2 1  3,5)
X 2 ~ N (  X 2  35,1;  X2 2  7,8)
X 3 ~ N (  X 3  107,8;  X2 3  1,2)
μ  10,3 35,1 107,8
 3,5  2,2  0,5
Σ X   2,2 7,8
0,9 
  0,5 0,9
1,2 
 X ~ N μ, ΣX 
Deseja-se simular 500 valores destas 3 variáveis respeitando suas correlações.
A simulação conjunta, mesmo para uma distribuição gaussiana, não é muito simples.
Por outro lado, a simulação de variáveis independentes é bastante simples de ser
implementada quando se conhece a distribuição destas variáveis:
Y ~ N ( ,  2 )
P(Z  z ) ~ U (0,1)  Y  z  
Solução: para simular variáveis correlacionadas, primeiramente gera-se variáveis
independentes através de uma transformação por componentes principais,
simulam-se estas componentes independentemente uma das outras e depois
aplica-se a transformação inversa.
17
Exemplo 2
 3,5  2,2  0,5
Σ X   2,2 7,8
0,9 
  0,5 0,9
1,2 
λ  8,8629 2,5819 1,0552
TCP

0,1182 
 0,3867 0,9146
α   0,9126 0,3980  0,0934
  0,1324  0,0717 0,9886 
Por definição, todas as componentes principais tem média zero e a variância de
cada componente é dada pelo autovalor correspondente.
Além disso, no caso de distribuições gaussianas, a forma da distribuição é
preservada após a transformação, ou seja, as componentes também têm
distribuição gaussiana:
2
PC1 ~ N (   0;  PC
 1  8,8629)
1
2
PC2 ~ N (   0;  PC
 2  2,5819)
2
2
PC3 ~ N (   0;  PC
 3  1,0552)
3
PC  αX  1(5001) μ '
 
'
X  α 1 PC  1(5001) μ '
α 1  α
X  αPC  1(5001) μ '
18
Exemplo 2
Neste exemplo, deseja-se gerar 500 valores para cada variável X.
O processo inicia-se com a geração de números aleatórios de uma distribuição
uniforme contínua entre 0 e 1:
U1
0.9789
0.7383
0.9224
0.5247
0.9010
.
.
.
0.7961
0.4524
0.5987
0.6906
0.2216
U2
0.1352
0.7855
0.3914
0.0030
0.0007
.
.
.
0.7392
0.9962
0.3815
0.6003
0.5514
U3
0.1669
0.3370
0.2114
0.6637
0.7718
.
.
.
0.0535
0.5041
0.5824
0.7568
0.0853
19
Exemplo 2
A partir de U calcula-se Z através da relação: P(Z  zij )  uij
Exemplo: P(Z  z1,1 )  0,9789  z1,1  2,0312
U1
Z1
2.0312
0.9789
0.7383
0.6381
0.9224
1.4215
0.5247
0.0619
0.9010
1.2872
.
.
.
0.7961
0.8277
-0.1195
0.4524
0.5987
0.2500
0.6906
0.4975
-0.7667
0.2216
U2
Z2
-1.1020
0.1352
0.7855
0.7909
-0.2756
0.3914
-2.7474
0.0030
-3.1950
0.0007
.
.
.
0.7392
0.6408
0.9962
2.6664
-0.3016
0.3815
0.6003
0.2542
0.5514
0.1293
U3
Z3
-0.9663
0.1669
-0.4208
0.3370
-0.8016
0.2114
0.6637
0.4226
0.7718
0.7448
.
.
.
-1.6114
0.0535
0.5041
0.0103
0.5824
0.2081
0.7568
0.6960
-1.3702
0.0853
20
Exemplo 2
A partir de Z calcula-se PC multiplicando cada zij pelo desvio padrão da
componente principal i, ou seja, a raiz quadrada do autovalor i: pcij  zij 
Exemplo: pc1,1  z1,1 
PC1
Z1
6.0470
2.0312
0.6381
1.8996
1.4215
4.2320
0.0619
0.1842
1.2872
3.8320
.
.
.
0.8277
2.4641
-0.1195
-0.3558
0.2500
0.7443
0.4975
1.4812
-0.7667
-2.2824
PC2
Z2
-1.7708
-1.1020
0.7909
1.2708
-0.2756
-0.4428
-2.7474
-4.4145
-3.1950
-5.1337
.
.
.
0.6408
1.0296
2.6664
4.2844
-0.3016
-0.4846
0.2542
0.4085
0.1293
0.2077
i
1  2,0312 8,8629  6,0470
PC3
Z3
-0.9926
-0.9663
-0.4208
-0.4322
-0.8016
-0.8234
0.4226
0.4341
0.7448
0.7651
.
.
.
-1.6114
-1.6553
0.0103
0.0106
0.2081
0.2138
0.6960
0.7149
-1.3702
-1.4075
21
Exemplo 2
Finalmente, aplica-se a rotação inversa definida pela matriz de autovetores e
soma-se a média correspondente a cada uma das variáveis X:
xij  1,1 pc1, j  1,2 pc2, j  1,3 pc3, j  i onde kw é o k-ésimo elemento do autovetor w
Exemplo: x1,1  1,1 pc1,1  1, 2 pc2,1  1,3 pc3,1  1 
 0,3867 6,0470 0,91461,7708 0,1182 0,9926 10,3  10,9015
PC1
X1
Z1
10.9015
2.0312
6.0470
12.1458
0.6381
1.8996
11.4342
1.4215
4.2320
0.0619
0.1842
6.3850
1.2872
3.8320
7.1769
.
.
.
11.9989
0.8277
2.4641
14.0822
-0.1195
-0.3558
10.1698
0.2500
0.7443
11.3309
0.4975
1.4812
-0.7667
-2.2824
9.4410
PC2
X2
Z2
28.9694
-1.1020
-1.7708
33.9125
0.7909
1.2708
31.1386
-0.2756
-0.4428
33.1343
-2.7474
-4.4145
29.4882
-3.1950
-5.1337
.
.
.
33.4157
0.6408
1.0296
37.1289
2.6664
4.2844
34.2079
-0.3016
-0.4846
33.8440
0.2542
0.4085
37.3971
0.1293
0.2077
PC3
X3
Z3
106.1450
-0.9663
-0.9926
107.0300
-0.4208
-0.4322
106.4574
-0.8016
-0.8234
108.5212
0.4226
0.4341
108.4171
0.7448
0.7651
.
.
.
105.7635
-1.6114
-1.6553
107.5504
0.0103
0.0106
107.9476
0.2081
0.2138
108.2814
0.6960
0.7149
106.6958
-1.3702
-1.4075
22
Exemplo 2
É importante notar que, por se tratar de uma simulação, os valores de média,
variância e covariância estimados a partir dos valores simulados não são
exatamente os que se desejava:
μ  10,3 35,1 107,8
 3,5  2,2  0,5
Σ X   2,2 7,8
0,9 
  0,5 0,9
1,2 
X  10,28 35,11 107,83
 3,37  2,16  0,59

ˆ    2,16 8,19
Σ
0
,
93
X


 0,59 0,93
1,40 
23
Exemplo 2 em R
>mX<-c(10.3,35.1,107.8)
>covX<-cbind(c(3.5,-2.2,-0.5),c(-2.2,7.8,0.9),c(-0.5,0.9,1.2))
>auto<-eigen(covX) #gera autovalores e autovetores a partir da matriz de covariância
>n<-500
>U1<-runif(n,0,1) #gera n valores aleatórios de uma distribuição uniforme contínua entre 0 e 1
>U2<-runifn,0,1)
>U3<-runif(n,0,1)
>U<-cbind(U1,U2,U3)
>pairs(U,upper.panel = NULL)
>Z<-qnorm(U) #em R há uma função que gera valores de uma normal padrão: Z<-matrix(rnorm(3*n),nrow=n,ncol=3)
>colnames(Z)<-c("Z1","Z2","Z3")
>pairs(Z,upper.panel = NULL)
>PC<-Z*matrix(rep(1,n))%*%t(sqrt(auto$values))
>colnames(PC)<-c("PC1", "PC2", "PC3")
>pairs(PC,upper.panel = NULL)
>X<-t(auto$vectors%*%t(PC)+t(matrix(rep(1,n))%*%mX))
>colnames(X)<-c("X1", "X2", "X3")
>pairs(X,upper.panel = NULL)
>round(colMeans(X),2)
>round(cov(X))
24
Exemplo 2 em R
X
ˆ 
Σ
X
X1
X2
X3
10.28211 35.19518 107.77577
X1
X2
X3
X1 3.8273337 -2.524998 -0.6134409
X2 -2.5249980 7.808794 1.0919343
X3 -0.6134409 1.091934 1.1782311
25
Exemplo 3
Exemplo Slide 48 (15EstNaoParam.ppt)
Exemplo: Um pesquisador deseja saber se duas regiões de uma mesma imagem
apresentam a mesma distribuição de valores (desconhecida). Para testar esta
hipótese, amostrou-se 15 pontos independentes de cada região. Os valores
observados são apresentados na tabela abaixo. O que se conclui a partir destes
valores?
Valor FRAA FRAB |Dif|
Dobs = 4/15
54
0
1/15
1/15
Região A
81
78
61
89
69
58
64
84
89
83
88
56
87
95
75
Região B
56
55
76
54
83
97
85
66
78
80
61
69
71
55
91
55
0
3/15
3/15
56
1/15
4/15
3/15
58
2/15
4/15
2/15
61
3/15
5/15
2/15
64
4/15
5/15
1/15
66
4/15
6/15
2/15
69
5/15
7/15
2/15
71
5/15
8/15
3/15
75
6/15
8/15
2/15
76
6/15
9/15
3/15
78
7/15
10/15
3/15
80
7/15
11/15
4/15
81
8/15
11/15
3/15
83
9/15
12/15
3/15
84
10/15
12/15
2/15
85
10/15
13/15
3/15
87
11/15
13/15
2/15
88
12/15
13/15
1/15
89
14/15
13/15
1/15
91
14/15
14/15
0
95
15/15
14/15
1/15
97
15/15
15/15
0
KDobs = 4
H0: As duas amostras provêm da mesma
população
H1: As duas amostras provêm de
populações diferentes (bilateral)
Pela simulação:
Valor-P = P(KD > KDobs) = 64,3%
KDcrít 5% = 8
Conclusão: aceita-se H0, ou seja, as duas
amostras provêem da mesma população,
adotando-se 5% de significância
26
Exemplo 3 em R
>regA<-c(81,78,61,89,69,58,64,84,89,83,88,56,87,95,75)
>regB<-c(56,55,76,54,83,97,85,66,78,80,61,69,71,55,91)
>min<-min(regA,regB)
>max<-max(regA,regB)
>dif<-rep(0,max-min+1)
>for (i in min:max) dif[i-min+1]<-abs(length(which(regA <= i))-length(which(regB <= i)))
>KDobs<-max(dif)
>KDobs
>regAB<-c(regA,regB)
>n<-10000
>pKD<-rep(0,16)
>ValorP<-0
>for (k in 1:n) {
>regAB<-sample(regAB)
>regAt<-regAB[1:length(regA)]
>regBt<-regAB[length(regA)+1:length(regAB)]
>dif<-rep(0,max-min+1)
>for (i in min:max) dif[i-min+1]<-abs(length(which(regAt <= i))-length(which(regBt <= i)))
>KD<-max(dif)
>if (KD <= KDobs) ValorP<- ValorP+1
>pKD[KD+1]<- pKD[KD+1]+1 }
>pKDcum <-rev(cumsum(rev(pKD))/n)
>KDcrit<-min(which(pKDcum < 0.05))-1
>KDcrit
>ValorP<-ValorP/n
>ValorP
27