使用雙變量橢圓尋找覆蓋

Question

我在估計問題上掙扎。在幾個例子中，我已經展示瞭如何使用由二元正態分佈生成的點向量來計算雙變量橢圓。代碼工作正常，但覆蓋率（生成的或真實的ps（p1，p2）包含在估計的橢圓中的次數）我看起來似乎非常低。我還應該說，與新版本相比，R的舊版本給出了顯着不同的結果。我現在使用R 3.0.1。這裏是能夠重現問題的代碼。

library(MASS) 
    set.seed(1234) 
    x1<-NULL 
    x2<-NULL 
    k<-1 
    Sigma2 <- matrix(c(.72,.57,.57,.46),2,2) 
    Sigma2 
    rho <- Sigma2[1,2]/sqrt(Sigma2[1,1]*Sigma2[2,2]) 
    eta<-replicate(300,mvrnorm(k, mu=c(-1.01,-2.39), Sigma2)) 
    p1<-exp(eta)/(1+exp(eta)) # true p's 
    n<-60 
    x1<-replicate(300,rbinom(k,n,p1[1,])) 
    x2<-replicate(300,rbinom(k,n,p1[2,])) 

    rate1<-x1/60 # Estimated p's 
    rate2<-x2/60 
    library(car) 
    ell <- dataEllipse(rate1, rate2, levels=c(0.05, 0.95)) 
    library(sp) 
    within<-point.in.polygon(p1[1,], p1[2,], ell$`0.95`[,1], ell$`0.95`[,2]) 
    mean(within) # coverage 

Answer 1

的錯誤在於行：

x1<-replicate(300,rbinom(k,n,p1[1,])) 
x2<-replicate(300,rbinom(k,n,p1[2,])) 

k=1因爲，呼叫rbinom(k,n,p1[1,])生成單個隨機偏離，並且僅使用在p1[1,]所述第一概率。你正在複製這個呼叫300次，所以每個偏差都使用相同的概率。因此，與p1相比，rate1和rate2佔用的參數空間要小得多。

x1<-replicate(300,rbinom(k,n,p1[1,])) 
x2<-replicate(300,rbinom(k,n,p1[2,])) 

rate1<-x1/60 # Estimated p's 
rate2<-x2/60 
library(car) 
plot.new() 
ell <- dataEllipse(rate1, rate2, levels=c(0.05, 0.95), plot.points=T, pch=NA) 
library(sp) 
within<-point.in.polygon(p1[1,], p1[2,], ell$`0.95`[,1], ell$`0.95`[,2]) 
mean(within) 

plot(p1[1,which(within==1)], p1[2,which(within==1)], col="blue", ylim=c(0,1),xlim=c(0,1)) 
points(p1[1,which(within==0)], p1[2,which(within==0)], col="green") 

ell <- dataEllipse(rate1, rate2, levels=c(0.05, 0.95), plot.points=T, pch=NA, add=T) 

正確的代碼給出了適當的覆蓋率（95％左右）：通過您的數據橢圓繪製p1想象這

x1<-rbinom(300,n,p1[1,]) 
x2<-rbinom(300,n,p1[2,]) 
rate1<-x1/60 # Estimated p's 
rate2<-x2/60 
library(car) 
plot.new() 
ell <- dataEllipse(rate1, rate2, levels=c(0.05, 0.95), plot.points=T, pch=NA) 
library(sp) 
within<-point.in.polygon(p1[1,], p1[2,], ell$`0.95`[,1], ell$`0.95`[,2]) 
mean(within) 

plot(p1[1,which(within==1)], p1[2,which(within==1)], col="blue", ylim=c(0,1),xlim=c(0,1)) 
points(p1[1,which(within==0)], p1[2,which(within==0)], col="green") 
ell <- dataEllipse(rate1, rate2, levels=c(0.05, 0.95), plot.points=T, pch=NA, add=T)