Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
library(ggplot2)
library(BayesLogit)
N <- 100
X <- data.frame(X1=rnorm(N),X2=rnorm(N))
logistic <- function (x) 1/(1+exp(-x))
nu <- 2*X$X1
pi <- 1/(1+exp(-nu))
X$y <- rbinom(N,1,pi)
ggplot(X,aes(x=X1,y=X2,shape=y)) + geom_point()
Xa <- data.frame(X1=rnorm(20,0.5,0.5),X2=rnorm(20,2,0.3),y=-1)
Xb <- data.frame(X1=rnorm(5,1,0.5),X2=rnorm(5,-2,0.3),y=-1)
Xall <- rbind(X,Xa,Xb)
ggplot(Xall,aes(x=X1,y=X2,color=y)) + geom_point()
r <- logit(X$y, as.matrix(X[c("X1","X2")]))
nua <- r$beta %*% c(0,2)
nub <- r$beta %*% c(-0.25,-2)
pa <- logistic(nua)
pb <- logistic(nub)
c(mean(pa),mean(pb))
plot(pa,pb)
ea <- function(pa,pb,na,nb) {
n <- na+nb
x <- 0
for (ra in 0:na) {
for (rb in 0:nb) {
r <- ra+rb
if (0 < r & r < n) {
a<- (ra*(nb-rb)-rb*(na-ra))/(r*(n-r))
}
else {
a <- 0.5
}
x <- x + dbinom(ra,na,pa)*dbinom(rb,nb,pb)*a
}
}
x
}
amatrix <- function(na,nb) {
a <- matrix(0.0,na+1,nb+1)
n <- na + nb
for (ra in 0:na) {
for (rb in 0:nb) {
r <- ra+rb
if (0 < r & r < n) {
a[ra+1,rb+1] <- (ra*(nb-rb)-rb*(na-ra))/(r*(n-r))
}
else {
a[ra+1,rb+1] <- 0
}
}
}
a
}
ea <- function(pa,pb,na,nb,a) {
da <- dbinom(0:na,na,pa)
db <- dbinom(0:nb,nb,pb)
(da %*% a %*% db)[[1]]
}
mean(pa)-mean(pb)
aa <- function(pa,pb,na,nb) {
a <- amatrix(na,nb)
mapply(ea,pa,pb,na,nb)
}
mean(aa(pa,pb,10,3))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment