Project Euler 46

Project Euler 46.

Christian Goldbach apparently proposed, that every odd composite number, can be written as the sum of a prime, and twice a square. As in:

9 = 7 + 2 * 12

I had to look up what a composite number is. It is basically a number that is not a prime.

The conjecture was false. And we should now find the smallest, odd, composite number, that can not be written as the sum of a prime and twice a square.

I was not able to find a shortcut, and apparently no one else has.

So I brute forced it.

I started with loading the trusted numbers, and purrr libraries:

library(numbers)
library(purrr)

Then I guessed, correctly as it turned out, that the number I’m looking for, is less than 100000.

That means, that the prime that can take part in the sum, must also be less than 100000. So – generating two vectors, one of primes below 100000, and candidate odd composite numbers below 100000:

prim <- Primes(100000)
test <- seq(3,100000,by=2)

1 is not included. I’m not sure if it is ruled out for not being a composite. But it is trival that it cannot be a sum of a prime and any thing else (positive). So I’m starting at 3.

Actually test does not contain only composite numbers. Getting rid of the primes is the next step:

test <- setdiff(test, prim)

setdiff() is rather useful. It takes two sets – vectors, and removes all the elements in the last from the first.

Next up is getting candidate squares.
As I’m squaring them, I don’t need to look at squares of numbers larger than 707. Since 708 squared and doubled is larger than 100000.

squares <- 1:707
squares <- 2*squares**2

Next I’m writing a function to test if any x can be expressed as the sum of twice a square and a prime.

x – squares will give me a vector with the result of substracting each of the elements in squares from x. I do not need to consider values in squares that are larger than x. Subtracting those from x would give a negative number – and that is not a prime.

To do that:

x – squares[squares<x]

Next, are the elements in that vector prime? The easy way is to test if they are in the list of primes. I use the %in% operator for that:

(x – squares[squares<x]) %in% prim

That will return a logical vector, with a value TRUE or FALSE for each element in the vector resulting from the subtraction, depending on whether it is in the list of primes. If any element is TRUE, we have a result.

The any() function does that. If any element in the vector is true, any() will return true.

Putting all of that into a function:

testing <- function(x){
  any((x - squares[squares<x]) %in% prim)
}

I now have a function that returns true, if a number x can be written as the sum of a prime and twice a square.

Finally I filter that through the discard() function – basically getting rid of all the numbers in my test-vector, that can be written as this sum. And finally passing those values to the min() function. And I have my answer:

answer <- test %>%
  discard(function(x) testing(x)) %>%
  min()

Project Euler 38

Multiplying 192 with 1, 2 and 3:

192 x 1 = 192

192 x 2 = 384

192 x 3 = 576

Concatenating the products, we get 192384576. That number is 1-9 pandigital (as in all digits 1 to 9 is represented in the number exactly one time).

The same can be achieved by multiplying 9 with 1, 2, 3, 4 and 5. That gives us the pandigital concatenation 918273645.

We should now find the largest 1-9 pandigital number, formed as a concatenation in the described way, by multiplying an integer with 1, 2,…n, where n > 1.

Testing with 918273645. Nope, that is not the solution. That would be too easy 🙂

OK. Lets call the integer we’re multiplying, m.

And as we’re using R, lets call the numbers we are multiplying it with, for n, as in the vector 1:n.

The concatenated product will have to have 9 as its first digit. The first digits comes from multiplying m with 1. Therefore, the first digit in m must be 9.

We need to end up with a total of 9 digits in the concatenated product. If m has two digits, multiplying with 1 will give us two digits in the result. Multiplying with 2 will give us 3 digits in the result (any two digit number larger than 50 will, multiplied by 2, give a three digit result. And we know that m should begin with 9).

Multiplying with 3 will give another 3 digits in the result, bringing us to at total of 8. Multiplying with 4 brings the total number of digits in the concatenated result to 11. There can be only 9.

Chosing a three digit m gives similar problems.

But a four digit m, with 9 as the first digit, gives us four digits multiplying by 1 and five when we multiply by 2. Giving us a total of 9 digits in the result.

What we know so far:

  • The integer we have to multiply has four digits, and begins with 9.
  • We will have to multiply it by 1 and 2. Or 1:2.

That should be possible to brute-force.

Multiply all four digit integers where the first digit is 9 by 1:2.

Collapse the results to a single string.

