fr <- function(zzz){ # xx[1:18] <- zzz[1:18] xx[21] <- zzz[19] xx[19:20] <- 0.0 xx[22] <- 0.0 xx[23:24] <- zzz[20:21] sigmasq <- zzz[22] We have an additional parameter -- the variance term of the log-normal kappasq <- 100.0 The variance term of the Coordinate priors sumzsquare <- sum(zzz[1:21]**2) Sum of the squared coordinates # i <- 1 sse <- 0.0 while (i <= np) { ii <- 1 while (ii <= i) { j <- 1 dist_i_ii <- 0.0 while (j <= ndim) { # # Calculate distance between points # dist_i_ii <- dist_i_ii + (xx[ndim*(i-1)+j]-xx[ndim*(ii-1)+j])*(xx[ndim*(i-1)+j]-xx[ndim*(ii-1)+j]) j <- j + 1 } Note that this is the sum of squared differences of the logs if (i != ii)sse <- sse + (log(sqrt(TT[i,ii])) - log(sqrt(dist_i_ii)))**2 ii <- ii + 1 } i <- i + 1 } Here we have to include the variance term and the logs of the priors sse <- (((np*(np-1))/2.0)/2.0)*log(sigmasq)+(1/(2*sigmasq))*sse+(1/(2*kappasq))*sumzsquare return(sse) Note that above we have multiplied the log-likelihood by -1 so that the optimizers will find a minimum }The number of parameters is equal to 22 (nparam=22) now and further down in the problem before the optimizers are called zzz[22] is initialized to 0.2.
fr <- function(zzz){ # xx[1:8] <- zzz[1:8] There are 26 parameters: 14*2=28-3+1=26 xx[9:10] <- 0.0 The 14 points in 2 dimensions = 28 coordinates xx[11:19] <- zzz[9:17] The 5th point -- 490 Angstroms is the origin and xx[20] <- 0.0 the 2nd Dimension of 600 Angstroms is set to 0.0 xx[21:28] <- zzz[18:25] Hence, 25 coordinates plus the variance term = 26 parameters sigmasq <- zzz[26] kappasq <- 100.0 sumzsquare <- sum(zzz[1:25]**2) # i <- 1 sse <- 0.0 while (i <= np) { ii <- 1 while (ii <= i) { j <- 1 dist_i_ii <- 0.0 while (j <= ndim) { # # Calculate distance between points # dist_i_ii <- dist_i_ii + (xx[ndim*(i-1)+j]-xx[ndim*(ii-1)+j])*(xx[ndim*(i-1)+j]-xx[ndim*(ii-1)+j]) j <- j + 1 } if (i != ii)sse <- sse + (log(sqrt(TT[i,ii])) - log(sqrt(dist_i_ii)))**2 ii <- ii + 1 } i <- i + 1 } sse <- (((np*(np-1))/2.0)/2.0)*log(sigmasq)+(1/(2*sigmasq))*sse+(1/(2*kappasq))*sumzsquare return(sse) }The number of parameters is equal to 26 (nparam=26) and further down in the program before the optimizers are called zzz[26] is initialized to 0.2. This program also has much more efficient code for putting the coordinates estimated by the optimizers back into a matrix for plotting purposes:
# # Put Coordinates into matrix for plot purposes # xmetric <- rep(0,np*ndim) dim(xmetric) <- c(np,ndim) zdummy <- NULL zdummy[1:(np*ndim)] <- 0 Note the parentheses around "np*ndim" -- R does weird things if you write [1:np*ndim] zdummy[1:8] <- rhomax[1:8] zdummy[11:19] <- rhomax[9:17] I transfer the coordinates from the optimizer back into a vector np*ndim long with the points stacked zdummy[21:28] <- rhomax[18:25] # i <- 1 while (i <= np) { j <- 1 while (j <= ndim) { xmetric[i,j] <- zdummy[ndim*(i-1)+j] Now it is easy to put the coordinates back into an np by ndim matrix j <- j +1 } i <- i + 1 }
# InvCheck <- xhessian%*%XXInv checkinverse <- sum(abs(InvCheck-diag(nparam))) #and report checkinverse.
fr <- function(zzz){ # # A is set to the origin, (0.0, 0.0) and the 2nd dimension coordinate of # B is set to 0.0 Recall how dissimilar the letters "A" and "B" are # # (26+10)*2 = 72 # 72 - 3 + 1 = 70 There are 70 parameters: 36*2=72-3+1=70 # xx[1:2] <- 0.0 # A xx[3] <- zzz[1] xx[4] <- 0.0 # B 2nd dimension coordinate xx[5:72] <- zzz[2:69] sigmasq <- zzz[70] kappasq <- 100.0 sumzsquare <- sum(zzz[1:69]**2) etc. etc.
# InvCheck <- xhessian%*%XXInv checkinverse <- sum(abs(InvCheck-diag(nparam))) #and report checkinverse.