Creare un grafico in R con la retta di regressione e le distanze dei punti


x <- c(173, 169, 176, 166, 161, 164, 160, 158, 180, 187)
y <- c(80, 68, 72, 75, 70, 65, 62, 60, 85, 92)

# plottiamo i punti e la retta di regressione
mod1 <- lm(y ~ x)
plot(x, y, xlim=c(min(x)-5, max(x)+5), ylim=c(min(y)-10, max(y)+10))
abline(mod1, lwd=2)



# calcoliamo (e approssimiamo) i residui e i valori predicted
res <- signif(residuals(mod1), 5)
pre <- predict(mod1)

# tracciamo i segmenti di distanza tra i punti e la retta di regressione
segments(x, y, x, pre, col="red")

# aggiungiamo i valori dei residui come labels dei punti
library(calibrate)
textxy(x, y, res, cx=0.7)


R is a cool image editor!

Ecco qualche funzione per ricreare alcuni dei più famosi effetti da applicare sulle immagini, disponibili nella gran parte degli editor (Photoshop incluso)
Per poter caricare e modificare le immagini, è necessario installare la library rimage.
Per caricare una immagina, il codice è il seguente:

y <- read.jpeg("path")

Per plottare l'immagine, usiamo invece:

plot(y)

---

Original image

---

Sepia tone

rgb2sepia <- function(img){
 iRed <- img[,,1]*255
 iGreen <- img[,,2]*255
 iBlue <- img[,,3]*255
 
 oRed <- iRed * .393 + iGreen * .769 + iBlue * .189
 oGreen <- iRed * .349 + iGreen * .686 + iBlue * .168
 oBlue <- iRed * .272 + iGreen * .534 + iBlue * .131
 
 qw <- array( c(oRed/255 , oGreen/255 , oBlue/255), dim=c(dim(iRed)[1],dim(iRed)[2],3) )
 
 imagematrix(qw, type="rgb")
}
 
plot(rgb2sepia(y))

---

Negative

rgb2neg <- function(img){
 iRed <- img[,,1]
 iGreen <- img[,,2]
 iBlue <- img[,,3]
 
 oRed <- (1 - iRed)
 oGreen <- (1 - iGreen)
 oBlue <- (1 - iBlue)
 
 qw <- array( c(oRed, oGreen, oBlue), dim=c(dim(iRed)[1],dim(iRed)[2],3) )
 
 imagematrix(qw, type="rgb")
}
 
plot(rgb2neg(y))

---

Pixelation

pixmatr <- function(a, n){
 aa <- seq(1,dim(a)[1],n)
 ll <- seq(1,dim(a)[2],n)
 
 for(i in 1:(length(aa)-1) ){
  for(j in 1:(length(ll)-1) ){
   sub1 <- a[aa[i]:(aa[i+1]-1),ll[j]:(ll[j+1]-1)]
   k <- mean(sub1)
   sub1m <- matrix( rep(k, n*n), n, n)
   a[aa[i]:(aa[i+1]-1),ll[j]:(ll[j+1]-1)] <- sub1m
   }
  }
 
 for(j in 1:(length(ll)-1) ){
  sub1 <- a[max(aa):dim(a)[1],ll[j]:(ll[j+1]-1)]
  k <- mean(sub1)
  sub1m <- matrix( rep(k, nrow(sub1)*ncol(sub1)), nrow(sub1), ncol(sub1))
  a[max(aa):dim(a)[1],ll[j]:(ll[j+1]-1)] <- sub1m
 }
 
 for(i in 1:(length(aa)-1) ){
  sub1 <- a[aa[i]:(aa[i+1]-1),max(ll):dim(a)[2]]
  k <- mean(sub1)
  sub1m <- matrix( rep(k, nrow(sub1)*ncol(sub1)), nrow(sub1), ncol(sub1))
  a[aa[i]:(aa[i+1]-1),max(ll):dim(a)[2]] <- sub1m
 }
 
 sub1 <- a[max(aa):dim(a)[1], max(ll):dim(a)[2]]
 k <- mean(sub1)
 sub1m <- matrix( rep(k, nrow(sub1)*ncol(sub1)), nrow(sub1), ncol(sub1))
 a[max(aa):dim(a)[1], max(ll):dim(a)[2]] <- sub1m
 
a
}
 