Keep all the mumbers that are pandigital.

Take the maximum of it.

I don’t need to check all the numbers. I know that the first four digits must be larger than 9182, since this is the first four digits in the example. The highest number is not 9999, it is 9876.

t <- 9183:9876

Using the purrr library:

library(purrr)
answer <- t %>%
  keep(function(x) length(unique(setdiff(unlist(strsplit(paste(x*1:2,collapse=""),split="")),0)))==9)

Ok, lets stop there. What am I doing?

I pass all the possible four digit integers to keep(). keep() applies a function to each of the integers, and keeps the ones where the function evaluates to true.

The function in question is an anonymous function of x:

length(unique(setdiff(unlist(strsplit(paste(x*1:2,collapse=“”),split=“”)),0)))==9

x, the integer from t, is multiplied wiht 1:2. That results in a vector where the first element is x multiplyed by 1, and the second x multiplied by 2.

That is then put into the paste() function, with the parameter collapse=“”. The result is that the vector from just before, is concatenated to just one string.

That string is splitted by strsplit into individual characters. The result is a list. Which is bloody annoying to work with.
So I put everything into unlist(), that returns a vector.

Setdiff(x,0) removes any zeroes from that vector.

Having removed the zeroes, unique() returns all the unique values in the vector. If there is 9, the number is 1 to 9 pandigital.

And if it is, the integer is kept by keep.

It doesn’t really give me the result. It returns the four digit integers beginning with 9, that, multiplied by 1:2 and concatenated gives a pandigital result.

But the largest of them is the integer I should multiply by 1:2. So lets expand the code a bit:

answer <- t %>%
  keep(function(x) length(unique(setdiff(unlist(strsplit(paste(x*1:2,collapse=""),split="")),0)))==9) %>%
  max() %>%
  (function(x) x*1:2) %>%
  (function(x) paste(x, collapse="")) %>%
  as.numeric()

max() gives me the largest of the four digit integers. The next function mulitplies it by 1:2. After that, the result is collapsed. And finally converted to numeric by as.numeric(). Not really necessary, but I’m doing it anyway.

Lessons learned

When using pipes to do stuff, anonymous functions are useful. But should be put into parantheses.

Project Euler 32

Project Euler 32

More pandigitality. Or whatever it’s called.

39 x 186 = 7254

And that is interesting, because 39 186 7254 is pandigital. As in all the digits 1 to 9 is represented exactly once.

The task is: Find all the products, where multiplicand (39), multiplier (186) and product (7254) is 1 through 9 pandigital.

First question I need to answer is:

How many digits can there be in the multiplicand and multiplier?

If there is 1 digit in the multiplicand, there has to be 4 digits in the multiplier, and 4 in the product.

If there is only 3 in the multiplier, we would have to have 5 digits in the product, which is impossible when we multiply a 1 digit number with a 3 digit number.

And if there was 5 digits in the multiplier, we would get 5 digits in the product, bringing the total number of digits in the identity to 11.

Therefore: If there is 1 digit in the multiplicand, there has to be 4 digits in the multiplier. We will get some products with 5 digits, but we’ll filter them out.

What if there are 2 digits in the multiplicand?

Then there has to be 3 digits in the multiplier, and 4 in the the product.

If there are 2 digits in the multiplier, the highest number of digits we can get in the product, is 4. And then we would only have 8 digits in total.

If on the other hand we have 2 digits in the multiplicand, and 4 digits in the multiplier, we would get 5 digits in the product, bringing the total number of digits to 11.

So. There must be 1 or 2 digits in the multiplicand. And 4 or 3 digits respectively in the multiplier.

And there must never be more than 4 digits in the product.

Lets begin by making all possible combinations of multiplicands and multipliers.

multiplicand <- 1:9
multiplier <- 1023:9876
res1 <- expand.grid(multiplicand,multiplier)

multiplicand <- 10:98
multiplier <- 102:987
res2 <- expand.grid(multiplicand,multiplier)

testcases <- rbind(res1,res2)

All possibilities for 1 digit in the multiplicand, combined with all possibilities for a 4-digit multiplier.

Combined (rbind) with all possibilities for 2-digit multiplicands, and all possibilites for 3-digit multipliers.

Next I’ll need a function that returns TRUE if its input it 1 throug 9 pandigital:

