melt() – tidy data

Just a short note to help me remember the melt()-function.

Lets create some messy data:

ID <- 1:10
T1 <- runif(10,10,20)
T2 <- runif(10,20,30)
T3 <- runif(10,30,40)
df <- data.frame(ID=ID, T1=T1, T2=T2, T3=T3)

This is heavily inspired by a practical problem a student came to us with. There is 10 different patients, at time = T1, there is a certain value measured on the patient. At time = T2 the same value is measured. And again at time = T3, where T1<T2<T3.

We would now like to plot the development of those values as a function of time (or just T1,T2 and T3).

How to do that?

Using reshape2, and making sure we have the tidyverse packages loaded:

library(tidyverse)
library(reshape2)
clean <- melt(df, id.vars=c("ID"), value.names=c("T1", "T2", "T3"))
head(clean)
##   ID variable    value
## 1  1       T1 13.93662
## 2  2       T1 10.30468
## 3  3       T1 18.14351
## 4  4       T1 17.58294
## 5  5       T1 18.13877
## 6  6       T1 10.21993

Nice, we now have tidy data.

Note that melt() also takes a variable “variable.name”, in case we have a different sort of mess.

Now it is easy to plot:

plot <- ggplot(clean, aes(x=variable, y=value, group=ID)) +
  geom_line()
plot

plot of chunk unnamed-chunk-4

Neat.

Euler 100

Project Euler – problem 100

Back to the hopeless examples of probabilities from school.
In a bag there are 15 black balls and six white ones. Project Euler talks about discs, math-teachers has always used balls as examples, and they where always white and black. So I’ll stick with that.

It you draw two balls from the bag, there is a 50/50 chance of drawing 2 black balls:

(15/21)*(14/20)
## [1] 0.5

I’m told that the next set of balls in the bag with that property, is 85 black balls and 35 white ones:
(85/120)

(85/120)*(84/119)
## [1] 0.5

Find the mix of black and white balls, that gives a probability of 50/50 of drawing 2 black balls, given that there should be more than 10¹² = 1000000000000 balls in the rather large bag.

That should be straight-forward.
Lets call the number of black balls b and the number of white balls w. And lets define the total number of balls in the back as n=w+b
The probability of drawing two black balls is:

(b/n)((b-1)/(n-1)) = ½

n = w + b > 10¹²

Two equations with two unknowns.
The probability can be rearranged:

(b/n)((b-1)/(n-1)) = ½ <=>

b(b-1) / n(n-1) = ½ <=>

(b² – b) / (n² – n) = ½ <=>

b² – b = ½(n² – n) <=>

2b² – 2b = n² – n <=>

2b² – 2b – n² + n = 0

Hm. Maybe it is not that simple after all. First of all I don’t know if n is 100000000000 or 100000000001. That actually makes a pretty big difference:

1000000000001**2 - 1000000000000**2
## [1] 1.999978e+12

Second of all, I need to find integer solutions. An analytical solution might not give integer results. And I can’t have one third of a ball in the bag.

Googling “finding integer solutions to equations” give, as the first result, a link to the wikipedia article on “Diophantine equations”.
Which apparently are equations that should have integer solutions.

All right, a couple of the problems I’ve tackled earlier, and quite a lot of Project Euler problems I’ve given up on appears to be about solving these Diophantine equations.

So. Nice. The last link of the wikipedia page is to https://www.alpertron.com.ar/QUAD.HTM.
I should probably read up on the methods. But that will have to wait.

The point is, that this Diophantine equation can be solved by:

b~n+1~ = 3b~n~ + 2n~n~ -2

n~n+1~ = 4b~n~ + 3n~n~ -3

The idea is that we have a solution (b~n~, n~n~). And these two equations allows us to calculate the next solution, (b~n+1~, n~n+1~)

Lets try that, we was given that (15,21) was a solution. The next should be (85,120). Do we get that?

b <- 15
n <- 21
b_n <- 3*b + 2*n -2
n_n <- 4*b + 3*n -3
print(paste(b_n, n_n, sep=","))
## [1] "85,120"

Qap’la, it works. Nice. Now I just need to run through this until n~n+1~ gets above 10¹².

b <- 15
n <- 21
while(n<10**12){
  b_n <- 3*b + 2*n -2
  n_n <- 4*b + 3*n -3
  b <- b_n
  n <- n_n
}
answer <- b

Lessons learned:

  1. Solving Diophantine equations is at the heart of a lot of these problems. I’ve learned a new tools to handle them!
  2. If you want to subscript stuff in RMarkdown, you place a ~ on each side of what you want subscripted.

Other stuff to note: Maybe it is time someone wrote a new solver for Diophantine equations. The one I found is 19 years old. Something to do in Shiny perhaps?

