n_steps =10000delta =0.05x =rep(NA, n_steps)x[1] =0.75# initialize# Posterior is proportional to prior * likelihoodpi_x <-function(x) {if (x >0& x <1) dnorm(x, 0.75, 0.05) *dbinom(60, 100, x) else0}for (n in2:n_steps){ current = x[n -1]# propose next state from "jumping" distribution proposed = current +rnorm(1, mean =0, sd = delta)# compute acceptance probability accept =min(1, pi_x(proposed) /pi_x(current))# simulate next state x[n] =sample(c(current, proposed), 1, prob =c(1- accept, accept)) }# display the first few stepsdata.frame(step =1:n_steps, x) |>head(10) |>kbl() |>kable_styling()
step
x
1
0.7500000
2
0.6448298
3
0.6389480
4
0.6389480
5
0.6114167
6
0.6114167
7
0.6388756
8
0.6503970
9
0.6503970
10
0.6427185
# trace plot of first 100 stepsplot(1:100, x[1:100], type ="o", xlab ="step",ylab ="x", main ="First 100 steps")
# simulated values of xhist(x, breaks =50, freq =FALSE,xlab ="theta", ylab ="pi(x|y = 60)", main ="Conditional Distribution")# plot of theoretical conditional density of x# Normal(0.75, 0.05) is approximately Beta(55.5, 18.5)# So conditional given Y = 60 is approximately Beta(55.5 + 60, 18.5 + 40)x_plot =seq(0, 1, 0.0001)lines(x_plot, dbeta(x_plot, 55.5+60, 18.5+40), col ="seagreen", lwd =2)
Multivariable
Attempts 4 free throws and makes 3 in the first game, and attempts 6 free throws and makes 4 in the second game.
# posterior is proportion to product of prior and likelihoodpi_xz <-function(x, z) {if (z >0& x >0& x <1) {dexp(z, rate =1/2.5) *dnorm(x, 0.75, 0.05) *dpois(4, z) *dbinom(3, 4, x) *dpois(6, z) *dbinom(4, 6, x) } else {0 }}pi_xz(0.7, 5) /pi_xz(0.75, 2.5)
[1] 1.641296
n_steps =11000delta =c(0.05, 0.5) # x, zxz =data.frame(x =rep(NA, n_steps),z =rep(NA, n_steps))xz[1, ] =c(0.75, 2.5) # initializefor (n in2:n_steps){ current = xz[n -1, ]# proposed next state proposed = current +rnorm(2, mean =0, sd = delta)# compute acceptance probability accept =min(1, pi_xz(proposed$x, proposed$z) /pi_xz(current$x, current$z))# simulate next state accept_ind =sample(0:1, 1, prob =c(1- accept, accept)) xz[n, ] = proposed * accept_ind + current * (1- accept_ind)}# display the first few stepsdata.frame(step =1:n_steps, xz) |>head(10) |>kbl() |>kable_styling()
step
x
z
1
0.7500000
2.500000
2
0.7862239
2.751374
3
0.7283005
2.306634
4
0.6954620
2.549367
5
0.6445811
3.424299
6
0.6915942
3.870056
7
0.6915942
3.870056
8
0.6929044
4.835510
9
0.6978748
4.782784
10
0.7428304
4.611797
# Trace plot of first 100 stepsggplot(xz[1:100, ] %>%mutate(label =1:100),aes(x, z)) +geom_path() +geom_point(size =2) +geom_text(aes(label = label, x = x +0.005, y = z +0.01)) +labs(title ="Trace plot of first 100 steps")
Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
ggplot(xz, aes(x, z)) +stat_density_2d(aes(fill =after_stat(level)),geom ="polygon", color ="white") +scale_fill_viridis_c()
Random image
I’m not sure this is working correctly.
beta =5L =20Nsteps =10000print_iter =1000image_initial =matrix(sample(c(0, 1), L ^2, replace =TRUE),nrow = L)image_initial =matrix(rep(0, L ^2), nrow = L)image_current = image_initialfor (n in1:Nsteps) {# proposed pixel to flip x_proposed =sample(1:L, 1) y_proposed =sample(1:L, 1)# find neighbors of proposed pixel# interior pixels have four neighbors each# unique, max/min, and which below are just# to correct for corners and edges neighbor_proposed =unique(rbind(c(x_proposed, min(y_proposed +1, L)),c(x_proposed, max(y_proposed -1, 1)),c(min(x_proposed +1, L), y_proposed),c(max(x_proposed -1, 1), y_proposed) ) ) neighbor_proposed = neighbor_proposed[-which( (( neighbor_proposed[1, ] == x_proposed) & (neighbor_proposed[2, ] == y_proposed) )), ]# for current configuration, number of neighbors proposed state matches sum_current =sum(image_current[x_proposed, y_proposed] == image_current[neighbor_proposed])# for proposed configuration, number of neighbors proposed state matches sum_proposed =sum((1- image_current[x_proposed, y_proposed]) == image_current[neighbor_proposed]) accept_prob =min(1, exp(beta * (sum_proposed - sum_current))) color_current = image_current[x_proposed, y_proposed] image_current[x_proposed, y_proposed] =sample(c(1- color_current, color_current), 1,prob =c(accept_prob,1- accept_prob))if ( n %% print_iter ==0 ) {image(1:L, 1:L, t(image_current),xlab ="", ylab ="",zlim =c(0, 1), xaxt ="n", yaxt ="n",col =c("white", "black"))grid(L, L) }}
Optimization: Knapsack problem
1000 items, number 1, …, 1000
weight of item \(i^2\) is \(i^2\)
value of item \(i\) is \(i\log(i)\)
knapsack threshold is 1/3 of the total weight of all items
# test case# n = 4# w = c(1, 5, 3, 4)# v = c(15, 10, 9, 5)# L = 8# optimal solution is (1, 0, 1, 1) with value 29n =1000# number of itemsw = (1:n) ^2v = (1:n) *log(n)L =sum(w) /3
Run 100000 repetitions and print after each set of 10000 reps
plot(10001:Nsteps, wn[-(1:10000)], type ="l")abline(h = L, col ="orange")
Approximate optimal solution is highest value among states explored
# on which step did the max occuroptimal_step =which(vn ==max(vn))optimal_step
[1] 5369
# which items are includeditems_included =which(xn[max(optimal_step), ] ==1)# weight and value for items includedknapsack =data.frame(items_included,weight = w[items_included],value = v[items_included])DT::datatable(knapsack)
# total weightsum(knapsack$weight)
[1] 111267652
# as fraction of weight thresholdsum(knapsack$weight) / L
[1] 0.9999085
# total valuesum(knapsack$value)
[1] 1298092
# as fraction of total value of all itemssum(knapsack$value) /sum(v)
[1] 0.3754605
Cryptography
This code is very old and adapted from Dobrow, Introduction to Stochastic Processes with R.The code can be way simplified.
Text functions
# ascii(char) returns the numerical ascii value for charascii <-function(char) {strtoi(charToRaw(char), 16L) #get 'raw' ascii value} # charIndex takes in a character and returns its 'char value'# defined as a=1, b=2, ..., z=26, space=27# this matches the matrix in Austen countschar_index <-function(char) { aValue <-ascii(char)if (aValue ==32) {27 } else {#ascii sets "a" as 97, so subtract 96 aValue -96 }}
Function to decrypt the coded message given a particular cipher
# Decrypts code according to current_cipherdecrypt <-function(code, current_cipher) { out <- code# for each character in the message, decode it according to the current_cipherfor (i in1:nchar(code)) { char_ind <-char_index(substr(code, i, i))if (char_ind <27) {# change the ith character to the character determined by the current_ciphersubstr(out, i, i) <-rawToChar(as.raw(current_cipher[char_ind] +96)) } } out }
Secret message that we want to decode
# the scrambled messagecoded_message <-decrypt(message, sample(1:26))coded_message
[1] "adv cnqw qw x wcdiz xtt xmdec ndv oz tqpr sdc ptqffrg ceiarg efwqgr gdva xag qg tqhr cd cxhr x oqaecr bewc wqc iqsnc cnrir qtt crtt zde ndv q mryxor cnr fiqayr dp x cdva yxttrg mrtxqi qa vrwc fnqtxgrtfnqx mdia xag ixqwrg da cnr ftxzsideag vxw vnrir q wfrac odwc dp oz gxzw ynqttqa dec oxkqa irtxkqa xtt yddt xag xtt wnddcqas wdor mmxtt decwqgr dp cnr wynddt vnra x ydeftr dp sezw vnd vrir ef cd ad sddg wcxicrg oxhqas cidemtr qa oz arqsnmdinddg q sdc qa dar tqcctr pqsnc xag oz odo sdc wyxirg wnr wxqg zdeir odlqa vqcn zdei xeacqr xag eaytr qa mrtxqi"
Function to score the cipher based on how well the decoded message resembles real text
mat <-read.table("AustenCount.txt", header = F)logmat <-log(mat +1)# Computes the score of the decoded message using the given codescore <-function(code) { p <-0# For each pair of letters in the decoded message# query the transition matrix for the probability of that pairfor (i in1:(nchar(code) -1)){ p <- p + logmat[char_index(substr(code, i, i)),char_index(substr(code, i +1, i +1))] }# return the sum of these probabilities p}
Run 10000 steps of the MCMC algorithm
(Code records scores of ciphers we have already seen so we don’t have to rescore them if we revisit the cipher.)
# instantiate a map to hold previously computed codes' scoresmap <-new.env(hash = T, parent =emptyenv())# we begin with (a-> a, ..., z-> z) function for decrypting the coded messagecurrent_cipher <-1:27# calculate the score for current_cipher and store it in the mapcurrent_score <-score(decrypt(coded_message, current_cipher))map[[paste(current_cipher, collapse='')]] <- current_scoremax_score = current_scoremax_score_cipher = current_cipher# run 10000 iterations of the Metropolis-Hastings algorithmfor (iteration in1:10000) {# sample two letters to swap swaps <-sample(1:26, 2) old_cipher <- current_cipher# swap two letters current_cipher[swaps[1]] <- old_cipher[swaps[2]] current_cipher[swaps[2]] <- old_cipher[swaps[1]]# if we have already scored this decoding,# retrieve score from our mapif (exists(paste(current_cipher, collapse =''), map)) { new_score <- map[[paste(current_cipher, collapse ='')]] } else {# if we have not already scored this decoding,# calculate it and store it in the map new_score <-score(decrypt(coded_message, current_cipher)) map[[paste(current_cipher, collapse ='')]] <- new_score }if (new_score > max_score) { max_score = new_score max_score_cipher = current_cipher }# decide whether to accept current cipherif (runif(1) >exp(new_score - current_score)) { current_cipher <- old_cipher } else { current_score <- new_score }# print out our decryption every 1000 iterationsif ((iteration %%1000) ==0) {print(c(iteration, decrypt(coded_message, current_cipher))) }}
[1] "1000"
[2] "now tris is a stomy all acout row by like fot klipped tumned upside down and id lige to tage a binute qust sit mifrt treme ill tell you row i cehabe tre pminhe ok a town halled celaim in west priladelpria comn and maised on tre playfmound was wreme i spent bost ok by days hrillin out bavin melavin all hool and all srootinf sobe ccall outside ok tre shrool wren a houple ok fuys wro weme up to no food stamted baginf tmoucle in by neifrcomrood i fot in one little kifrt and by bob fot shamed sre said youme bozin witr youm auntie and unhle in celaim"
[1] "2000"
[2] "now tris is a stomy all acout row by life got flipped tumned upside down and id live to tave a binute just sit migrt treme ill tell you row i cehabe tre pminhe of a town halled celaim in west priladelpria comn and maised on tre playgmound was wreme i spent bost of by days hrillin out bakin melakin all hool and all srooting sobe ccall outside of tre shrool wren a houple of guys wro weme up to no good stamted baving tmoucle in by neigrcomrood i got in one little figrt and by bob got shamed sre said youme boxin witr youm auntie and unhle in celaim"
[1] "3000"
[2] "now this is a story all about how my life got flipped turned upside down and id live to tave a minute qust sit right there ill tell you how i became the prince of a town called belair in west philadelphia born and raised on the playground was where i spent most of my days chillin out maxin relaxin all cool and all shooting some bball outside of the school when a couple of guys who were up to no good started maving trouble in my neighborhood i got in one little fight and my mom got scared she said youre mokin with your auntie and uncle in belair"
[1] "4000"
[2] "now this is a story all about how my life got flipped turned upside down and id live to tave a minute qust sit right there ill tell you how i became the prince of a town called belair in west philadelphia born and raised on the playground was where i spent most of my days chillin out makin relakin all cool and all shooting some bball outside of the school when a couple of guys who were up to no good started maving trouble in my neighborhood i got in one little fight and my mom got scared she said youre moxin with your auntie and uncle in belair"
[1] "5000"
[2] "now this is a story all about how my life got flipped turned upside down and id live to tave a minute just sit right there ill tell you how i became the prince of a town called belair in west philadelphia born and raised on the playground was where i spent most of my days chillin out makin relakin all cool and all shooting some bball outside of the school when a couple of guys who were up to no good started maving trouble in my neighborhood i got in one little fight and my mom got scared she said youre moxin with your auntie and uncle in belair"
[1] "6000"
[2] "now this is a story all about how my life got flipped turned upside down and id live to tave a minute just sit right there ill tell you how i became the prince of a town called belair in west philadelphia born and raised on the playground was where i spent most of my days chillin out makin relakin all cool and all shooting some bball outside of the school when a couple of guys who were up to no good started maving trouble in my neighborhood i got in one little fight and my mom got scared she said youre moxin with your auntie and uncle in belair"
[1] "7000"
[2] "now this is a story all about how my life got flipped turned upside down and id live to tave a minute qust sit right there ill tell you how i became the prince of a town called belair in west philadelphia born and raised on the playground was where i spent most of my days chillin out makin relakin all cool and all shooting some bball outside of the school when a couple of guys who were up to no good started maving trouble in my neighborhood i got in one little fight and my mom got scared she said youre moxin with your auntie and uncle in belair"
[1] "8000"
[2] "now this is a story all about how my life got flipped turned upside down and id live to tave a minute just sit right there ill tell you how i became the prince of a town called belair in west philadelphia born and raised on the playground was where i spent most of my days chillin out makin relakin all cool and all shooting some bball outside of the school when a couple of guys who were up to no good started maving trouble in my neighborhood i got in one little fight and my mom got scared she said youre moxin with your auntie and uncle in belair"
[1] "9000"
[2] "now this is a story all about how my life got flipped turned upside down and id live to tave a minute just sit right there ill tell you how i became the prince of a town called belair in west philadelphia born and raised on the playground was where i spent most of my days chillin out makin relakin all cool and all shooting some bball outside of the school when a couple of guys who were up to no good started maving trouble in my neighborhood i got in one little fight and my mom got scared she said youre moxin with your auntie and uncle in belair"
[1] "10000"
[2] "now this is a story all about how my life got flipped turned upside down and id live to tave a minute qust sit right there ill tell you how i became the prince of a town called belair in west philadelphia born and raised on the playground was where i spent most of my days chillin out makin relakin all cool and all shooting some bball outside of the school when a couple of guys who were up to no good started maving trouble in my neighborhood i got in one little fight and my mom got scared she said youre moxin with your auntie and uncle in belair"
Message decoded using highest scoring cipher explored (should be pretty close to the real message; we can use our knowledge to fill in the rest).
decrypt(coded_message, max_score_cipher)
[1] "now this is a story all about how my life got flipped turned upside down and id live to tave a minute qust sit right there ill tell you how i became the prince of a town called belair in west philadelphia born and raised on the playground was where i spent most of my days chillin out makin relakin all cool and all shooting some bball outside of the school when a couple of guys who were up to no good started maving trouble in my neighborhood i got in one little fight and my mom got scared she said youre moxin with your auntie and uncle in belair"