isPan <-  function(x){ (length(unique(setdiff(unlist(strsplit(x,split="")),0)))==9)&&(nchar(x)==9) }
isPan <- Vectorize(isPan)

I also need it to be vectorized. This is not at perfect way to do it. But I can’t figure out another way.

And now we’ll use the lovely pipes to pass all the test cases through a series of functions:

library(dplyr)

answer <- testcases %>%
  mutate(product= Var1*Var2) %>%
  filter(product <= 9999) %>%
  mutate(concat = paste(Var1, Var2, product, sep="")) %>%
  filter(isPan(concat)) %>%
  select(product) %>%
  unique() %>%
  sum()

First I make a new column with the mutate() function, called product, that is the multiplication of Var1 and Var2 (multiplicand and multiplier respectively)

Next I filter out all the product that has more than 4 digits.

Next I make another new column, where I paste together the multiplicand, the multiplier and the product.

I filter, using my isPan() function. If the concatenation is pandigital, I keep it. Otherwise I throw the row away.

Now I only have rows where the concatenation of multiplicand, multiplier and product is pandigital.

But I’m actually only interested in the product, so I select the column “product”

Some of the products are repeated, unique takes care of that.

And then I just need to sum them.

Lessons learned

Vectorized functions are essential for this methodology. But can be difficult to write.

Project Euler 44

Project Euler 44 – Pentagon numbers

A pentagonal number is generated by the formula:

P~n~ = n(3n-1)/2

I must now find the pair of pentagonal numbers P~j~ and P~k~ where the sum P~j~ + P~k~ is also pentagonal.

And where P~j~ – P~k~ is pentagonal.

And where D = |P~j~ – P~k~| is minimised.

The answer is D.

First I’ll make a function generating pentagonal numbers:

pent <- function(x){
  x*(3*x - 1)/2
}

Nice and vectorized.

I can check it by calculating the first ten pentagonal numbers:

t <- 1:10
pent(t)
##  [1]   1   5  12  22  35  51  70  92 117 145

And see that they are the same as the first ten pentagonal numbers listed in the problemt text.

Next, lets make a nice long list of pentagonal numbers:

t <- 1:10000
t <- pent(t)

Next I’ll make all the pair-wise combinations of these numbers:

s <- expand.grid(t,t)

And then I’ll filter those combinations:

library(dplyr)

answer<- s %>%
  filter((Var1 + Var2) %in% t) %>%
  filter(abs(Var1 - Var2) %in% t) %>%
  transmute(D = Var1 - Var2) %>%
  abs() %>%
  min()

First – only retain the combinations where the sum of the two pentagonal numbers is itself in the list of pentagonal numbers.

Then, only retain the combinations where the absolute difference is also in the list of pentagonal numbers.

Then calculate a new row D, as the difference.

Pass that to abs() to get the absolute difference.

And finally to min() to get the minimum.

The final two steps are actually unnecessary, since I only find one pair of pentagonal numbers that meets the requirement.

Lessons learned

Hm. Not really. Why is the first I find the solution? Noone else seems to know.

Stacked barcharts with ggplot2

Barcharts to percentages

Maybe not the most correct headline.

Anyway. Assume we have five cases, with two values for a given treatment, each labelled “AVO” and “INT”. How to get a barchart, where each cases is one bar, and the bars are coloured based on the percentages of the two values?

Lets make some sample data:

set.seed(47)
cases <- 1:5
cases <- c(cases, cases)
treat1 <- rep("AVO", 5)
treat2 <- rep("INT", 5)
treat <- c(treat1, treat2)
val <- runif(10)

df <- data.frame(cases=cases, treatment = treat, value = val, stringsAsFactors = FALSE)

How does that look?

df
##    cases treatment     value
## 1      1       AVO 0.9769620
## 2      2       AVO 0.3739160
## 3      3       AVO 0.7615020
## 4      4       AVO 0.8224916
## 5      5       AVO 0.5735444
## 6      1       INT 0.6914124
## 7      2       INT 0.3890619
## 8      3       INT 0.4689460
## 9      4       INT 0.5433097
## 10     5       INT 0.9248920

OK, nice and tidy data.

What I want is a bar for case 1, filled with percentages. One color with 0.976/(0.976+0.691) and another color for 0.691/(0.976+0.691).

library(ggplot2)
ggplot(df)+
  geom_col(aes(y=value, x=cases, fill=treatment), position = "fill")

plot of chunk unnamed-chunk-3

