Draw some random variables

Let’s draw from a bivariate nomral distribution and construct the minimum of the two variables:

#####################################################################################
## Load libraries ##
#####################################################################################
library(MASS)
library(stats4)
library(tidyverse)
library(broom)
library(scales)
library(viridis)
library(forcats)

# Set seed
set.seed(11062017)
m1<-  0                  # mean of y1
m2<-  1                  # mean of y2
sigma1=1                 # sigma y1
sigma2=0.5               # sigma y2
r1=0                     # sigma_y1_y2
rho<-r1/(sigma1*sigma2)  # correlation
Sigma=matrix(c(sigma1^2,r1,r1,sigma2^2),2,2)
N=1000
y=mvrnorm(N,c(m1,m2),Sigma)
z=apply(y,1,FUN=min)
df<-data.frame(id=1:N,z=z,y1=y[,1],y2=y[,2])

# Function for the distribtuoin of Z = min( y1, y2)
fmin<-function(x,mu1,mu2,sigma1,sigma2,rho){
  f1<- dnorm(x,mu1,sigma1)*pnorm(rho*(x-mu1) / (sigma1*sqrt(1-rho*rho))-
                                   (x-mu2) / (sigma2*sqrt(1-rho*rho))) 
  f2<- dnorm(x,mu2,sigma2)*pnorm(rho*(x-mu2) / (sigma2*sqrt(1-rho*rho))-
                                   (x-mu1) / (sigma1*sqrt(1-rho*rho))) 
  return(f1+f2)
}


# plot a draw
# use gather so the plots can be overlaid
ggplot(data=df %>% gather(variable,value,-id), aes(x=value,color=variable,fill=variable))+
  geom_density(alpha=0.1)+
  stat_function(fun=fmin, args=list(m1,m2,sigma1,sigma2,rho),color="black",linetype=2,size=0.75)+
  theme_minimal()+
  scale_color_viridis(discrete=T,option="C",end=0.8)+
  scale_fill_viridis(discrete=T,option="C",end=0.8)+
    labs(title="Distribution of bivariate normal variables (y1,y2) and z=min(y1,y2)",
         subtitle="shaded areas kernel fit to draw of N=1,000,\nblack line theoretical distribution of z = min(y1,y2)")

We might be interested in seeing how well we can recover the means \(\mu_1\) and \(\mu_2\). We can do this via maximum likelihood.

#log likelihood function
ll<- function(params) {
  params <- as.list(params)
  R<-fmin(z,params$mu1,params$mu2,params$s1,params$s2,params$rho1)
  return(-sum(log(R)))
}

# maximum likelihood
out<-optim(par= c(mu1 = 0,mu2=1, s1=1, s2=0.5,rho1=0),
            fn=ll,
            hessian=T,
            control=list(fnscale=1) )

# combine estimated results with true parameter values
# stack up with broom::tidy
out2<-cbind(tidy(out),se.param=diag(solve(out$hessian)),true=c(m1,m2,sigma1,sigma2,rho))

knitr::kable(out2,row.names=FALSE)

Finally, we can use purrr to run some simulations for different sample sizes.

# function for returning a data frame of size N
draw_samp<-function(N,m1,m2,Sigma){
  y=mvrnorm(N,c(m1,m2),Sigma)
  z=apply(y,1,FUN=min)
  return(data.frame(id=1:N,z=z,y1=y[,1],y2=y[,2]))
}

# draw samples of 10, 100, 1,000 and 10,000
df0<-data.frame(N=c(10,100,1000,10000)) %>% 
  mutate(df=map(N,draw_samp,m1,m2,Sigma)) %>% unnest(df)

# make a plot
ggplot(data=df0 %>% gather(variable,value,-id,-N), aes(x=value,color=variable,fill=variable))+
  geom_density(alpha=0.1)+
  stat_function(fun=fmin, args=list(m1,m2,sigma1,sigma2,rho),color="black",linetype=2,size=0.75)+
  theme_minimal()+
  scale_color_viridis(discrete=T,option="C",end=0.8)+
  scale_fill_viridis(discrete=T,option="C",end=0.8)+
  facet_wrap(~forcats::fct_reorder(paste("N =",comma(N)),N))+
  labs(title="Distribution of bivariate normal variables (y1,y2) and z=min(y1,y2)",
       subtitle="shaded areas kernel fit to draw of N=1,000,\nblack line theoretical distribution of z = min(y1,y2)")

Ground-breaking? No, considering these results were published before I was born. But still, could be handy.