I apologize for the poor R coding and can't guarantee its all correct. It also has a couple of modifications .
#Entaglement, Bells inequality, and hidden variables
#Simulation of the experiment described by N. David Mermin in Physics Today April 1985 pages 38-47.
# http://lilith.fisica.ufmg.br/~fqii/Mermin-PhysToday85.pdf
#Two detectors A and B and one source C.
#The detectors have 3 settings
#The source produces 2 "particles" one received by A and the other received by B
#The "particles" are identical and the detectors flash red or green
#Here the "particles" are also identical, but the detectors record 1 or -1
#The particles are described by an instruction that provides what will be recorded by the detector for each setting
#The possible instructions are the set of full possible combinations or a subset
#The setting for each detector is random and independent of the other detector
#Expected results from Quantum mechanics experiments
# 1) When the settings are the same they flash the same colors or in this case the detector measurements are the same [1,1] or [-1,-1]
# 2) When the runs are evaluated independent of the settings of the detectors the pattern of the flashing (light colors) is random.
# Half the time lights flash the same color and half the time they flash different colors.
#Result (1) is a consequence of the experiment using the same instruction and can simply be explained by a common source making the two "particles" completely correlated
#Result (2) Bells inequal states that in this type of instruction set experiment, the same colors flash at least 5/9 of the time.
#
#Modifications
#p_setting probably the setting is changed to be the same for A and B, or alternatively the instruction is changed to [1,1,1] or [2,2,2]
#p_random probability the result is random independent of the instruction, i.e., the instruction is not used
#p_different probability the results are changed to be different
#
#set up all possible combinations of the instructions
state_instruction <- matrix(NA,8,3)
state_instruction[1,] <- c(1,1,1)
state_instruction[2,] <- c(1,1,-1)
state_instruction[3,] <- c(1,-1,1)
state_instruction[4,] <- c(1,-1,-1)
state_instruction[5,] <- c(-1,1,1)
state_instruction[6,] <- c(-1,1,-1)
state_instruction[7,] <- c(-1,-1,1)
state_instruction[8,] <- c(-1,-1,-1)
n <- 1000000
tt <- Entanglement(n,state_instruction,1:8,p_setting=0,p_random=0,p_different=0)
#tt <- Entanglement(n,state_instruction,1:8,p_setting=0,p_random=0,p_different=2.25/9)
#tt <- Entanglement(n,state_instruction,c(2,4,6),p_setting=0,p_random=0,p_different=0) #Use a subset of the instructions
#Report the results
#verify result (1)
#The number of runs when the settings are the same between detectors A and B
#dim(tt$state[tt$A_setting[]==tt$B_setting[],])
#The number of runs where the measurements are the same when the detector settings are the same
#tt.equal <-tt$state[tt$A_setting[]==tt$B_setting[],]
#dim(tt.equal[tt.equal[,1]==tt.equal[,2],])[1]
#The number of runs where the measurements are different when the detector settings are the same
#dim(tt.equal[tt.equal[,1]!=tt.equal[,2],])[1]
#Report the results
#verify result (2)
#The number of runs where the measurements are the same when the detector settings are the same
dim(tt$state[tt$state[,1]==tt$state[,2],])[1]
#The number of runs where the measurements are different when the detector settings are the same
dim(tt$state[tt$state[,1]!=tt$state[,2],])[1]
#Proportion the same
dim(tt$state[tt$state[,1]==tt$state[,2],])[1]/n
#minimum
5/9
Entanglement <- function(n=1000000,state_instruction,which_instructions,p_setting, p_random, p_different){
n_settings <- dim(state_instruction)[2]
n_instructions <- dim(state_instruction)[1]
#Run the simulation
#Create the settings for each detector
A_setting <- sample(1:n_settings,n, replace=T)
B_setting <- sample(1:n_settings,n, replace=T)
#Creat the instructions
instructions <- sample(which_instructions,n, replace=T) #Use all the possible instructions
#change the result based on hidden variable
setting_change <- runif(n)
B_setting[setting_change <=p_setting] <- A_setting[setting_change <=p_setting]
no_instruction <- runif(n)
different <- runif(n)
#Determine the measurement of each detector
state <- matrix(NA,n,2)
for(i in 1:n){
state[i,1] <- state_instruction[instructions[i],A_setting[i]] #state for detector A in sample i
state[i,2] <- state_instruction[instructions[i],B_setting[i]] #state for detector B in sample i #Multiply this by -1 to get spin up spin down
if (no_instruction[i] <= p_random) {
state[i,1] <- sample(c(1,-1),1, replace=T)
state[i,2] <- sample(c(1,-1),1, replace=T)
}
if (no_instruction[i] <= p_different) {
state[i,2] <- -1*state[i,1]
}
}
return(list(state=state,A_setting=A_setting,B_setting=B_setting))
}