rgb2pix <- function(img,n){
 iRed <- img[,,1]*255
 iGreen <- img[,,2]*255
 iBlue <- img[,,3]*255
 
 oRed <- pixmatr(iRed,n)
 oGreen <- pixmatr(iGreen,n)
 oBlue <- pixmatr(iBlue,n)
 
 qw <- array( c(oRed/255 , oGreen/255 , oBlue/255), dim=c(dim(iRed)[1],dim(iRed)[2],3) )
 
 imagematrix(qw, type="rgb")
}
 
plot(rgb2pix(y, 6))
plot(rgb2pix(y, 10))

---

Remove red

rgb2blu <- function(img){
 iRed <- img[,,1]
 iGreen <- img[,,2]
 iBlue <- img[,,3]
 
 oRed <- matrix(0, dim(iRed)[1], dim(iRed)[2])
 oGreen <- iGreen
 oBlue <- iBlue
 
 qw <- array( c(oRed, oGreen, oBlue), dim=c(dim(iRed)[1],dim(iRed)[2],3) )
 
 imagematrix(qw, type="rgb")
}
 
plot(rgb2blu(y))

---

Remove green

rgb2vio <- function(img){
 iRed <- img[,,1]
 iGreen <- img[,,2]
 iBlue <- img[,,3]
 
 oRed <- iRed
 oGreen <- matrix(0, dim(iRed)[1], dim(iRed)[2])
 oBlue <- iBlue
 
 qw <- array( c(oRed, oGreen, oBlue), dim=c(dim(iRed)[1],dim(iRed)[2],3) )
 
 imagematrix(qw, type="rgb")
}
 
plot(rgb2vio(y))

---

Remove blue

rgb2yel <- function(img){
 iRed <- img[,,1]
 iGreen <- img[,,2]
 iBlue <- img[,,3]
 
 oRed <- iRed
 oGreen <- iGreen
 oBlue <- matrix(0, dim(iRed)[1], dim(iRed)[2])
 
 qw <- array( c(oRed, oGreen, oBlue), dim=c(dim(iRed)[1],dim(iRed)[2],3) )
 
 imagematrix(qw, type="rgb")
}
 
plot(rgb2yel(y))

---

Adjust brightness

rgb2bri <- function(img, n){
 iRed <- img[,,1]
 iGreen <- img[,,2]
 iBlue <- img[,,3]
 
 oRed <- iRed + (iRed * n)
 oGreen <- iGreen + (iGreen * n)
 oBlue <- iBlue + (iBlue * n)
 
 qw <- array( c(oRed, oGreen, oBlue), dim=c(dim(iRed)[1],dim(iRed)[2],3) )
 
 imagematrix(qw, type="rgb")
}
 
plot(rgb2bri(y, +0.5))
plot(rgb2bri(y, -0.5))

---

Truncate colors into bands (posterize)

rgb2ban <- function(img, n){
 iRed <- img[,,1]*255
 iGreen <- img[,,2]*255
 iBlue <- img[,,3]*255
 
 band_size <- trunc(255/n)
 
 oRed <- band_size * trunc(iRed / band_size)
 oGreen <- band_size * trunc(iGreen / band_size)
 oBlue <- band_size * trunc(iBlue / band_size)
 
 qw <- array( c(oRed/255, oGreen/255, oBlue/255), dim=c(dim(iRed)[1],dim(iRed)[2],3) )
 
 imagematrix(qw, type="rgb")
}
 
plot(rgb2ban(y, 5))
plot(rgb2ban(y, 10))

---

Solarize

