the strange occurrence of the one bump

When answering an X validated question on running an accept-reject algorithm for the Gamma distribution by using a mixture of Beta and drifted (bt 1) Exponential distributions, I came across the above glitch in the fit of my 10⁷ simulated sample to the target, apparently displaying a wrong proportion of simulations above (or below) one.

a=.9
g<-function(T){
  x=rexp(T)
  v=rt(T,1)<0
  x=c(1+x[v],exp(-x/a)[!v])
  x[runif(T)<x^a/x/exp(x)/((x>1)*exp(1-x)+a*(x<1)*x^a/x)*a]}

It took me a while to spot the issue, namely that the output of

  z=g(T)
  while(sum(!!z)<T)z=c(z,g(T))
  z[1:T]

was favouring simulations from the drifted exponential by truncating. Permuting the elements of z before returning solved the issue (as shown below for a=½)!

weird bug

Le Monde puzzle [#6]

A simple challenge in Le Monde this week: find the group of four primes such that any sum of three terms in the group is prime and the overall sum is minimised. Here is a quick exploration by simulation, using the schoolmath package (with its imperfections):

A=primes(start=1,end=53)[-1]
lengthA=length(A)

res=4*53
for (t in 1:10^4){

 B=sample(A,4,prob=1/(1:lengthA))
 sto=is.prim(sum(B[-1]))
 for (j in 2:4)
 sto=sto*is.prim(sum(B[-j]))

 if ((sto)&(sum(B)<res)){
 res=sum(B)
 sol=B}
 }
}

providing the solution 5 7 17 19.

A subsidiary question in the same puzzle is whether or not it is possible to find a group of five primes such that any sum of three terms is still prime. Running the above program with the proper substitutions of 4 by 5 does not produce any solution, even when increasing the upper boundary in A. So it is most likely that the answer is no.

bug in schoolmath

Neil Gunther has pointed out on his blog that the prime number decomposition R package schoolmath contains mistakes in the function primes, listing 1 as a prime number but also including decomposable numbers like 133 in its list of prime numbers:

> primes(100,140)
[1] 101 107 111 113 123 129 131 137
> primes(50,140)
[1]  51  53  59  61  67  71  73  79  83  89  97 101 103 107 109 113 127 131 133
[20] 137 139
> is.prim(primes(133)
[1] FALSE
> is.prim(primes(200,300))
[1] FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE TRUE
[13] TRUE TRUE TRUE TRUE TRUE TRUE
> sum(1-is.prim(primes(1,1000)))
[1] 10
> data(primlist)
> sum(1-is.prim(primlist[1:25000]))
[1] 3309

This is rather annoying and I hope it gets quickly fixed!