Series of homeworks from the course "Shape Analysis" teched by Prof. PaedDr. RNDr. Stanislav Katina, Ph.D. at Masaryk University - Science faculty.
%load_ext rmagic
%%R
seal = function(X){
return (cbind(rbind(X, X[1,])))
}
cplot = function(...){
# clean plot
plot(..., asp=1, axes=FALSE, xlab='', ylab='')
}
rotate = function(X, a){
T = matrix(c(
cos(a), sin(a),
-sin(a), cos(a)),
byrow=TRUE, nrow=2)
return(X %*% T)
}
# tetragon
X = rbind(c(-1,0), c(0,1), c(1,0), c(0,-1))
# plot tetragon
plot(seal(X), type='b', asp=1, pch=4, axes=FALSE, ylab='', xlab='')
lines(seal(rotate(X, 1)), type='b', asp=1, pch=5, lty=2)
lines(seal(rotate(X, -1)), type='b', asp=1, pch=6, lty=3)
legend('bottomleft', c('original', 'clockwise', 'counter-clockwise'), lty=1:3)
%%R
rotate = function(X, a){
# rotation matrix
T = matrix(c(
cos(a), sin(a),
-sin(a), cos(a)),
byrow=TRUE, nrow=2)
return(X %*% T)
}
translate = function(X, t){
# translation by given vector
X = cbind(X, rep(1, dim(X)[1]))
T = matrix(c(
1,0,t[1],
0,1,t[2],
0,0,1
), byrow=TRUE, nrow=3
)
T = t(T)
return((X %*% T)[,c(1,2)])
}
scale = function(X, k){
# scale by constant or vector
if(typeof(k) == "double"){
k = rep(k, dim(X)[1])
}
T = matrix(c(
k[1],0,
0,k[2]
), byrow=TRUE, nrow=2
)
return(X %*% T)
}
cs = function(X){
# centroid size
X = opt_translation(X)
return(sqrt(sum(diag(X %*% t(X)))))
}
centroid = function(X){
# centroids
return(colMeans(X))
}
opt_translation = function(X){
# optimal translation
return(translate(X, -centroid(X)))
}
opt_scaling = function(X, Y){
# optimal scaling from X to Y
k = rep(cs(Y) / cs(X), 2)
return(scale(X, k))
}
opt_rotation = function(X, Y){
# optimal rotation from X to Y
s = svd(t(X) %*% Y)
G = s$u %*% s$v
return(X %*% G)
}
%%R -w 800
fplot = function(X1, X2, sub=''){
plot(seal(X1), type='b', asp=1, axes=FALSE, ylab='', xlab='', pch=1:4, sub=sub)
lines(seal(X2), type='b', lty=2, asp=1, axes=FALSE, ylab='', xlab='', pch=1:4)
}
X1 = rbind(c(-1,0), c(0,1), c(1,0), c(0,-1))
X2 = rbind(c(-1,0), c(0,0.7), c(1,0), c(0,-0.9))
# transform X2
X2 = translate(scale(rotate(X2, 1), 0.5), c(0.2,0.2))
par(mfrow=c(2,2))
# original coordinates
fplot(X1, X2, sub='original coordinates')
# center coordinates
fplot(X1, opt_translation(X2), sub='centered coordinates')
# centered and scaled coordinates
fplot(X1, opt_scaling(opt_translation(X2), X1), sub='center and scaled coordinates')
# centered, scaled and rotated coordinates
fplot(X1, opt_rotation(opt_scaling(opt_translation(X2), X1), X1), sub='centered, scaled and rotated coordinates')
%%R
library(MASS)
apply3d = function(X, fun){
# apply function on matrices in third dimension
for(i in 1:dim(X)[3]){
X[,,i] = fun(X[,,i])
}
return(X)
}
procrust = function(M, threshold=0.00001){
# center and scale all matrices
M = apply3d(M, function(x) opt_translation(x))
M = apply3d(M, function(x) scale(x, 1./cs(x)))
# Procrust cycle
Xp = M[,,1]
while(TRUE){
last_Xp = Xp
# rotate to reference matrix
M = apply3d(M, function(x) opt_rotation(x, Xp))
# new reference matrix as a mean
Xp = apply(M, c(1,2), mean)
# break
cat('Procrust error:', sum(abs(last_Xp - Xp)), '\n')
if(sum(abs(last_Xp - Xp)) < threshold){
break
}
}
return(M)
}
gen_similar = function(Y, m=10){
# generate similar configuration matrices
M = array(0, dim=c(dim(Y)[1], dim(Y)[2], m))
for(k in 1:dim(Y)[1]){
M[k,,] = t(mvrnorm(n=m, mu=Y[k,], Sigma=diag(1,2) * 0.1^2))
}
return(M)
}
flatten_matrix = function(X){
return (cbind(as.vector(X[,1,]), as.vector(X[,2,])))
}
# original tetragons
X1 = rbind(c(-1,0), c(0,1), c(1,0), c(0,-1))
X2 = rbind(c(-1,0), c(0,2), c(1,0), c(0,-1))
M = gen_similar(X1, m=1000)
N = gen_similar(X2, m=1000)
%%R -w 800
fplot = function(M, ...){
plot(flatten_matrix(M), type = 'p', col = 'gray', pch=16, asp=1, axes=FALSE, ylab='', xlab='', ...)
# mean coordinates
X = apply(M, c(1,2), mean)
lines(seal(X), type='b', pch=1:4)
abline(v=0, col='red')
}
par(mfrow=c(1,2))
fplot(M, sub='original coordinates', xlim = c(-1.5, 1.5), ylim = c(-1.5, 2.5))
fplot(N, sub='original coordinates', xlim = c(-1.5, 1.5), ylim = c(-1.5, 2.5))
%%R -w 800
PM = procrust(M)
PN = procrust(N)
par(mfrow=c(1,2))
fplot(PM, sub='Procrustes coordinates', ylim = c(-0.6, 0.7))
fplot(PN, sub='Procrustes coordinates', ylim = c(-0.6, 0.7))
Procrust error: 0.1011805 Procrust error: 0.0002462531 Procrust error: 6.217558e-07 Procrust error: 0.1635371 Procrust error: 0.0002536229 Procrust error: 3.983855e-07
%%R -w 800
align_with_x = function(M, point1, point2){
# rotate shapes so that mean landmarks 2 and 4 lies on x-axis
v = rowMeans(M[point1,,]) - rowMeans(M[point2,,])
angle = atan(v[1] / v[2])
return(apply3d(M, function(x) rotate(x, angle)))
}
par(mfrow=c(1,2))
fplot(align_with_x(PM, 2, 4), sub='Rotated Procrustes coordinates', ylim = c(-0.6, 0.7))
fplot(align_with_x(PN, 2, 4), sub='Rotated Procrustes coordinates', ylim = c(-0.6, 0.7))
Let $SD(x)$ denote standard deviation in $x$ direction, $SD(y)$ sd in $y$ direction and $d$ distance from centroid. Then $k = \frac{\frac{1}{2} (SD(x) + SD(y))}{\sqrt{1 - d^2}}$ is supposed to be constant for all landmarks.
%%R
procrust_res = function(M){
sdx = apply(M[,1,],1,sd)
sdy = apply(M[,2,],1,sd)
cx = rowMeans(M[,1,])
cy = rowMeans(M[,2,])
c = apply(M, 2, mean)
d = sqrt((cx - c[1])^2 + (cy - c[2])^2)
sqrt_d = sqrt(1-d^2)
df = data.frame(d=d, sdx=sdx, sdy=sdy, '1-sqrt_d'=sqrt_d, k=0.5 * (sdx + sdy) / sqrt_d)
colnames(df) = c('d', 'SD(x)', 'SD(y)', '1-sqrt(d)', 'k')
return(df)
}
procrust_res(PN)
d SD(x) SD(y) 1-sqrt(d) k 1 0.3942879 0.02949443 0.03048001 0.9189870 0.03263073 2 0.6718438 0.02080161 0.02060487 0.7406929 0.02795118 3 0.3952984 0.02945890 0.02913935 0.9185527 0.03189705 4 0.4806373 0.02711922 0.02805762 0.8769195 0.03146060
%%R -w 1000
shear = function(X, ax=0, ay=0){
T = matrix(c(
1, tan(ay),
tan(ax), 1),
byrow=TRUE, nrow=2)
return(X %*% T)
}
# tetragon
X = rbind(c(0,0), c(0,1), c(1,1), c(1,0))
# shearing
par(mfrow=c(1,4))
plot(seal(X), type='b', asp=1, axes=FALSE, xlim=c(0,1.5),ylim=c(0,1.5), ylab='', xlab='', sub='original coordinates')
plot(seal(shear(X, ax=0.5)), type='b', asp=1, axes=FALSE, xlim=c(0,1.5), ylim=c(0,1.5), ylab='', xlab='', sub='horizontal shearing')
plot(seal(shear(X, ay=0.5)), type='b', asp=1, axes=FALSE, xlim=c(0,1.5), ylim=c(0,1.5), ylab='', xlab='', sub='vertical shearing')
plot(seal(shear(X, ax=0.5, ay=0.5)), type='b', asp=1, axes=FALSE, xlim=c(0,1.5), ylim=c(0,1.5), ylab='', xlab='', sub='horizontal and vertical shearing')
%%R
library(rgl)
shear3d_x = function(X, ax){
T = diag(3)
T[2,1] = tan(ax)
return(X %*% T)
}
shear3d_xy = function(X, ax, ay){
T = diag(3)
T[2,1] = tan(ax)
T[1,2] = tan(ay)
return(X %*% T)
}
shear3d_xyz = function(X, ax, ay, az){
T = diag(3)
T[3,1] = tan(ax)
T[3,2] = tan(ay)
T[2,3] = tan(az)
return(X %*% T)
}
save_3d_image = function(X, p, filename){
open3d()
rgl.viewpoint(theta=0,phi=0,fov=30)
bg3d('white')
par3d(windowRect=c(100,100,600,600))
# draw
spheres3d(X, radius=0.1, col='black', add=TRUE)
segments3d(X[p,], col='black', lwd=5, add=TRUE)
# save as a file
snapshot3d(filename)
}
# cube
X1 = c(-1,-1,-1)
X2 = c(-1,-1,1)
X3 = c(-1,1,-1)
X4 = c(-1,1,1)
X5 = c(1,-1,-1)
X6 = c(1,-1,1)
X7 = c(1,1,-1)
X8 = c(1,1,1)
X = rbind(X1,X2,X3,X4,X5,X6,X7,X8)
# connect vertices
p = c()
for(i in 1:7){
for(j in (i+1):8){
if(sum(abs(X[i,] - X[j,])) == 2){
p = c(p, c(i,j))
}
}
}
save_3d_image(X, p, 'original.png')
save_3d_image(shear3d_x(X, 0.5), p, 'x.png')
save_3d_image(shear3d_xy(X, 0.5, 0.5), p, 'xy.png')
save_3d_image(shear3d_xyz(X, 0.5, 0.5, 0.5), p, 'xyz.png')
from IPython.display import Image
Image(filename='original.png', width=300, height=300)
Image(filename='x.png', width=300, height=300)
Image(filename='xy.png', width=300, height=300)
Image(filename='xyz.png', width=300, height=300)
%%R -w 1000
reflection = function(X, ax=1){
if(ax==2){
T = matrix(c(
1, 0,
0, -1),
byrow=TRUE, nrow=2)
} else {
T = matrix(c(
-1, 0,
0, 1),
byrow=TRUE, nrow=2)
}
return(X %*% T)
}
X = rbind(c(-1,0), c(-1,1), c(1,1.2), c(1,0))
par(mfrow=c(1,2))
plot(seal(X), type='b', asp=1, axes=FALSE, ylab='', xlab='', sub='original coordinates', pch=1:4)
abline(v=0, col='gray')
plot(seal(reflection(X, ax=1)), type='b', asp=1, axes=FALSE, ylab='', xlab='', sub='reflected coordinates', pch=1:4)
abline(v=0, col='gray')
%%R -w 800
Q = matrix(0,4,4)
Q[1,4] = 1
Q[2,3] = 1
Q[3,2] = 1
Q[4,1] = 1
X = rbind(c(-1,0), c(-1,1), c(1,1.2), c(1,0))
XR = reflection(X, ax=1)
print(Q)
par(mfrow=c(1,2))
plot(seal(X), type='b', asp=1, axes=FALSE, ylab='', xlab='', sub='original coordinates', pch=1:4)
abline(v=0, col='gray')
plot(seal(Q %*% XR), type='b', asp=1, axes=FALSE, ylab='', xlab='', sub='reflected and relabeled coordinates', pch=1:4)
abline(v=0, col='gray')
[,1] [,2] [,3] [,4] [1,] 0 0 0 1 [2,] 0 0 1 0 [3,] 0 1 0 0 [4,] 1 0 0 0
After writing out model IM1 we get equivalent system of equations \begin{align} y &\sim Sw + 1_k c + xa \\ 0 &\sim 1_k w \\ 0 &\sim x^T a \end{align} model IM1a can be breaken down into the same system. The difference between IM1 a IM1a is only in interchange of rows/columns and they are equivalent to each other.
See first question. Matrix $X_p$ looks like that to fulfill conditions $0 = 1_k w$ a $0 = x^T a$.
Assuming full rank we get
\begin{align} \hat{\beta} &= (X_p^T X_p)^{-1} X_p^T y_p \\ \hat{\beta^*} &= (X_p^{*T} X_p^*)^{-1} X_p^{*T} y_p \end{align}model = lm(yp ~ Xp)
beta = summary(model)$coefficients
where $H$ is called hat matrix or projective matrix
%%R -w 800
S_matrix = function(x,y){
# distane matrix
k = length(x)
S = 1/12 * abs(x %*% t(rep(1,k)) - rep(1,k) %*% t(x))^3
}
IM1 = function(x, y){
k = length(x)
# distance matrix
S = S_matrix(x,y)
# design matrix
Xp = matrix(0,k+2,k+2)
Xp[3:(k+2),3:(k+2)] = S
Xp[k+1,1:k] = rep(1,k)
Xp[k+2,1:k] = x
Xp[1:k,k+1] = rep(1,k)
Xp[1:k,k+2] = x
# yp
yp = c(y, c(0,0))
# linear model
model = lm(yp ~ Xp - 1)
# estiamtion of y
yh = predict(model)[1:k]
return(yh)
}
# sinus function
n = 101
x = seq(0,1,length=n)
y = sin(2*pi*x^3)^3 + 0.1 * rnorm(n)
# IM1 model
par(mfrow=c(1,2))
yh = IM1(x, y)
plot(x, y, sub='IM1 model')
lines(x, yh, col='blue')
# spline method
sp = spline(x,y,method='natural')
plot(x, y, sub='spline(x,y,method=natural)')
lines(sp$x, sp$y, col='red')
%%R
library(rgl)
save_3d = function(S, x, filename){
slim = range(S)
ncols = 50
colorlut = terrain.colors(ncols + 1) # height color lookup table
foo = trunc(ncols * (S - slim[1]) / (slim[2] - slim[1])) + 1
col = colorlut[foo] # assign colors to heights for each point
open3d()
rgl.viewpoint(theta=0,phi=30,fov=60)
par3d(windowRect=c(100,100,600,600))
rgl.surface(x,x,S, color=col, back='lines')
bg3d("white")
decorate3d()
snapshot3d(filename)
}
S = 10 * S_matrix(x,y)
save_3d(S, x, 'basal.png')
from IPython.display import Image
Image(filename='basal.png')
%%R
S = S_matrix(x,y)
image(x, x, S, asp=1, col=terrain.colors(30), axes=FALSE, xlab='', ylab='')
contour(x,x,S,add=TRUE)
abline(a=0,b=1, col='red', lw=2)
%%R -w 800 -h 800
sqrt_matrix = function(a){
# hack positive definiteness by replacing negative values with zeros
a.eig <- eigen(a)
eig_values = pmax(a.eig$values, 0)
a.sqrt <- a.eig$vectors %*% diag(sqrt(eig_values)) %*% solve(a.eig$vectors)
return(a.sqrt)
}
PRM1 = function(x, y, lam, method=1, return_model=FALSE){
# method == 0 -> estimation by hat matrix
# method == 1 -> estimation with lm function
k = length(x)
# matrix Sp
S = S_matrix(x,x)
Sp = matrix(0,k+2,k+2)
Sp[3:(k+2),3:(k+2)] = S
# matrix X
R = sqrt_matrix(Sp)
# matrix Xdm
Xdm = cbind(rep(1,k), x, S)
if(method == 0){
# matrix H
H = Xdm %*% solve(t(Xdm) %*% Xdm + lam * Sp) %*% t(Xdm)
# estimate y
yh = H %*% y
return(yh)
} else if(method == 1){
# matrix Xp
Xp = rbind(Xdm, sqrt(lam) * R)
# vector yp
yp = c(y, rep(0, k+2))
# linear model
model = lm(yp ~ Xp - 1)
if(return_model){
return(model)
}else{
yh = predict(model)[1:k]
return(yh)
}
}
}
polyfit = function(x, y, d=12){
return(predict(lm(y ~ poly(x, d, raw=TRUE))))
}
lam = 4.774251e-6
par(mfrow=c(2,2))
plot(x, y, sub='smooth.spline')
sp = smooth.spline(x,y,all.knots = TRUE, cv=TRUE)
lines(sp$x, sp$y, col='blue')
plot(x, y, sub='PRM1')
lines(x, PRM1(x,y,lam=lam,method=0), col='blue')
plot(x, y, sub='12 degree polynomial')
lines(x, polyfit(x,y,d=12), col='blue')
plot(x, y, sub='PRM1a')
lines(x, PRM1(x,y,lam=lam,method=1), col='blue')
%%R -w 1000 -h 700
library(psych)
MSS = function(x,y,lam,fun){
# return MSS error for estimation PRM1 (assuming you know the original function)
yh = PRM1(x,y,lam)
# real function values
yr = fun(x)
return(mean((yh - yr)^2))
}
OCV = function(x, y, lam){
k = length(x)
model = PRM1(x, y, lam, return_model=TRUE)
hat = influence(model)$hat[1:k]
yh = predict(model)[1:k]
return(mean((yh - y)^2 / (1 - hat)^2))
}
GCV = function(x, y, lam){
k = length(x)
model = PRM1(x, y, lam, return_model=TRUE)
hat = influence(model)$hat[1:k]
yh = predict(model)[1:k]
return(k*(sum((yh - y)^2) / tr(diag(k) - hat)^2))
}
plot_cv = function(title, x, y, objfun, yfun,...){
fun = function(lam) objfun(x,y,lam,...)
# find minimum lambda
f = function(lam) objfun(x,y,exp(lam), ...)
opt = optimize(f=f, lower=-20, upper=0)
opt_lam = exp(opt$minimum)
# plot surrounding lambdas
lambdas = seq(opt_lam /5, opt_lam * 5, length.out=30)
error = lapply(lambdas, fun)
# plot lambdas
plot(lambdas, error, type='l', ylab=paste(title, '(lambda)'), xlab='lambda',
sub=paste('lambda(MSS)=', sprintf("%.7f",opt_lam)))
points(opt_lam, fun(opt_lam), pch=19, lwd=5)
# plot spline
plot(x, y, sub=paste('penalized regression spline', title))
lines(x, PRM1(x,y,opt_lam,method=1), col='blue')
lines(x, yfun(x), col='red')
legend('topleft', c('E[y]', 'odhad y'), col=c('red', 'blue'), lty=1)
}
# MSS
fun = function(x) sin(2*pi*x^3)^3
par(mfcol=c(2,3))
plot_cv('MSS', x, y, MSS, fun, fun)
plot_cv('OCV', x, y, OCV, fun)
plot_cv('GCV', x, y, GCV, fun)
%%R -w 800
cumchord = function(X){
# cumulative chord distance
return(cumsum(sqrt(apply((X-rbind(X[1,],X[-(nrow(X)),]))^2,1,sum))))
}
resample = function(X, method, eps=1e-8){
x = X[,1]
y = X[,2]
i = 0
while(TRUE){
i = i + 1
# cumulative chordal distance
CH <- cumchord(cbind(y,x))
# if distance is small enough, stop
dCH = diff(CH)
error = sum(abs(dCH - mean(dCH)))
if(error < eps) break
# new splines
if(method == 'spline'){
x = spline(CH, x, method="natural", n=50)$y
y = spline(CH, y, method="natural", n=50)$y
} else if(method == 'approx'){
x = spline(CH, x, n=50)$y
y = spline(CH, y, n=50)$y
}
}
cat(i, 'iterations\n')
return(cbind(x,y))
}
interpolation = function(X, n=50){
CH <- cumchord(cbind(X[,2],X[,1]))
x = spline(CH, X[,1], n=n)$y
y = spline(CH, X[,2], n=n)$y
return(cbind(x,y))
}
# load data from file
X = read.table("symphysis.txt",header=TRUE)
# draw resampled curve
par(mfrow=c(1,3))
plot(interpolation(X), type='l', asp=1, axes=FALSE, xlab='', ylab='', sub='interpolation')
points(X[,1], X[,2])
plot(resample(X, method='spline'), type='o', asp=1, axes=FALSE, xlab='', ylab='', sub='spline', col='blue')
plot(resample(X, method='approx'), type='o', asp=1, axes=FALSE, xlab='', ylab='', sub='approx', col='red')
7 iterations 7 iterations
%%R
# phi distance function
phi = function(x){
d = sqrt(sum(x^2))
if(d == 0){
return(0)
}else{
return(d^2 * log(d))
}
}
S_matrix3 = function(X,Y){
# return distance matrix between rows of X a Y
m = dim(X)[1]
n = dim(Y)[1]
S = matrix(, nrow = m, ncol = n)
for(i in 1:m){
for(j in 1:n){
S[i,j] = phi(X[i,] - Y[j,])
}
}
return(S)
}
S = S_matrix3(X,X)
save_3d(S, seq(0,2000,length=dim(S)[1]), 'basal3d.png')
from IPython.display import Image
Image(filename='basal3d.png')
%%R
S = S_matrix3(X,X)
x = seq(0,2000,length=dim(S)[1])
image(x, x, S, asp=1, col=terrain.colors(30), axes=FALSE, xlab='', ylab='')
contour(x,x,S,add=TRUE)
abline(a=0,b=1, col='red', lw=2)
%%R -w 1000
net = function(X, kx=20, margin=0.2){
x1 = min(X[,1])
x2 = max(X[,1])
y1 = min(X[,2])
y2 = max(X[,2])
# add 1/5 margin to all directions
x_range = x2 - x1
y_range = y2 - y1
x1 = x1 - margin * x_range
x2 = x2 + margin * x_range
y1 = y1 - margin * y_range
y2 = y2 + margin * y_range
x_grid = seq(x1,x2,length=kx)
xmin = min(x_grid)
xmax = max(x_grid)
hy = (xmax - xmin) / (kx - 1)
y_grid = seq(y1,y2,by=hy)
ky = length(y_grid)
Xcp <- as.matrix(expand.grid(x_grid,y_grid))
return(Xcp)
}
draw_net = function(X, Xcp, kx=20, net_type='lines', type='b', ...){
# show the shape
X = cbind(rbind(X, X[1,]))
plot(X, asp=1, type=type, xlim=range(Xcp[,1]), ylim=range(Xcp[,2]), axes=FALSE, xlab='', ylab='', pch=19, ...)
# show net
if(net_type == 'lines'){
ky = dim(Xcp)[1] / kx
for(i in 1:ky){
lines(Xcp[1:kx + (i-1)*kx,], lwd=0.3)
}
for(i in 1:kx){
lines(Xcp[(1:ky)*kx-i+1,], lwd=0.3)
}
} else if(net_type == 'points'){
points(Xcp, pch=19, cex=0.2)
}
}
# shapes
X1 = rbind(c(-0.5,-0.5),c(-0.5,0.5),c(0.5,0.5),c(0.5,-0.5))
X2 = rbind(c(-0.5,-0.65),c(-0.5,0.65),c(0.5,0.65),c(0.5,-0.65))
# TPS net
par(mfrow=c(1,4))
draw_net(X1, net(X1, kx=20), kx=20, net_type='points')
draw_net(X2, net(X2, kx=20), kx=20, net_type='points')
draw_net(X1, net(X1, kx=20), kx=20, net_type='lines')
draw_net(X2, net(X2, kx=20), kx=20, net_type='lines')
%%R -w 800
shear = function(X, ax=0, ay=0){
T = matrix(c(
1, tan(ay),
tan(ax), 1),
byrow=TRUE, nrow=2)
return(X %*% T)
}
X = rbind(c(-0.5,-0.5), c(-0.5,0.5), c(0.5,0.5), c(0.5, -0.5))
XN = net(X)
par(mfrow=c(1,2))
draw_net(shear(X, ax=0.5), shear(XN, ax=0.5))
draw_net(shear(X, ay=0.5), shear(XN, ay=0.5))
%%R -w 1000
extrapolation = function(X,Y,l){
return(Y + l * (Y - X))
}
IM3_Xp = function(X,Y){
k = dim(X)[1]
S = S_matrix3(X,X)
Xp = rbind(cbind(S, rep(1,k), X),
c(rep(1,k), rep(0,3)),
cbind(t(X), matrix(0,2,3)))
return(Xp)
}
IM3 = function(X,Y,Xcp){
k = dim(X)[1]
if(is.vector(Y)){
# IM2 model
Yp = c(Y,rep(0,3))
} else {
# IM3 model
Yp = rbind(Y,matrix(0,3,2))
}
Xp = IM3_Xp(X,Y)
# linear model a its coefficients
model = lm(Yp ~ Xp - 1)
beta = coef(model)
# new transformed matrix Xcp
Scp = S_matrix3(Xcp,X)
Xpcp = rbind(cbind(Scp, rep(1,k), Xcp),
c(rep(1,k), rep(0,3)),
cbind(t(X), matrix(0,2,3)))
Ypcp = Xpcp %*% beta
return(Ypcp[1:dim(Xcp)[1],])
}
plot_extrapolation = function(X,Y,l, only_net=FALSE, kx=20, ...){
Z = extrapolation(X,Y,l)
# shapes
if(!only_net){
cplot(seal(X), type='b', lty=2, sub=paste('l =', l))
lines(seal(Z), type='b')
# arrows
for(i in 1:dim(Z)[1]){
arrows(X[i,1], X[i,2], Z[i,1], Z[i,2], code=2, length=0.1, col='blue')
}
}
# net
Zcp = IM3(X, Z, net(X, kx=kx))
draw_net(Z, Zcp, sub=paste('l =', l), kx=kx, ...)
}
X1 = rbind(c(0,-0.5), c(-0.5,0), c(0,0.5), c(0.5,0))
X2 = rbind(c(0,-0.51), c(-0.51,0), c(0,0.4), c(0.51,0))
par(mfrow=c(1,2))
plot_extrapolation(X1, X2, 0)
plot_extrapolation(X1, X2, 1)
%%R -w 800
par(mfrow=c(1,3))
plot_extrapolation(X1, X2, 1, only_net=TRUE)
plot_extrapolation(X1, X2, 4, only_net=TRUE)
plot_extrapolation(X1, X2, 8, only_net=TRUE)
%%R -w 800
library(psych)
bending_energy = function(X,Y){
k = dim(X)[1]
L = IM3_Xp(X,Y)
Be = solve(L)[1:k,1:k]
return(tr(t(Y) %*% Be %*% Y))
}
X = rbind(c(-0.5,-0.5), c(-0.5,0.5), c(0.5,0.5), c(0.5,-0.5))
X1 = rbind(X, c(0,0))
X2 = rbind(X, c(0.3,0.3))
Xcp = net(X1)
be = bending_energy(X1,X2)
par(mfrow=c(1,2))
draw_net(X1, Xcp, type='p', cex=2, sub=paste('Bending energy', 0))
draw_net(X2, IM3(X1,X2,Xcp), type='p', cex=2, sub=paste('Bending energy', be))
%%R -w 1000
par(mfrow=c(1,3))
# sipka
cplot(X1, pch=19, cex=1.5, sub='quint')
arrows(-0.5, -0.5, 0.5, 0.5, col='red', length=0.2, lwd=2)
# bending energy 1d
x = seq(-0.5,0.5,length=100)
y = sapply(x, function(e) bending_energy(X1, rbind(X, c(e,e))))
plot(x,y,type='l',xlab='quint diagonal', ylab='bending energy')
# bending energy 2d
B = matrix(0,100,100)
for(i in 1:100){
for(j in 1:100){
B[i,j] = bending_energy(X1, rbind(X, c(x[i], x[j])))
}
}
image(x, x, B, asp=1, col=terrain.colors(30), axes=FALSE, xlab='', ylab='', sub='bending energies')
contour(x,x,B,add=TRUE)
%%R -w 1000 -h 800
# IM3 model is programmed in a way to work with all dimensions of y
IM2 = IM3
PRM2 = function(X, y, Xcp=NULL, lam=NULL, return_model=FALSE){
k = dim(X)[1]
# matrix
yp = c(y, rep(0, k+3))
# matrix Sp
S = S_matrix3(X,X)
Sp = matrix(0,k+3,k+3)
Sp[4:(k+3),4:(k+3)] = S
# matrix X
R = sqrt_matrix(Sp)
# matrix Xdm
Xdm = cbind(rep(1,k), X, S)
# matrix Xp
Xp = rbind(Xdm, sqrt(lam) * R)
# linear model a its coefficients
model = lm(yp ~ Xp - 1)
if(return_model){
return(model)
}
beta = coef(model)
# new transformed matrix Xp
l = dim(Xcp)[1]
Scp = S_matrix3(Xcp,X)
Xdm_cp = cbind(rep(1,l), Xcp, Scp)
ypcp = Xdm_cp %*% beta
return(ypcp[1:l])
}
GCV_PRM2 = function(X, y, lam){
k = dim(X)[1]
model = PRM2(X, y, lam=lam, return_model=TRUE)
hat = influence(model)$hat[1:k]
yh = predict(model)[1:k]
return(k*(sum((yh - y)^2) / tr(diag(k) - hat)^2))
}
opt_lambda = function(X, y){
# find minimum lambda for PRM2
f = function(lam) GCV_PRM2(X,y,exp(lam))
opt = optimize(f=f, lower=-20, upper=0)
opt_lam = exp(opt$minimum)
return(opt_lam)
}
implode.grid = function(Xcp, y){
# long vectors useful for image format
xx = unique(Xcp[,1])
yy = unique(Xcp[,2])
z = matrix(y, length(xx), length(yy))
return(list(x=xx,y=yy,z=z))
}
fplot = function(X,Xcp,ycp,...){
r = implode.grid(Xcp, ycp)
image(r$x, r$y, r$z, xlab='x1.coords', ylab='x2.coords',...)
contour(r$x, r$y, r$z, add=TRUE)
lines(X, type='p', pch=20, cex=1.5)
axis(1, at=unique(X[,1]))
axis(2, at=unique(X[,2]))
}
# load data
t = read.table('VCG01.txt')
X = cbind(t[['x1.coords']], t[['x2.coords']])
y = t[['y1']]
par(mfrow=c(2,1))
# IM2 model
Xcp = net(X, kx=100, margin=0.03)
ycp = IM2(X, y, Xcp)
fplot(X, Xcp, ycp, sub='IM2')
# PRM2 model
opt_lam = opt_lambda(X,y)
cat('Optimal lambda =', opt_lam)
ycp = PRM2(X, y, Xcp, lam=opt_lam)
fplot(X, Xcp, ycp, sub=paste('PRM2 (lambda =', opt_lam, ')'))
Optimal lambda = 3.177167e-07
%%R -w 1000
draw_circle_net = function(X, Xcp, ycp, r=0.9, ...){
# assign NaN values outside circle
g = implode.grid(Xcp, ycp)
Z = g$z
for(i in 1:length(g$x)){
for(j in 1:length(g$y)){
if(g$x[i]^2 + g$y[j]^2 > r^2){
Z[i,j] = NA
}
}
}
image(g$x, g$y, Z, asp=1, xlab='', ylab='', axes=FALSE, ...)
contour(g$x, g$y, Z, add=TRUE)
# show original points
points(X, pch=19, cex=1.5)
# show net
d = Xcp[,1]^2 + Xcp[,2]^2 <= r^2
points(Xcp[d,], pch=19, cex=0.2)
# circle
x = r * cos(seq(0,2*pi,0.01))
y = r * sin(seq(0,2*pi,0.01))
lines(x,y)
}
# load data
t = read.table('EEG12.txt')
X = cbind(t[['x1.coords']], t[['x2.coords']])
y = t[['y']]
par(mfrow=c(1,2))
Xcp = net(X, kx=100, margin=0.1)
# IM2 model
ycp = IM2(X, y, Xcp)
draw_circle_net(X, Xcp, ycp, r=0.9, sub='IM2')
# PRM2 model
opt_lam = opt_lambda(X,y)
cat('Optimal lambda =', opt_lam)
ycp = PRM2(X, y, Xcp, lam=opt_lam)
draw_circle_net(X, Xcp, ycp, r=0.9, sub=paste('PRM2 (lambda =', opt_lam, ')'))
Optimal lambda = 7.485854e-09
%%R -w 800
g = read.table('gorilas.txt', header=TRUE)
sex = g[,1]
g = g[,2:length(g)]
g = as.matrix(g)
M = array(0, c(8,2,59))
for(i in 1:dim(g)[1]){
x = g[i,][1:8]
y = g[i,][9:16]
M[,1,i] = x
M[,2,i] = y
}
# permute data for easier plotting
M = M[c(2,3,4,5,1,6,7,8),,]
fplot = function(M, ...){
for(i in 1:dim(M)[3]){
if(i==1){
plot(seal(M[,,i]), type='b', col='gray', pch=19, asp=1, axes=FALSE, ylab='', xlab='',
xlim=range(as.vector(M[,1,])), ylim=range(as.vector(M[,2,])), ...)
} else {
lines(seal(M[,,i]), type='b', col='gray', pch=19)
}
points(M[,,i], pch=19)
}
}
align_with_x = function(M, point1, point2){
# rotate shapes so that mean landmarks lie on x-axis
v = rowMeans(M[point1,,]) - rowMeans(M[point2,,])
angle = atan2(v[2], v[1])
return(apply3d(M, function(x) rotate(x, -angle)))
}
PM = procrust(M)
par(mfrow=c(1,2))
fplot(PM, sub='Procrust')
PM = align_with_x(PM, point1=1, point2=5)
fplot(PM, sub='Procrust a rotace')
Procrust error: 0.1448269 Procrust error: 3.802056e-05 Procrust error: 3.162222e-08
%%R -w 1000
par(mfrow=c(2,2))
fplot(PM)
# mean shapes
XF = apply(PM[,,sex=='f'], c(1,2),mean)
XM = apply(PM[,,sex=='m'], c(1,2),mean)
cplot(seal(XF), col='red', type='b')
lines(seal(XM), col='blue', type='b')
# extrapolation
plot_extrapolation(XM, XF, 3, only_net=TRUE, col='red', kx=40)
plot_extrapolation(XF, XM, 3, only_net=TRUE, col='blue', kx=40)
%%R -w 1000 -h 400
# helper function to substract means
center_mean = function(x, mu=NULL) {
if(is.null(mu)){
mu = colMeans(x)
}
x_mean = rep(1, nrow(x)) %*% t(mu)
return(x - x_mean)
}
# confidence ellipsoid
confidence_ellipsoid = function(d, alpha=0.05){
n = dim(d)[1]
k = dim(d)[2]
# svd
S = cov(d)
eig = eigen(S)
# ellipsis length
ds = sqrt(eig$values * (n-1) * k / (n-k) * qf(1-alpha, k, n-k))
# non-rotated ellipsis
a = seq(0, 2*pi, len=100)
x = ds[1] * cos(a)
y = ds[2] * sin(a)
M = cbind(x,y)
# rotate it to direction of eigenvectors
M = M %*% t(eig$vectors)
# add means
mu = colMeans(d)
M = center_mean(M, -mu)
return(M)
}
# show two componenta
sc = function(T, i, j, evals, category){
T = T[,c(i,j)]
xlab = paste('PC', i, ' (', round(100*evals[i], 2), "%", ')', sep="")
ylab = paste('PC', j, ' (', round(100*evals[j], 2), "%", ')', sep="")
color = rep(NA, length(category))
color[category == 'm'] = 'blue'
color[category == 'f'] = 'red'
lim = max(abs(T)) * 1.2
plot(T, asp=1, xlab=xlab, ylab=ylab, col=color, pch=19, xlim=c(-lim,lim), ylim=c(-lim,lim))
m_male = apply(T[category == 'm',], 2, mean)
m_female = apply(T[category == 'f',], 2, mean)
points(rbind(m_male, m_female), pch=10, cex=3)
legend('bottomright', legend=c('male', 'female'), col = c('blue', 'red'), pch=19)
# highlight zero axes
abline(h=0, col='grey')
abline(v=0, col='grey')
# confidence ellipsoid
lines(confidence_ellipsoid(T[category == 'm',]), col='blue')
lines(confidence_ellipsoid(T[category == 'f',]), col='red')
}
pca_scores = function(D, category){
# center matrix
m = apply(D, 2, mean)
D = D - as.matrix(rep(1,dim(D)[1])) %*% t(as.matrix(m))
# covariance matrix
Sx = cov(D)
# eigenvalues
f = eigen(Sx)
evals = f$values
evec = f$vectors
# explained variance
evals = evals / sum(evals)
# skore
T = D %*% evec
# show components
par(mfrow=c(1,3))
sc(T, 1, 2, evals, category)
sc(T, 1, 3, evals, category)
sc(T, 2, 3, evals, category)
}
flatten = function(PM){
D = rbind(PM[,1,], PM[,2,])
D = t(D)
return(D)
}
D = flatten(PM)
pca_scores(D, sex)
%%R -w 800
pca_reduction = function(D, components=c(1)){
# centered matrix
mu = apply(D, 2, mean)
D = D - as.matrix(rep(1,dim(D)[1])) %*% t(as.matrix(mu))
# covariance matrix
Sx = cov(D)
# eigenvalues
f = eigen(Sx, symmetric=TRUE)
evals = f$values
evec = f$vectors
# PCA score
base = evec[,components]
Y = D %*% base
# back to original basis
X = Y %*% t(base) + rep(1, nrow(Y)) %*% t(mu)
return(X)
}
X = pca_reduction(D, components=c(1))
# transform back to 3dim
M = array(0, c(8,2,59))
for(i in 1:dim(X)[1]){
M[,1,i] = X[i,1:8]
M[,2,i] = X[i,9:16]
}
# extrapolation of original coordinates
par(mfrow=c(2,2))
plot_extrapolation(XM, XF, 3, only_net=TRUE, col='red', kx=40)
plot_extrapolation(XF, XM, 3, only_net=TRUE, col='blue', kx=40)
# extrapolation of new coordinates
MF = apply(M[,,sex=='f'], c(1,2),mean)
MM = apply(M[,,sex=='m'], c(1,2),mean)
plot_extrapolation(MM, MF, 3, only_net=TRUE, col='red', kx=40)
plot_extrapolation(MF, MM, 3, only_net=TRUE, col='blue', kx=40)
%%R
hotelling_test = function(X,Y,comp){
D = rbind(X,Y)
# matrix centering
mu = apply(D, 2, mean)
D = D - as.matrix(rep(1,dim(D)[1])) %*% t(as.matrix(mu))
# covariance matrix
Sx = cov(D)
# eigenvalues
f = eigen(Sx, symmetric=TRUE)
evals = f$values
evec = f$vectors
# score for mean difference
r = (colMeans(X) - colMeans(Y)) %*% evec
# statistic
n1 = dim(X)[1]
n2 = dim(Y)[1]
s = min(dim(X)[2], n1+n2-2)
Th2 = sum((r^2 / evals)[1:s])
cat('Contribution of principal components into statistic:\n', r^2 / evals, '\n')
Fh = n1*n2*(n1+n2-s-1) / ((n1+n2)*(n1+n2-2)*s) * Th2
p_val = pf(Fh, s, n1+n2-s-1)
return(1-p_val)
}
p_val = hotelling_test(flatten(PM[,,sex=='m']), flatten(PM[,,sex=='f']))
cat('p-value of Hotteling test ', p_val)
Contribution of principal components into statistic: 2.111833 0.5542956 0.0006460684 0.07493997 0.5802357 0.0276503 0.02825453 0.00345909 0.02343099 7.929654e-05 0.01540087 0.01534201 0.004033227 1.306697e-16 -2.746089e-15 -8.845693e-15 p-value of Hotteling test 0.01420651
because p-value is $0.014 < 0.05$, we reject hypothesis of identical gender dimorphism. The greatest influence on testing statistic based on $T_H$ has first component followed by the fourth and the second.