rgb2sol <- function(img){
 iRed <- img[,,1]*255
 iGreen <- img[,,2]*255
 iBlue <- img[,,3]*255
 
 for(i in 1:dim(iRed)[1]){
  for(j in 1:dim(iRed)[2]){
   if(iRed[i,j]<128) iRed[i,j] <- 255-2*iRed[i,j]
   else iRed[i,j] <- 2*(iRed[i,j]-128)
  }
 }
 
 for(i in 1:dim(iGreen)[1]){
  for(j in 1:dim(iGreen)[2]){
   if(iGreen[i,j]<128) iGreen[i,j] <- 255-2*iGreen[i,j]
   else iGreen[i,j] <- 2*(iGreen[i,j]-128)
  }
 }
 
 for(i in 1:dim(iBlue)[1]){
  for(j in 1:dim(iBlue)[2]){
   if(iBlue[i,j]<128) iBlue[i,j] <- 255-2*iBlue[i,j]
   else iBlue[i,j] <- 2*(iBlue[i,j]-128)
  }
 }
 
 qw <- array( c(iRed/255, iGreen/255, iBlue/255), dim=c(dim(iRed)[1],dim(iRed)[2],3) )
 
 imagematrix(qw, type="rgb")
}
 
plot(rgb2sol(y))

Inversione veloce di matrici: Strassen #2

In maniera molto simile a quanto è stato fatto per creare una funzione per eseguire velocemente moltiplicazioni tra grandi matrici con l'algoritmo di Strassen (vedi post precedente), scriviamo adesso funzioni per calcolare velocemente l'inversa di una matrice.

Per evitare di scrivere pagine e pagine di commenti e formule, come ho fatto per la moltiplicazione tra matrici, questa volta posto direttamente il codice delle funzioni create (il ragionamento che sta alla base è del tutto analogo). Copiate e incollate tutto il codice per poterlo visualizzare correttamente.

Funzione strassenInv(A)

strassenInv <- function(A){

 div4 <- function(A, r){
  A <- list(A)
  A11 <- A[[1]][1:(r/2),1:(r/2)]
  A12 <- A[[1]][1:(r/2),(r/2+1):r]
  A21 <- A[[1]][(r/2+1):r,1:(r/2)]
  A22 <- A[[1]][(r/2+1):r,(r/2+1):r]
  A <- list(X11=A11, X12=A12, X21=A21, X22=A22)
  return(A)
 }

        if (nrow(A) != ncol(A)) 
          { stop("only square matrices can be inverted") }

 is.wholenumber <-
     function(x, tol = .Machine$double.eps^0.5)  abs(x - round(x)) < tol

 if ( (is.wholenumber(log(nrow(A), 2)) != TRUE) || (is.wholenumber(log(ncol(A), 2)) != TRUE) )
   { stop("only square matrices of dimension 2^k * 2^k can be inverted with Strassen method") }

 A <- div4(A, dim(A)[1])

 R1 <- solve(A$X11)
 R2 <- A$X21 %*% R1
 R3 <- R1 %*% A$X12
 R4 <- A$X21 %*% R3
 R5 <- R4 - A$X22
 R6 <- solve(R5)
 C12 <- R3 %*% R6
 C21 <- R6 %*% R2
 R7 <- R3 %*% C21
 C11 <- R1 - R7
 C22 <- -R6
 
 C <- rbind(cbind(C11,C12), cbind(C21,C22))

 return(C)
}

Funzione strassenInv2(A)

strassenInv2 <- function(A){

 div4 <- function(A, r){
  A <- list(A)
  A11 <- A[[1]][1:(r/2),1:(r/2)]
  A12 <- A[[1]][1:(r/2),(r/2+1):r]
  A21 <- A[[1]][(r/2+1):r,1:(r/2)]
  A22 <- A[[1]][(r/2+1):r,(r/2+1):r]
  A <- list(X11=A11, X12=A12, X21=A21, X22=A22)
  return(A)
 }

 strassen <- function(A, B){
  A <- div4(A, dim(A)[1])
  B <- div4(B, dim(B)[1])
  M1 <- (A$X11+A$X22) %*% (B$X11+B$X22)
  M2 <- (A$X21+A$X22) %*% B$X11
  M3 <- A$X11 %*% (B$X12-B$X22)
  M4 <- A$X22 %*% (B$X21-B$X11)
  M5 <- (A$X11+A$X12) %*% B$X22
  M6 <- (A$X21-A$X11) %*% (B$X11+B$X12)
  M7 <- (A$X12-A$X22) %*% (B$X21+B$X22)

  C11 <- M1+M4-M5+M7
  C12 <- M3+M5
  C21 <- M2+M4
  C22 <- M1-M2+M3+M6
 
  C <- rbind(cbind(C11,C12), cbind(C21,C22))
  return(C)
 }

        if (nrow(A) != ncol(A)) 
          { stop("only square matrices can be inverted") }

 is.wholenumber <-
     function(x, tol = .Machine$double.eps^0.5)  abs(x - round(x)) < tol

 if ( (is.wholenumber(log(nrow(A), 2)) != TRUE) || (is.wholenumber(log(ncol(A), 2)) != TRUE) )
   { stop("only square matrices of dimension 2^k * 2^k can be inverted with Strassen method") }

 A <- div4(A, dim(A)[1])

 R1 <- strassenInv(A$X11)
 R2 <- strassen(A$X21 , R1)
 R3 <- strassen(R1 , A$X12)
 R4 <- strassen(A$X21 , R3)
 R5 <- R4 - A$X22
 R6 <- strassenInv(R5)
 C12 <- strassen(R3 , R6)
 C21 <- strassen(R6 , R2)
 R7 <- strassen(R3 , C21)
 C11 <- R1 - R7
 C22 <- -R6
 
 C <- rbind(cbind(C11,C12), cbind(C21,C22))

 return(C)
}

