**infoRύχος » R**, and kindly contributed to R-bloggers)

This is one of my favourites: in 1840 the German mathematician Dirichlet proved an elegant theorem, known as “Dirichlet’s approximation theorem“. The proof is surprisingly simple, but the usefulness of the proposition in some fields of mathematics, such as Diophantine analysis is remarkable. It goes as follows:

Let a be a real number and N be a positive integer. Then, there exist two integers p and n with 1 <= n < N such that |na – p| < 1/N

It is clear that without loss of generality, the value of a can be assumed in (0,1). In this case also 1 <= p < N. So we are actually looking for integers 1 <= n,p <= N-1, or (N-1)^2 pairs. At least one of them, says Dirichlet, satisfies the above inequality. Now let’s see how many there are, or how far from “Exactly one” is this “at least one”. This is something Dirichlet could not have done by hand. The following R code is for N = 3 to 100. It is quite slow, probably because of the many for-loops.

rep <- 10^4 a <- runif(rep) freq.1 <- c() freq.2 <- c() for (N in 3:100) { freq.n <- c() for (k in 1:rep) { count.n <- c() for (j in 1:(N-1)) { count.n[j] <- sum((abs(j*a[k]-(1:(N-1))) < 1/N)) } freq.n[k] <- sum(count.n) } freq <- c() freq.1[N-2] <- sum(freq.n == 1)/rep freq.2[N-2] <- sum(freq.n == 2)/rep } png(filename = "plot.png", width = 400, height = 300) plot(3:100, freq.1, type="l", xlim=c(2,101), main="frequency plot", xlab = "integer N", ylab="frequency", ylim=c(0,1), lwd=2) lines(3:100, freq.2, type="l", cex=1.5, lty=2, lwd=2) legend(0, 1, c("exactly 1", "exactly 2"), cex=1, lty=c(1,2), col=rep("black",2), lwd=rep(2,2)) dev.off()

It seems that the two frequencies converge; in about 48% of the cases there exists exactly one pair (n,p) satisfying the inequality and in about 39% of the cases there exist exactly two. I wonder if it is possible to work with these two probabilities theoretically.

**leave a comment**for the author, please follow the link and comment on his blog:

**infoRύχος » R**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...