what is the probability distribution of this random walk where you have a chance to add 1 or go to one of your factors.

Question

You start at $1$.

If you're at $1$ or a prime number you can only $+1$

Let's say you're on the number $12$ there is a $50$% chance to go to $12+1$ or $13$ and a $50$% chance to go to one the factors of $12$ not including $1$ or $12$ $(2,3,4,6)$ so there is a $12.5$% to go to $4$.

A possible walk would be $1,2,3,4,2,3,4,5,6,3,4,5,6,7,8,4,5,...$ after $4$ you have a $50$% chance to go to $5$ and a $50$% chance to go to $2$ because the only factor that isn't $1$ or $4$ is $2$.

$P_n(x)$ is the probability after $n$ steps you land on $x$.

My question is what is

$$f(k) = \lim_{n \to \infty} P_n(x)?$$

Answer 1

I doubt there is going to be anything approaching a closed form, but it is easy enough to find the limiting distribution.

You can say

• State $1$ is transient - you can never get there after the first step
• All the other states are recurrent
• The chain is aperiodic (e.g. from $6$ you can go to $2$ or $3$, which breaks any periodicity)
• for large $x$ you have $f(x+1)\approx f(x)$ when $x$ is prime and $f(x+1)\approx \frac12f(x)$ otherwise
• All of states $2$ upwards are ergodic with stable probabilities which fall quickly

So we can calculate these probabilities accurately even if we ignore the very low probability tail. For example, using R

nontrivialfactor <- function(n){
  which( n %% (1:(n-1)) == 0)[-1]
  }
maxn <- 500
transitionmatrix <- matrix(0,nrow=maxn,ncol=maxn)
for (n in 1:(maxn-1)){
  pf <- properfactor(n)
  numpf <- length(pf)
  if (numpf == 0){
    transitionmatrix[n,n+1] <- 1
    }else{
    transitionmatrix[n,n+1] <- 1/2
    transitionmatrix[n,pf]  <- 1/(2*numpf)
    }
  }
matrixpower <- transitionmatrix
for (i in 1:20){
  matrixpower <- matrixpower %*% matrixpower
  }

where (as checks)

> nontrivialfactor(60)  # factors of 60 other than 1 and 60 
 [1]  2  3  4  5  6 10 12 15 20 30
> transitionmatrix[1:10,1:10]    # top left corner of transition matrix
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]    0 1.00 0.00 0.00 0.00    0  0.0    0  0.0   0.0
 [2,]    0 0.00 1.00 0.00 0.00    0  0.0    0  0.0   0.0
 [3,]    0 0.00 0.00 1.00 0.00    0  0.0    0  0.0   0.0
 [4,]    0 0.50 0.00 0.00 0.50    0  0.0    0  0.0   0.0
 [5,]    0 0.00 0.00 0.00 0.00    1  0.0    0  0.0   0.0
 [6,]    0 0.25 0.25 0.00 0.00    0  0.5    0  0.0   0.0
 [7,]    0 0.00 0.00 0.00 0.00    0  0.0    1  0.0   0.0
 [8,]    0 0.25 0.00 0.25 0.00    0  0.0    0  0.5   0.0
 [9,]    0 0.00 0.50 0.00 0.00    0  0.0    0  0.0   0.5
[10,]    0 0.25 0.00 0.00 0.25    0  0.0    0  0.0   0.0
> sum(matrixpower[1,] )    # should sum to 1
[1] 1

and we get the start of the stable distribution, i.e. $f(x)$, below.

> t(t(matrixpower[1,1:20]))  
              [,1]
 [1,] 0.0000000000
 [2,] 0.1591341850
 [3,] 0.2040682692
 [4,] 0.2198374810
 [5,] 0.1140962021
 [6,] 0.1150784653
 [7,] 0.0584909263
 [8,] 0.0586471712
 [9,] 0.0293821620
[10,] 0.0147199935
[11,] 0.0073744718
[12,] 0.0073768769
[13,] 0.0036902425
[14,] 0.0036904680
[15,] 0.0018453091
[16,] 0.0009227109
[17,] 0.0004613836
[18,] 0.0004613856
[19,] 0.0002306963
[20,] 0.0002306966

This pattern is much as expected. $4$ is the most likely value, which makes sense as the first composite number, as well as being the median, and the mean sum(matrixpower[1,]*(1:maxn)) is $4.703738$.