funzione strassenInv3(A)

strassenInv3 <- function(A){

 div4 <- function(A, r){
  A <- list(A)
  A11 <- A[[1]][1:(r/2),1:(r/2)]
  A12 <- A[[1]][1:(r/2),(r/2+1):r]
  A21 <- A[[1]][(r/2+1):r,1:(r/2)]
  A22 <- A[[1]][(r/2+1):r,(r/2+1):r]
  A <- list(X11=A11, X12=A12, X21=A21, X22=A22)
  return(A)
 }

 strassen <- function(A, B){
  A <- div4(A, dim(A)[1])
  B <- div4(B, dim(B)[1])
  M1 <- (A$X11+A$X22) %*% (B$X11+B$X22)
  M2 <- (A$X21+A$X22) %*% B$X11
  M3 <- A$X11 %*% (B$X12-B$X22)
  M4 <- A$X22 %*% (B$X21-B$X11)
  M5 <- (A$X11+A$X12) %*% B$X22
  M6 <- (A$X21-A$X11) %*% (B$X11+B$X12)
  M7 <- (A$X12-A$X22) %*% (B$X21+B$X22)

  C11 <- M1+M4-M5+M7
  C12 <- M3+M5
  C21 <- M2+M4
  C22 <- M1-M2+M3+M6
 
  C <- rbind(cbind(C11,C12), cbind(C21,C22))
  return(C)
 }

 strassen2 <- function(A, B){
  A <- div4(A, dim(A)[1])
  B <- div4(B, dim(B)[1])
  M1 <- strassen((A$X11+A$X22) , (B$X11+B$X22))
  M2 <- strassen((A$X21+A$X22) , B$X11)
  M3 <- strassen(A$X11 , (B$X12-B$X22))
  M4 <- strassen(A$X22 , (B$X21-B$X11))
  M5 <- strassen((A$X11+A$X12) , B$X22)
  M6 <- strassen((A$X21-A$X11) , (B$X11+B$X12))
  M7 <- strassen((A$X12-A$X22) , (B$X21+B$X22))

  C11 <- M1+M4-M5+M7
  C12 <- M3+M5
  C21 <- M2+M4
  C22 <- M1-M2+M3+M6

  C <- rbind(cbind(C11,C12), cbind(C21,C22))
  return(C)
 }

        if (nrow(A) != ncol(A)) 
          { stop("only square matrices can be inverted") }

 is.wholenumber <-
     function(x, tol = .Machine$double.eps^0.5)  abs(x - round(x)) < tol

 if ( (is.wholenumber(log(nrow(A), 2)) != TRUE) || (is.wholenumber(log(ncol(A), 2)) != TRUE) )
   { stop("only square matrices of dimension 2^k * 2^k can be inverted with Strassen method") }

 A <- div4(A, dim(A)[1])

 R1 <- strassenInv2(A$X11)
 R2 <- strassen2(A$X21 , R1)
 R3 <- strassen2(R1 , A$X12)
 R4 <- strassen2(A$X21 , R3)
 R5 <- R4 - A$X22
 R6 <- strassenInv2(R5)
 C12 <- strassen2(R3 , R6)
 C21 <- strassen2(R6 , R2)
 R7 <- strassen2(R3 , C21)
 C11 <- R1 - R7
 C22 <- -R6
 
 C <- rbind(cbind(C11,C12), cbind(C21,C22))

 return(C)
}