Euler 80

Problem 80 from Project Euler.

The problem tells us that if the square root of a natural number is not an integer, it is irrational.
Project Euler also claims that it is well known. I did not know it.

We are then told that the square root of 2 is 1.4142135623730950… And that the digital sum of the first 100 digits is 475.

The task is now to take the first 100 natural numbers. And find the total of the digital sums for the first 100 digits for all the irrational square roots.

Lets begin by figuring out how to handle that many digits. R does not support more than around 15 places after the decimal point.

The library Rmpfr can handle arbitrary precision:

library(Rmpfr)
## Loading required package: gmp
## 
## Attaching package: 'gmp'
## The following objects are masked from 'package:base':
## 
##     %*%, apply, crossprod, matrix, tcrossprod
## C code of R package 'Rmpfr': GMP using 64 bits per limb
## 
## Attaching package: 'Rmpfr'
## The following objects are masked from 'package:stats':
## 
##     dbinom, dnorm, dpois, pnorm
## The following objects are masked from 'package:base':
## 
##     cbind, pmax, pmin, rbind
a <- sqrt(mpfr(2,500))

The variable a now contains the square root of 2 with a precision of 500 bytes. I’m not quite sure how many decimal places that actually translates to. But testing the following code allows me to conclude with confidence that it is at least 100.

A thing to note here is, that

a <- mpfr(sqrt(2),500)

and

a <- sqrt(mpfr(2,500))

are not equal. In the first exampel sqrt(2) is evaluated before saving the value with the high precision. Start by converting the number 2 to a high precision representation, before doing math on it.

Next is writing a function that will return the digital sum of the first 100 digits of a number.

digitsum <- function(x){
  s <- 0
  for(i in 1:100){
    s <- s + floor(x)
    x <- (x - floor(x))*10
  }
  s
}

First s is initialized to 0. Then floor(x) gives us the first digit in x. We add that to s, and subtract it from x, and multiply by 10. Repeat that 100 times, and you get the sum of the first 100 digits in x.

Let us test that it works. Project Euler told us what the result for sqrt(2) is:

digitsum(a)
## 1 'mpfr' number of precision  500   bits 
## [1] 475

Nice, the correct result (not that that guarantees that I’ve done everything correctly).

Now, lets find all the irrational square roots we need to look at:

library(purrr)
t <- 1:100
s <- t %>%
  keep(~as.logical(sqrt(.x)%%1))

I need to practice this way of coding a bit more. t contains the first 100 natural numbers. I pass that to the keep()-function, and the predicate function takes the square root of each number, take the modulus 1, and convert it to a logical value. If the modulus of the square root is 0, the square root is an integer, and 0 i false. So we’re keeping all the non-integer squareroots.

Now I’ll convert all the natural numbers to the mpfr-class. The next line takes the square root. The third line calculate the digitalsum. And the final line gives us the result.

s <- mpfr(s,500)
r <- sqrt(s)
r <- digitsum(r)
sum(r)
## 1 'mpfr' number of precision  500   bits 
## [1] Censored

Lessons learned:
Rmpfr allows us to work with (more or less) arbitrary precision.
But we need to convert numbers to the relevant class before doing math on it.

Replacing values in a dataframe – to what a previous value was

Given a set of data, where some values indicate that they are the same as a previous value, how to replace them with the correct value.

Eg, this dataframe:

(m <- data.frame(i=c(1:10,NA), t=c("lorem", "do", "do", "Do", "ipsum", "do", "Do", "(do)", "dolor", NA, "test"), stringsAsFactors=F))
##     i     t
## 1   1 lorem
## 2   2    do
## 3   3    do
## 4   4    Do
## 5   5 ipsum
## 6   6    do
## 7   7    Do
## 8   8  (do)
## 9   9 dolor
## 10 10  <NA>
## 11 NA  test

How to replace the first three “do”s with “lorem” and the next set of “do”s with “ipsum”

Using fill() from the tidyr package is straight forward. It takes a vector, locates all NA, and replaces them with the last, non-NA value.
Simple enough, change all the variations of “do” to NA, run fill(). Done.
One problem, there might be NAs in the dataset, that we do not want to affect.
Solution – there might be a more elegant one, but this works:

  1. Change the NAs to something that do not occur in the data
  2. Change to variations of “do” to NA
  3. Use the fill()-function
  4. Change the NAs from step 1 back to NA
library(tidyr)
rpl <- "replacement"
m[is.na(m$t),]$t <- rpl
doset <- c("do", "Do", "(do)")
m[(m$t %in% doset),]$t <- NA

