###################### ## MATH5040Assn14.r ## ###################### ## Curtis Miller ## ## 12/10/14 ## ## MATH 5040 ## ## Assignment 14 ## ###################### ## sigma (an element of omega) will be represented with a square matrix ## consisting only of -1 and 1, which I will construct later. s = 25 # size of the matrix ## The S function mentioned in the problem, which is the ## sum of the neighbors of a vertex v (represented by x ## and y) in matrix sigma S <- function(sigma, x, y) { ## I use the formula: ## S(sigma, v) = sum(w~v; sigma(w)) return(sigma[x, y %% s + 1] + sigma[x %% s + 1, y] + sigma[x, (y - 2) %% s + 1] + sigma[(x - 2) %% s + 1,y]) } ## The P function gives the probability of moving from ## sigma to another matrix that differs only at a ## single vertex, w (given by coords. x, y). Beta needs ## to be specified. P <- function(sigma, x, y, beta) { ## I use the formula: ## exp(beta*S(sigma, w))/(exp(beta*S(sigma, w)) + exp(-beta*S(sigma, w))) ## This is the probability of getting a 1 at a certain vertex holding all other vertices constant, as in the Gibbs sampler return(exp(beta * S(sigma, x = v[1], y = v[2])) / (exp(beta * S(sigma, x = v[1], y = v[2])) + exp(-1 * beta * S(sigma, x = v[1], y = v[2])))) } ## This function is for simulating the Glauber dynamics for the Gibbs ## distribution. It takes a sigma, picks a point at random, changes that ## point with probability P(sigma, rho) (rho differs from sigma at that ## point, then returns the resulting matrix. This is meant to be used ## in a loop repeatedly. Beta is needed as well. simulate.Ising <- function(sigma, beta) { # Random vertex v = sample(1:s, size = 2, replace = T) p = P(sigma, x = v[1], y = v[2], beta = beta) sigma[v[1], v[2]] <- sample(c(1, -1), size = 1, prob = c(p, 1-p)) return(sigma) } ## This function has one job: repeat simulate.Ising n times, then ## return the resulting matrix. repeat.Ising <- function(sigma, beta, n, verbose = F) { for (i in 1:n) { sigma <- simulate.Ising(sigma, beta = beta) # Makes verbose if (i / (n / 100) == floor(i / (n / 100)) & verbose == T) cat('On i =', i, '\n') } return(sigma) } ## This function is like repeat.Ising, only instead of returning ## a matrix, it keeps a running total of times the upper qxq ## submatrix is all equal and returns the probability that that ## event occurs. repeat.count <- function(sigma, beta, n, q, verbose = F) { qxq.equal = 0 for (i in 1:n) { sigma <- simulate.Ising(sigma, beta = beta) # Running total of all 1's in the upper qxq corner if (abs(sum(sapply(1:q, function (x) (sum(sigma[x, 1:q]) / q))) / q) == 1) qxq.equal = qxq.equal + 1 # Makes verbose if (i / (n / 100) == floor(i / (n / 100)) & verbose == T) cat('On i =', i, '\n') } return(qxq.equal / n) } ## A function for printing a torus matrix, with blanks for -1 print.sigma <- function(sigma) { # Turn -1's to zeroes, then print with 0's showing up as ' ' print(as.table(sapply(1:s, function (y) sapply(1:s, function (x) if (sigma[x,y] == -1) sigma[x,y] <- 0 else sigma[x,y] <- 1))), zero.print = ' ') } # Run the simulations print.sigma(repeat.Ising(matrix(rep(1, times = s ^ 2), nrow = s, ncol = s), beta = -2/5, n = 1000000)) print.sigma(repeat.Ising(matrix(rep(1, times = s ^ 2), nrow = s, ncol = s), beta = 0, n = 1000000)) print.sigma(repeat.Ising(matrix(rep(1, times = s ^ 2), nrow = s, ncol = s), beta = 2/5, n = 1000000)) ## Here are the results of one run of the code: ## > print.sigma(repeat.Ising(matrix(rep(1, times = s ^ 2), nrow = s, ncol = s), beta = -2/5, n = 1000000)) ## A B C D E F G H I J K L M N O P Q R S T U V W X Y ## A 1 1 1 1 1 1 1 1 1 1 ## B 1 1 1 1 1 1 1 1 1 1 1 1 ## C 1 1 1 1 1 1 1 1 1 1 ## D 1 1 1 1 1 1 1 1 1 1 1 1 1 ## E 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ## F 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ## G 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ## H 1 1 1 1 1 1 1 1 1 1 1 1 ## I 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ## J 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ## K 1 1 1 1 1 1 1 1 1 1 ## L 1 1 1 1 1 1 1 1 1 1 1 1 1 ## M 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ## N 1 1 1 1 1 1 1 1 1 1 1 ## O 1 1 1 1 1 1 1 1 1 1 1 1 1 ## P 1 1 1 1 1 1 1 1 1 1 ## Q 1 1 1 1 1 1 1 1 1 1 1 1 1 ## R 1 1 1 1 1 1 1 1 1 1 ## S 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ## T 1 1 1 1 1 1 1 1 1 1 1 1 1 ## U 1 1 1 1 1 1 1 1 1 1 1 1 ## V 1 1 1 1 1 1 1 1 1 ## W 1 1 1 1 1 1 1 1 1 1 1 1 1 ## X 1 1 1 1 1 1 1 1 1 1 1 1 1 ## Y 1 1 1 1 1 1 1 1 1 1 1 ## > print.sigma(repeat.Ising(matrix(rep(1, times = s ^ 2), nrow = s, ncol = s), beta = 0, n = 1000000)) ## A B C D E F G H I J K L M N O P Q R S T U V W X Y ## A 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ## B 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ## C 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ## D 1 1 1 1 1 1 1 1 1 1 1 ## E 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ## F 1 1 1 1 1 ## G 1 1 1 1 1 1 1 1 1 ## H 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ## I 1 1 1 1 1 1 1 1 1 1 1 1 ## J 1 1 1 1 1 1 1 1 1 ## K 1 1 1 1 1 1 1 1 1 1 ## L 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ## M 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ## N 1 1 1 1 1 1 1 1 1 1 ## O 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ## P 1 1 1 1 1 1 1 1 1 1 ## Q 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ## R 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ## S 1 1 1 1 1 1 1 1 1 1 1 1 1 ## T 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ## U 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ## V 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ## W 1 1 1 1 1 1 1 1 1 ## X 1 1 1 1 1 1 1 1 1 1 1 1 1 ## Y 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ## > print.sigma(repeat.Ising(matrix(rep(1, times = s ^ 2), nrow = s, ncol = s), beta = 2/5, n = 1000000)) ## A B C D E F G H I J K L M N O P Q R S T U V W X Y ## A 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ## B 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ## C 1 1 1 1 1 1 1 1 1 1 1 1 ## D 1 1 1 1 1 1 ## E 1 1 1 1 1 1 1 1 1 ## F 1 1 1 1 1 1 1 1 1 1 ## G 1 1 1 1 1 1 1 1 1 1 ## H 1 1 1 1 1 1 1 1 ## I 1 1 1 1 1 1 1 1 ## J 1 1 1 1 1 1 ## K 1 1 1 1 1 1 1 1 1 1 ## L 1 1 1 1 1 1 1 1 ## M 1 1 1 1 1 ## N ## O 1 1 1 1 1 1 ## P 1 1 1 1 1 1 1 1 1 1 ## Q 1 1 1 1 1 1 1 1 1 1 1 1 ## R 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ## S 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ## T 1 1 1 1 1 1 1 ## U 1 1 1 1 1 1 ## V 1 1 1 1 1 1 1 1 1 1 ## W 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ## X 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ## Y 1 1 1 1 1 1 1 1 1 1 1 1 ## Let' now do the count, and calculate a standard error. ## We will do this by running repeat.Ising 10 times, then ## taking the mean and standard deviation of the sample. ## We will repeat only 10 times, with n = 1000000 m = 10 qxq.results <- sapply(1:m, function (x) repeat.count(matrix(rep(1, times = s ^ 2), nrow = s, ncol = s), beta = 2/5, n = 1000000, q = 5)) # The estimated probability mean(qxq.results) # The standard error sd(qxq.results)