Proviamo a fare qualche prova. Innanzitutto controlliamo che la funzione esegua correttamente i calcoli, confrontandoli con il risultato della funzione standard di R (solve):

A <- matrix(trunc(rnorm(512*512)*100), 512,512)

all( round(solve(A),8) == round(strassenInv(A),8) )
[1] TRUE

all( round(solve(A),8) == round(strassenInv2(A),8) )
[1] TRUE

all( round(solve(A),6) == round(strassenInv3(A),6) )
[1] TRUE

La funzione esegue correttamente i calcoli. Esiste però un problema di approssimazione, per cui le prime due funzioni riportano il risultato identico fino alla ottava cifra decimale, mentre il secondo fino alla sesta. Più che un problema di calcoli, è un problema di espressione dei numeri in formato binario e 32 bit, che determina questi errori.

Ora analizziamo i tempi di esecuzione. Riporto in tabella il risultato ottenuto eseguendo il seguente codice:

Time computation

A <- matrix(trunc(rnorm(512*512)*100), 512,512)
system.time(solve(A))
system.time(strassenInv(A))
system.time(strassenInv2(A))
system.time(strassenInv3(A))

A <- matrix(trunc(rnorm(1024*1024)*100), 1024,1024)
system.time(solve(A))
system.time(strassenInv(A))
system.time(strassenInv2(A))
system.time(strassenInv3(A))

A <- matrix(trunc(rnorm(2048*2048)*100), 2048,2048)
system.time(solve(A))
system.time(strassenInv(A))
system.time(strassenInv2(A))
system.time(strassenInv3(A))

A <- matrix(trunc(rnorm(4096*4096)*100), 4096,4096)
system.time(solve(A))
system.time(strassenInv(A))
system.time(strassenInv2(A))
system.time(strassenInv3(A))

I risultati sono abbastanza evidenti, e utilizzando una modificazione dell'algoritmo di Strassen per l'inversione di matrici, esiste un effettivo risparmio di tempo.

Nuovamente valgono le due raccomandazioni già fatte:
- il codice è da migliorare; se qualcuno vuole aiutarmi, sarò lieto di aggiornare il mio codice
- se ritenete utili queste funzioni e le usate per qualche lavoro, una citazione è sempre gradita

Moltiplicazione veloce tra matrici: algoritmo di Strassen

Nel pdf qui sotto ho implementato una prima versione dell'algoritmo di Strassen in R, utile per svolgere rapidamente moltiplicazioni tra matrici.


Il codice è sicuramente da rivedere, e perfezionare, però può essere un utile inizio per una implementazione migliore.
Vorrei ringraziare l'utente G. Grothendieck per avermi suggerito un modo per creare matrici quadrate più grandi partendo da matrici più piccole (capirete leggendo il pdf), su StackOverFlow.

Se utilizzate il mio codice, è sempre cortese specificare la fonte. Inoltre se qualcuno vuole collaborare, perfezionando il codice, è assolutamente ben accetto!

Conversione numeri decimali in IEEE-754-32bit

Per un pò di teoria sullo standard IEEE-754 rimando alla pagina di Wikipedia. Qui posto solamente la funzione per effettuare la conversione in R.


Innanzitutto scriviamo alcune funzioni per eseguire la conversione dei numeri decimali in numeri binari:


decInt_to_8bit <- function(x, precs) {
q <- c()
r <- c()
xx <- c()
for(i in 1:precs){
xx[1] <- x
q[i] <- xx[i] %/% 2
r[i] <- xx[i] %% 2
xx[i+1] <- q[i]
}
rr <- rev(r)
return(rr)
}

devDec_to_8bit <- function(x, precs) {
nas <- c()
nbs <- c()
xxs <- c()
for(i in 1:precs)
{
xxs[1] <- x*2
nas[i] <- (xxs[i]) - floor(xxs[i])
nbs[i] <- trunc(xxs[i], 1)
xxs[i+1] <- nas[i]*2
}
return(nbs)
}


