Simulating a dependent binary/logit random variable with a particular mean probability

Question

Here is the simplest statement of the problem: I need to generate a binary random variable (RV) with a given probability. However, this binary RV is dependent on another normal RV.

I will make it more concrete: Assume that $X \sim Normal(20, 2)$ and that $Y$ is a binary RV such that:

$logit(p_i) = b_1*x_{1i} + \alpha$, where $p_i$ is the probability of success in a Bernoulli trial.

Further, the $E(Y)=0.4$. What should be $\alpha$ as a function of the $b_1$ and $E(x_{1i})$ to ensure that when I draw the binary RV, its mean is 0.4???

Here is what I am trying to do (albeit unsuccessfully) in R! code, with $b_1=4$:

inv.logit<-function(l) exp(l)/(1+exp(l)); # helper function;
num_obs = 10000

x<- rnorm(mean=20, sd=2, n=num_obs)

y.raw<- 4*x - mean(4*x) + logit(0.4)
mean(y.raw)

y<-sapply(y.raw, FUN=function(a) rbinom(n=1,size=1,p=inv.logit(a)))
mean(y)

The mean for the simulated RV is not 0.4.

NB: With Gaussian/Normal dependent RVs, it is quite easy to find the adjustment constant using expectations. What is a general way to find such adjustments for binary RVs, count RVs (e.g. Poisson, Negative Binomial)???

Answer 1

How about binary search? You want to find $\alpha$ such that:

\begin{align} 0.4 = \int \frac{exp(4 \cdot x + \alpha)}{1+exp(4 \cdot x + \alpha)}f(x)dx \end{align}

Since you know that the right-hand-side is monotonically increasing in $\alpha$, binary search will give you the answer pretty quickly. The idea of binary search is that you start with an interval wide enough that the answer is inside and then cut the interval in half. Look at the midpoint of the interval to see if the answer is in the upper or lower half. Discard whichever half of the interval does not contain the answer. Bisect again. Etc. This makes the length of the interval decrease as $1/2^i$ where $i$ is the number of times you have bisected.

Denote the value of the right-hand-side $RHS(\alpha)$. The algorithm goes like this:

1. Begin with interval $[a_0,b_0]$ for which $RHS(a_0)<0.4<RHS(b_0)$
2. Evaluate $RHS(\frac{a_i+b_i}{2})$
3. If $RHS(\frac{a_i+b_i}{2})>0.4$, set $a_{i+1}=a_i$ and $b_{i+1}=\frac{a_i+b_i}{2}$
4. If $RHS(\frac{a_i+b_i}{2})<0.4$, set $a_{i+1}=\frac{a_i+b_i}{2}$ and $b_{i+1}=b_i$
5. Goto step 2 unless $b_i-a_i<\text{Tolerance}$

You could replace step 5 with unless $\left|RHS(\frac{a_i+b_i}{2})-0.4 \right|<\text{Tolerance}$ if you want.

A couple of notes. First, I assume that when you say $X \sim N(20,2)$, the 2 is a variance. If you meant standard deviation of 2, then change the sqrt(2) in the R script below to just plain 2. Second, I'm evaluating the RHS via Monte Carlo --- I take the sample mean of $\frac{exp(4 \cdot x + \alpha)}{1+exp(4 \cdot x + \alpha)}$ after drawing $x$ 10,000 times from the requested normal distribution. If you want more precision, you can increase the number of replications beyond 10,000.

In this particular example, the search takes 9 iterations to give the answer (which is $\alpha=-81.50232$) and the RHS is 0.3997855.

Binary search is generally a good way to solve for $\alpha$ given $c$ in an equation like $c=f(\alpha)$ as long as $f(\cdot)$ is monotonic.

Here is some R code to implement. I'm sure the code is inefficient because I'm inept that way:

   # This script written in response to a Cross Validated question
# https://stats.stackexchange.com/questions/347876/simulating-a-dependent-binary-random-variable-with-a-particular-mean-probability
# 
# It finds the right amount to add to a logit model with one RHS variable x so that
# the mean of the LHS variable is 0.4 given an arbitrary distribution of x.
# 
# specifically:  
#   y_star = 4*x + alpha + epsilon
#   y = 1 if y_star > 0 and 0 otherwise
# x has some distribution (in the example normal(20,2))
# epsilon is distributed logistic
# alpha is to be found so that P(y=1) = 0.4
#
# This comment last updated 05/23/2018


set.seed(12344321)

# Draw x:
x <- rnorm(n=10000,mean=20,sd=sqrt(2))

# Set LHS:
LHS <- 0.4

# Write a little function to evaluate P(y=1; alpha)
RHS <- function(alpha=0){
  return(mean(1-1/(1+exp(4*x+alpha))))
}

# Test it
RHS(alpha=-80)

# Now, we have to find a starting interval
# You could do something iterative, but I am going to just
# solve for alpha in the case of x really big and x really small:

top <- log(LHS/(1-LHS)) - 4*quantile(x,probs=0.01)
bot <- log(LHS/(1-LHS)) - 4*quantile(x,probs=0.99)

top <- unname(top)
bot <- unname(bot)

# Verify that the correct alpha is trapped between bot and top
# die if not
RHS(alpha=top)
RHS(alpha=bot)

if( (topval <- RHS(alpha=top)) < LHS) stop("top of initial interval not high enough")
if( (botval <- RHS(alpha=bot)) > LHS) stop("bottom of initial interval not low enough")

# Set tolerance --- this is how close you want RHS to get to 0.4 before it
# is close enough.  Also set MAXIT, the number of iterations after which
# the search gives up.

TOL <- 0.001
MAXIT <- 1000

# Here is the binary search loop
iter <- 0
print(cbind(iter,bot,top,botval,topval))
while(abs( (midval <- RHS(alpha=(bot+top)/2))-LHS)>TOL & iter<MAXIT){
  # if midval > LHS, then midval becomes new top
  if(midval>LHS){
    top <- (bot+top)/2
    topval <- midval
  } 
  # otherwise midval becomes new bottom
  else{
    bot <- (bot+top)/2
    botval <- midval
  }
  iter <- iter + 1
  print(cbind(iter,bot,top,botval,topval))
}

alpha <- (top+bot)/2
RHS <- RHS(alpha=alpha)
cbind(alpha,RHS)