We’re using geom_col. That is a shorthand for geom_bar with stat=“identity”, and is the shortcut to expressing values in the barchart.

The y-value is the value from df, the x-value is the cases. And we’re colouring the bars based on the value in treatment. The position=“fill” stack the bars, and normalize the heights. If we wanted to just stack the elements, we would use position=“stack”:

ggplot(df)+
  geom_col(aes(y=value, x=cases, fill=treatment), position = "stack")

plot of chunk unnamed-chunk-4

This is a great, like really great, cheatsheet for ggplot:

https://www.rstudio.com/wp-content/uploads/2015/03/ggplot2-cheatsheet.pdf

Udgivet i R

Project Euler 61

Project Euler 61

OK, this was a bit of a challenge.

Find six four-digit numbers that are cyclical. As in: The last two digits of the first number is equal to the first two digits of the second number, the last two digits for the second number is equal to the first two digits in the second number. Etc. The last two digits in the final sixth number must be equal to the first two digits in the first number.

Also, one number must be triangle, one square, one pentagonal, one hexagonal, one heptagonal and on octagonal.

Lets begin by making lists of the different kinds of numbers:

Four-digit triangle-numbers:

P3 <- function(x){
  x*(x+1)/2
}
L3 <- P3(45:140)

Square numbers:

P4 <- function(x){
  x*x
}
L4 <- P4(32:99)

Pentagonal:

P5 <- function(x){
  x*(3*x -1)/2
}
L5 <- P5(26:81)

Hexagonal

P6 <- function(x){
  x*(2*x-1)
}
L6 <- P6(23:70)

Heptagonal:

P7 <- function(x){
  x*(5*x-3)/2
}
L7 <- P7(21:63)

Octagonal:

P8 <- function(x){
  x*(3*x-2)
}
L8 <- P8(19:58)

Lets get them all into a list:

values <- list(L8, L7, L6, L5, L4, L3)

And this is where it gets complicated.

I take each value in L8, which is now values[[1]]

For each of them, I run through each of the other elements of values, 2:6.

For each values[[i]], I check if there are any of the numbers, where the first two digits is equal to the last two in the element in L8 I’ve come to.

If there is, I run through all of them (the ones where the first two digits is equal to the last two in the element in L8 I’ve come to.)

For each of these, I run through each of the remaining list elements, setdiff(2:6, b), where b is the iterator in 2:6

For each of those, I check that there are numbers where the first two digits etc etc etc.

When I get to the last number, I check that the last two digits are the same as the first two digits in the number in the first loop. Then I sum the numbers I’ve found, and get the answer.

It is not difficult. As such. But there are a lot of nested for-loops, and I have to check that there is a possible next number, and, if not, go to the next element.

It took quite a lot of time to get it right…

for(i in values[[1]]){
  for(b in 2:6){
    if(length(values[[b]][(values[[b]] %/% 100) == (i %% 100)])==0){
        next()
      }
    for(j in values[[b]][(values[[b]] %/% 100) == (i %% 100)]){
      for(c in setdiff(2:6, b)){
       if(length( values[[c]][(values[[c]] %/% 100) == (j %% 100)]  )==0){
         next()
       }
       for(k in values[[c]][(values[[c]] %/% 100) == (j %% 100)]){
         for(d in setdiff(2:6, c(b,c))){
           if(length( values[[d]][(values[[d]] %/% 100) == (k %% 100)]  )==0){
         next()
       }
           for(l in values[[d]][(values[[d]] %/% 100) == (k %% 100)]){
             for(e in setdiff(2:6, c(b,c,d))){
               if(length( values[[e]][(values[[e]] %/% 100) == (l %% 100)]  )==0){
         next()
       }
               for(m in values[[e]][(values[[e]] %/% 100) == (l %%100)]){
                  for(f in setdiff(2:6, c(b,c,d,e))){
                    if(length( values[[f]][(values[[f]] %/% 100) == (m %% 100)]  )==0){
                    next()
                     } else{
                    for(n in values[[f]][(values[[f]] %/% 100) == (m %% 100)]){
                          if((n %% 100) == (i %/% 100)){
                        answer <- sum(i,j,k,l,m,n)
                        break()
                      }}
                    }
                  }   
               }
             } 

           }
          }
       }
     }
    }
  }
}

Neat. But there must be an easier way.

Project Euler 31