La prima funzione permette di convertire numeri interi (o prima della virgola) in codice binario, mentre la seconda permette di eseguire la conversione dei numeri decimali (dopo la virgola).
Il parametro precs indica il numero di bit da utilizzare.

Ad esempio:


decInt_to_8bit(11, 8)
[1] 0 0 0 0 1 0 1 1


Il numero 11, in codice binario a 8 bit è: 00001011.


devDec_to_8bit(0.625, 8)
[1] 1 0 1 0 0 0 0 0


Il numero 0.625 in codice binario a 8 bit è 10100000. Aumentando la precisione, nel primo caso permettiamo il calcolo di numeri superiori a 256 (limite massimo per gli 8-bit), mentre nel secondo caso permette di aumentare la precisione dopo la virgola:


devDec_to_8bit(0.3, 8)
[1] 0 1 0 0 1 1 0 0
devDec_to_8bit(0.3, 16)
[1] 0 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0


Ritornando ai primi due esempi, possiamo eliminare gli zeri in abbondanza con le seguenti due funzioni:


remove.zero.aft <- function(a) {
n <- length(a)
for(i in n:1){
if (a[n]==0) a <- a[-n]
else return(a)
n <- n-1
}
}

remove.zero.bef <- function(a) {
n <- length(a)
for(i in 1:n){
if (a[1]==0) a <- a[-1]
else return(a)
}
}


Abbiamo quindi:


remove.zero.bef(decInt_to_8bit(11, 8))
[1] 1 0 1 1

remove.zero.aft(devDec_to_8bit(0.625, 8))
[1] 1 0 1


Possiamo unire le due funzioni per crearne una che converte numeri decimali con la virgola in numeri binari:


dec.to.nbit <- function(x,n) {
aa <- abs(trunc(x, 1))
bb <- abs(x) - abs(trunc(x))

q <- c()
r <- c()
xx <- c()
for(i in 1:n){
xx[1] <- aa
q[i] <- xx[i] %/% 2
r[i] <- xx[i] %% 2
xx[i+1] <- q[i]
}
rr <- rev(r)

nas <- c()
nbs <- c()
xxs <- c()
for(i in 1:n)
{
xxs[1] <- bb*2
nas[i] <- (xxs[i]) - floor(xxs[i])
nbs[i] <- trunc(xxs[i], 1)
xxs[i+1] <- nas[i]*2
}

bef <- paste(remove.zero.bef(rr), collapse="")
aft <- paste(remove.zero.aft(nbs), collapse="")
bef.aft <- c(bef, aft)
strings <- paste(bef.aft, collapse=".")
return(strings)
}


Abbiamo quindi:


dec.to.nbit(11.625,8)
[1] "1011.101"


---

Imparato quindi a convertire i numeri decimali in numeri binari, possiamo passare alla conversione in IEEE-754 single precision (ossia 32bit), utilizzando la seguente funzione:


dec.to.ieee754 <- function(x) {
aa <- abs(trunc(x, 1))
bb <- abs(x) - abs(trunc(x))

rr <- decInt_to_8bit(aa, 32)

ppc <- 24 - length(remove.zero.bef(rr))

nbs <- devDec_to_8bit(bb, ppc)

bef <- remove.zero.bef(rr)
aft <- remove.zero.aft(nbs)

exp <- length(bef) - 1
mantissa <- c(bef[-1], aft)

exp.bin <- decInt_to_8bit(exp + 127, 16)
exp.bin <- remove.zero.bef(exp.bin)

first <- c()
if (sign(x)==1) first=c(0)
if (sign(x)==-1) first=c(1)

ieee754 <- c(first, exp.bin, mantissa, rep(0, 23-length(mantissa)))
ieee754 <- paste(ieee754, collapse="")

return(ieee754)
}


Il numero 11.625 in IEEE-754 è:


dec.to.ieee754(11.625)
[1] "01000001001110100000000000000000"


Per utilizzare la funzione dec.to.ieee754 servono le funzioni scritte sopra: decInt_to_8bit, devDec_to_8bit, remove.zero.bef, remove.zero.aft.

Per sicurezza potete confrontare questo output con quello che si ottiene da questo convertitore online