Skip to content

Instantly share code, notes, and snippets.

@mikeguggis

mikeguggis/script.R

Last active Jul 28, 2018
Embed
What would you like to do?
Stack Exchange answer PR(X<Y|min(X,Y))
library(MASS)
library(mvtnorm)
# Set parameters for generating simulated data
set.seed(1)
mu = c(.1,1.5)
Sigma = matrix(c(1.5,.5,.5,.7),2,2)
# Set minumum
m = 1
# Set half of conditional set size
eps = .01
storage = matrix(, nrow = 0, ncol = 2)
while(dim(storage)[1] < 1000){
X = mvrnorm(100000,mu,Sigma)
# Calculate pairwise minimums
mins = pmin(X[,1],X[,2])
# Create subset for conditional statement
dex = mins < m + eps & mins > m - eps # min(x,y) = m
storage = rbind(storage,X[dex,])
}
# Calculate conditional distribution parameters
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 = (1 - Sigma[1,2]^2/(Sigma[2,2]*Sigma[1,1]))*Sigma[1,1]
s2yx = (1 - Sigma[1,2]^2/(Sigma[1,1]*Sigma[2,2]))*Sigma[2,2]
# Check it
out1 = mean(storage[,1]<storage[,2])
out2 = dnorm((m-mu[1])/sqrt(Sigma[1,1]))/sqrt(Sigma[1,1]) * (1-pnorm((m-muyx)/sqrt(s2yx))) /
(dnorm((m-mu[1])/sqrt(Sigma[1,1]))/sqrt(Sigma[1,1]) * (1-pnorm((m-muyx)/sqrt(s2yx))) +
dnorm((m-mu[2])/sqrt(Sigma[2,2]))/sqrt(Sigma[2,2]) * (1-pnorm((m-muxy)/sqrt(s2xy))))
c(out1,out2)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment