Skip to content
{{ message }}

Instantly share code, notes, and snippets.

# olooney/sim.R

Forked from mikeguggis/sim.R
Created Jul 9, 2018
Stack Exchange
 library(MASS) library(ggplot2) results <- NULL # Generate simulated data for ( seed in 1:30 ) { set.seed(seed+42) mu = c(0,-1+seed/10) Sigma = matrix(c(1.5,-.5,-.5,.7),2,2) #Sigma = matrix(c(0.5,0.2,0.2,1),2,2) X = mvrnorm(10000000,mu,Sigma) # Calculate pairwise minimums mins = pmin(X[,1],X[,2]) # Set minumum #m = 0.5 m = 1 # Create subsets for conditional statements epsilon = 0.0001 # approximate the conditions with a narrow strip dex = mins < m + epsilon & mins > m - epsilon # min(x,y) = m dexx = X[,1] < m + epsilon & X[,1] > m - epsilon # x = m dexy = X[,2] < m + epsilon & X[,2] > m - epsilon # y = m # conditional means muxy = mu+Sigma[1,2]*(m-mu)/Sigma[2,2] muyx = mu+Sigma[1,2]*(m-mu)/Sigma[1,1] # conditional variances s2xy = Sigma[1,1] - Sigma[1,2]^2/Sigma[2,2] s2yx = Sigma[2,2] - Sigma[1,2]^2/Sigma[1,1] # conditional z-scores zxy = (muxy-m)/sqrt(s2xy) zyx = (muyx-m)/sqrt(s2yx) # Check (1) #mean(X[dex,1]>X[dex,2]) practice <- mean(X[dexx,2]>m)/(mean(X[dexx,2]>m)+mean(X[dexy,1]>m)) # Check (7) #theory <- (1-pnorm(zxy))/(2-pnorm(zxy)-pnorm(zyx)) theory <- 1- (pnorm(zxy))/(pnorm(zxy)+pnorm(zyx)) if ( seed == 1 ) { n = 1000 # only show n points on the charts # scatter plot ggplot( data.frame(x=X[1:n,1], y=X[1:n,2]), aes(x,y)) + theme_bw() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("P(X
to join this conversation on GitHub. Already have an account? Sign in to comment