How many ways can you make 2 pounds given these sizes of coins:

1p, 2p, 5p, 10p, 20p, 50p, 1£ and 2£

Where 1£ is 100p.

We can take either 0 or 1 2£ coin.

answer <- 0

for(a in seq(0,200,by=200)){
  for(b in seq(0,(200-a),by=100)){
    for(c in seq(0,(200-a-b), by=50)){
      for(d in seq(0,(200-a-b-c),by=20)){
        for(e in seq(0,(200-a-b-c-d), by=10)){
          for(f in seq(0,(200-a-b-c-d-e), by=5)){
            for(g in seq(0,(200-a-b-c-d-e-f),by=2)){
              answer <- answer + 1
            }
          }
        }
      }
    }
 }
}

How to explain it…

I need to count up to 200p total.

a is the part of the 200p that is made up by a 2£ coint. That is either 0 or 200. Or the sequence from 0 to 200 taken in increments of 200.

What is left after considering 2£ coins, is 200 minus the value of a. The number of 1£ coins we can use to make up that amount is all the values from 0 to what is left. Ie 0 to 200-a in increments of 100. The part of the total 200p that is made up by 1£ coins is b.

After adding some 1£ coins, there are 200-a-b left. That can be made by 50p coins. The part of what is left i c, and that is values from 0 to 200-a-b in increments of 50, ie seq(0,(200-a-b), by=50)

Etc.

The trick is to make sure that the ranges of the inner loops is controlled by what happened in previous loops.

Otherwise I will loop over far too many values.

Project Euler – 28

Number spiral diagonals

We can make squares of numbers, by starting with 1, and add numbers in a spiral.

It is easier to show than descripe:

square <- matrix(c(21, 22, 23, 24, 25,20,  7,  8,  9, 10,19,  6,  1,  2, 11,18,  5,  4,  3, 12,17, 16, 15, 14, 13), 5,5,byrow=T)
square
##      [,1] [,2] [,3] [,4] [,5]
## [1,]   21   22   23   24   25
## [2,]   20    7    8    9   10
## [3,]   19    6    1    2   11
## [4,]   18    5    4    3   12
## [5,]   17   16   15   14   13

This is a 5 by 5 square. Start at 1, add the 2 to the right of that, and add numbers in a clockwise spiral.

If we look at all the numbers in the diagonals, we can see that the sum of those numbers, is 101. Ie, 21 + 7 + 1 + 3 + 13 + 25 + 9 + 5 + 17 = 101

We are now asked: What is the sum of the numbers in the diagonals in a 1001 by 1001 square?

Looking at an n by n square, where n is uneven, we can see that the north-eastern corner of the square will always be n2.

We know that the side is n. So the value of the north-western corner will always be n2 – (n-1).

Similar the south-western corner is n2 – 2(n-1). And the south-eastern is n2 – 3(n-1).

The first “square” is a special case. But I can now calculate the corners in all the squares with unequal n.

The sum of those corners is:

n2 + n2 -(n-1) + n2 – 2(n-1) + n2 – 3(n-1)

Which reduces to:

4*n2 – 6n + 6

If I calculate that for all unequal n from 3 to 1001, sum them, and add 1, I have my answer.

library(purrr)
answer <- seq(3,1001,by=2) %>%
  map(function(x) 4*x*x - 6*x + 6) %>%
  unlist() %>%
  sum() 
answer <- answer + 1

Lessons learned

Not really.

Euler 36

Some numbers are palindromic, ie they are the same read in both directions.
Some numbers are palindromic in more than one base.
The problem mentions 585, which is palindromic. In binary 585 is 1001001001. And that is palindromic as well.
Find all numbers below 1000000 that are palindromic in both base 10 and base 2.
What is their sum?

First I’m going to need a function that checks if a number is palindromic.

palindromic <- function(x){
  paste(rev(unlist(strsplit(as.character(x),""))), collapse="")==as.character(x)
}

This function converts the input to character, splits it in individual characters, and unlists it. Because strsplit() returns a list. Then that vector is reversed with rev(), and pasted together to a single string. That string is then compared to the input – converted to character.
It returns TRUE if the input is palindromic. And what is nice is that it doesn’t care if the input is a number or a string.

