Skip to content

Instantly share code, notes, and snippets.

@mikeguggis

mikeguggis/sim.R

Last active Jul 9, 2018
Embed
What would you like to do?
Stack Exchange
library(MASS)
# Generate simulated data
set.seed(1)
mu = c(1,1.5)
Sigma = matrix(c(1.5,-.5,-.5,.7),2,2)
X = mvrnorm(10000000,mu,Sigma)
# Calculate pairwise minimums
mins = pmin(X[,1],X[,2])
# Set minumum
m = .5
# Create subsets for conditional statements
dex = mins < m + .0001 & mins > m - .0001 # min(x,y) = m
dexx = X[,1] < m + .0001 & X[,1] > m - .0001 # x = m
dexy = X[,2] < m + .0001 & X[,2] > m - .0001 # y = m
# Calculate z values
muxy = mu[1]+Sigma[1,2]*(m-mu[2])/Sigma[2,2]
muyx = mu[2]+Sigma[1,2]*(m-mu[1])/Sigma[1,1]
s2xy = Sigma[1,1] - Sigma[1,2]^2/Sigma[2,2]
s2yx = Sigma[2,2] - Sigma[1,2]^2/Sigma[1,1]
zxy = (m-muxy)/sqrt(s2xy)
zyx = (m-muyx)/sqrt(s2yx)
# Check (1)
mean(X[dex,1]<X[dex,2])
mean(X[dexx,2]>m)/(mean(X[dexx,2]>m)+mean(X[dexy,1]>m))
# Check (7)
(1-pnorm(zxy))/(2-pnorm(zxy)-pnorm(zyx))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment