0

Eu tenho o seguinte conjunto de dados:

A = c(14, 11, 18,7,12,11,6)
B = c(14,23,17,16,15,20,16)

Gostaria de fazer uma simulação de diversas amostras com distribuição t de student dos dois grupos para comparar quantas vezes a amostra B é em média maior do que a A.

Eu tentei o seguinte código:

n<- 100000
amostra_A <- rt (n, mean(A),sd(A))
amostra_B <- rt (n, mean(B), sd(B))
p.b_superior <- sum(amostra_B > amostra_A)/n

O resultado deu que a amostra B seria maior somente 28% das vezes, o que não faz sentido, tendo em vista que os valores de B são maiores.

Tentei a modificação para uma distribuição normal e os valores ficaram mais próximos do esperado:

n<- 100000
amostra_A <- rnorm (n, mean(A),sd(A))
amostra_B <- rnorm (n, mean(B), sd(B))
p.b_superior <- sum(amostra_B > amostra_A)/n

Nesse caso, a amostra B aparece maior em 87,5% das vezes.

Entretanto, quando faço o t.test para as duas amostras, tenho um p-valor de 1%, o que indica que existe significância estatística para afirmar que as duas amostras são diferentes.

Pelo que li, a simulação de Monte Carlo chegaria a valores pŕóximos do p-valor, o que não ocorreu nem na simulação com uma normal e menos ainda na simulação com a distribuição t-student.

Imagino que no caso da distribuição t-student, eu esteja fazendo alguma coisa errada.

2
  • 1
    Os códigos amostra_A <- rt(n, mean(A),sd(A)) e amostra_B <- rt(n, mean(B), sd(B)) não fazem sentido para mim. De acordo com a documentação da função rt, o comando rt (n, mean(A), sd(A)) significa que serão geradas n amostras independente de uma t de Student com graus de liberdade iguais a mean(A) e parâmetro de não-centralidade igual a sd(A). Isso está correto? Commented 11/08/2023 às 0:46
  • Como o Marcus colocou, a função rt não utiliza informações de média e variância. rt(n, df, ncp) // df - degrees of freedom // ncp - non-centrality parameter δ Commented 11/08/2023 às 2:47

1 Resposta 1

0

Como foi mencionado nos comentários da questão original, a função rt depende de três parâmetros, tamanho da série (n), graus de liberdade (df) e não centralidade (ncp). Se você tivesse uma amostra maior de dados poderia usar alguma função para estimar esses parâmetros, como fGarch::stdFit().

No seu caso, pode usar a função rt para gerar uma distribuição com média zero, basta não definir o parâmetro ncp, e depois somar a média que desejar à série gerada. Para incorporar a variância (varx) pode usar a relação entre esta é os graus de liberdade (df): varx = df/(df-2). Dessa forma pode obter o valor do df para usar como parâmetro.

Varx =1.02  # valor definido para replicar aproximadamente a distribuição dos seus dados
n<- 100000
set.seed(432)
amostra_A <- rt (n, 2*Varx/(Varx-1))+mean(A)
amostra_B <- rt (n, 2*Varx/(Varx-1))+mean(B)
(p.b_superior <- sum(amostra_B > amostra_A)/n)
sapply(list(amostra_A , amostra_B),summary)

# para ggplot
data <- rbind(data.frame(valor=amostra_A,amostra="Amostra A"),
                    data.frame(valor=amostra_B,amostra="Amostra B"))
ggplot(data, aes(valor, fill = amostra)) + geom_density(alpha = 0.2)
> (p.b_superior <- sum(amostra_B > amostra_A)/n)
[1] 0.99995
> sapply(list(amostra_A , amostra_B),summary)
             [,1]     [,2]
Min.     6.375751 12.88123
1st Qu. 10.604696 16.59484
Median  11.291627 17.27721
Mean    11.284271 17.27938
3rd Qu. 11.961175 17.96108
Max.    16.251446 22.83622

inserir a descrição da imagem aqui

2
  • Muito obrigado pelos esclarecimentos. Só não entendi como vc chegou ao valor de Varx =1.02 . Poderia me esclarecer? Commented 31/12/2023 às 12:35
  • A fórmula é uma aproximação para o caso de uma série aleatória da distribuição t com média zero. O valor foi escolhido por tentativa e erro, procurando gerar e replicar a forma da distribuição de seus dados amostrais.
    – Robert
    Commented 10/01 às 10:30

Você deve fazer log-in para responder a esta pergunta.

Esta não é a resposta que você está procurando? Pesquise outras perguntas com a tag .