Next I’m going to need a function that converts a number to it’s binary representation.
There is one, intToBits(). But it returns all the individual bits. And in reverse order. Also, it returns the bits as characters. So. intToBits(x) to get the bits. as.integer() to convert them to integers. Then the vector is reversed. And collapsed with paste().
One more detail – leading zeroes are not allowed. Therefore I remove all leading zeroes with a regular expression.

binary <- function(x){
  unlist(map(x, function(x) sub("^[0]+", "", paste(rev(as.integer(intToBits(x))), collapse="")))         )
}

This is a bit annoying. paste() pastes everything together. When I pass a vector to it, it will not only split the individual elements apart, convert them, and then paste them together individually. It will also paste together every single element in the vector, and just return a very long string.
map() solvse this. I define the “split, convert, paste-together” part as an anonymous functtion, which is then passed to map along with the vector. map() applies the function to each element in the vector, and as if by magic, stuff is no longer pasted together. map() returns a list, so I’ll have to wrap it in unlist().

Now I have all the elements. I just need find a way to convert a binary number to its decimal representation.
strtoi() does that. It takes a character vector, and converts it to integer – given a certain base, in this case 2.
I’ll just use that as an anonymous function.

Piping everything together:

t <- 1:1000000
s <- t %>%
  keep(palindromic) %>%
  binary() %>%
  keep(palindromic) %>%
  (function(x) strtoi(x, base=2)) %>%
  sum()

Passing all numbers from 1 to a million to keep(), with the predicate function palindromic. There are faster ways to generate palindromes, I know.
Piping the result of that to binary, converting every palindromic number (decimal) to binary.
Then, once more, through keep, only keeping the binary numbers which are palindromic.
Now we have the numbers which are palindromic in both decimal and binary. They are converted to integers, and summed.

Lessons learned.
1. map() is good to know – The alternative here would be to use Vectorize().
2. How to convert from decimal to binary, and back.

Euler 42

A triangle number is given by t~n~ = ½(n(n+1))

Lets start by writing a function for that:

triangle <- function(n){
  (n*(n+1))/2
}

If we take a word, Project Euler uses SKY as an example, we can take each letter, assign to it the value of its position in the alphabet, and add together the values of the letters.
SKY becomes 19 + 11 + 25 = 55.

Lets write a function that finds that:

lettervalue <- function(chars){
  chars <- toupper(chars)
  t <- unlist(strsplit(chars,""))
  sum(match(t, LETTERS))
}

First I convert the character string in the input to uppercase letters. This might not be necessary.
Then I split the string in individual characters. And unlist() the result.
Then I pass it to the match() function. That gives me the position of a character as a numerical value, in the vector LETTERS.
Neat function.
Then I sum it.

This is nice and vectorized all the way through.

But what was the problem I should solve?
Project Euler provides a list of words. How many of these are, when converted to a sum with the function above, a triangle number?

I have no idea about how long a given word is. But lets find out.

First, lets read in the file.

wordlist <- scan("p042_words.txt", what=character(), sep=",", quiet=TRUE)

What are the numerical values as defined above?
My function lettervalue just adds everything together. Which is a bit unhelpful. If I use map() from the purrr-libray, I can apply the function to each element in wordlist. I’ll have to unlist() it afterwards. That gives me all the numerical values:

library(purrr)
values <- unlist(map(wordlist, lettervalue))

What is the largest of them?

max(values)
## [1] 192

That is the highest number I have to check is triangular.

Lets generate a list of the first 20 triangle numbers:

t <- 1:20
triangles <- triangle(t)
max(triangles)
## [1] 210

That was lucky, triangle-number 20 is 210, which is larger than the largest value I have in the set of words.
Now I just need to find out how many of the values in values are in the variable triangles:

answer <- sum(values%in%triangles)

Neat.
Lessons learned:
1. map() can be used to apply a function to a vector. Even functions that are not that well behaved.
And it’s pretty damn fast. An alternative would be one of the apply-functions:

library(microbenchmark)
map_test <- function(){map(wordlist, lettervalue)}
apply_test <- function() {lapply(wordlist, lettervalue)}
mbm <- microbenchmark(map_test, apply_test, times=1000)
library(ggplot2)
autoplot(mbm)
## Warning: Transformation introduced infinite values in continuous y-axis
## Warning: Removed 915 rows containing non-finite values (stat_ydensity).

plot of chunk unnamed-chunk-8

  1. match() is also pretty useful. I spend a couple of frustrated minutes trying to get a list to behave as a lookup-table. match() is much easier.