m <- m %>% fill(t)
m[(m$t == rpl),]$t <- NA
m
##     i     t
## 1   1 lorem
## 2   2 lorem
## 3   3 lorem
## 4   4 lorem
## 5   5 ipsum
## 6   6 ipsum
## 7   7 ipsum
## 8   8 ipsum
## 9   9 dolor
## 10 10  <NA>
## 11 NA  test

Done!

Oh, and by the way, this is my first post generated directly from RStudio!

Der skal mere liv her!

Og et par andre steder. Noget af det jeg bruger en del tid på, både privat og professionelt, er R. Som i det statistiske program R.

Jeg skal derfor snart have taget et kig på denne side:

Og få publiceret hvad jeg alligevel nusser rundt med her og andre steder.

Hide rows, based on value of cell – in Excel

So – you want to hide some rows on a worksheet, based on the value in a cell. Or more than one.
Here’s how to do that with VBA
Find the last row of the range that you want to apply the hiding to.
Get a range of rows, in this case starting at A7, and ending at “LastRow”.
For each value in that range, if the value i column A is equal to the value in cell G1 (that is Cells(1,7), And the value three columns over, eg in colum D (that is c.Offset(0,3)), is equal to the value in cell G2 (Cells(2,7), then set the entire row to be hidden, else set it to be shown.


Private Sub Worksheet_Change(ByVal Target As Range)
Dim LastRow As Long, c As Range
Application.EnableEvents = False
LastRow = Application.WorksheetFunction.CountA(Range("A7:A100000")) + 6
On Error Resume Next
For Each c In Range("A7:A" & LastRow)
If (c.Value = Cells(1, 7).Value And c.Offset(0, 3).Value = Cells(2, 7).Value) Then
c.EntireRow.Hidden = False
Else
c.EntireRow.Hidden = True
End If
Next
On Error GoTo 0
Application.EnableEvents = True
End Sub

Euler problem 41

Project Euler is a pretty good way to exercise your programming muscles. I tend to think that the hardest part is usually the math.

So when I figure one of them out in minutes, I’m pretty happy about it.

library(gtools)
library(numbers)
pandig <- function(n){
x <- permutations(n,n)
x <- apply(x,1, function(x) paste(x, collapse=""))
return(x)
}
y <- as.numeric(pandig(7))
max(y[isPrime(y)])

The surprising thing is that largest pandigital prime does not begin with either 8 or 9.

Bitcoins and crypto currency

These days, when at least one crypto currency has tanked completely, Bitcoin is coming under increasing pressure from authorities, revelations that there are organized manipulation of the trading etc etc etc.

I am reminded of something that happened in the Netherlands in 1637.

Read about it here. I am continually baffled by our inability to learn from history.

Practical project management

Project management is not easy. A handfull of practical hints from ~20 years of experience:

  • Dont ask for more ressources. Ask for the ressources you were promised, but never got.
  • Dont ask for an extension on deadlines. Ask that the ones you got in the first place are not changed.
  • Never expect anyone to have heard of the iron triangle of project management.
    Even if they actually have heard of it before.
  • Never expect anyone to have any understanding of what the project is about.
    Even if they are in charge of managing it.
  • When, at the beginning of a project, estimations of the necessary ressources are deemed irrelevant by the steering commitee, buy antacids. You will get an ulcer.
  • When someone is given a budget – never expect them to stay within it.

 

Command line tools are awesome

Linux have a lot of small tools, that only does one thing. But do it really well (compare to Windows, that has a lot of large tools that does everything rather badly).

This is really just a note to myself. These tools are really useful but they are not (yet) second nature to me. I often find myself in the situation, where I know there is a tool for something, that I have used several times before. But simply can’t remember what it was.

grep. Searches files for lines matching a regular expression. Useful parameters (or at least parameters I have a regular use for):

-c returns a count of the lines matching.
-n returns the linenumber (in the file) of the matching line.

tail. returns the last part of a file

-n 6 returns the last 6 lines of the file (standard 10)

head. Returns the first part of a file

-n 6 returns the first 6 lines of the file (standard 10)

cut. Removes sections of lines in a file (or other input)

-d x. Splits the line at x. Use ‘ ‘ for space
-f 1. Select the first field.

wc. Counts stuff in files.

-l. counts the lines in the file (or other input)
-w count the words in the file (or other input)

|. Piping. Takes the result of the command in front of it, and pass it to the command after it (and that is the direction. If you find examples on Stackoverflow that will only give the desired result if the direction is reversed, don’t be surprised if it does not work…)

cat. Prints one or more files to standard output (your screen).

But if we print to another output, eg with “> file.name”, we can concatenate several files.

find. Searches for files. “find .” finds everything. Pipe it to grep to search for something specific. eg “find . | grep ‘acta'” to find all files containing the string “acta”.

-print prints the complete filename.
-print0 prints the complete filename even if it includes a newline.