I wrote a code to calculate a series of summations and products for a maths assignment I'm doing. The equation is shown below.$$E(\sigma)=[1-(P(\text{$1^{st}$ rejects}) \times P(\text{$2^{nd}$ rejects}) \times \ldots \times P(\text{$n^{th}$ rejects}))] \times \frac{\frac{n(n+1)}{2} - \frac{(n-j)(n-j+1)}{2}}{j}$$where$$P(\text{$x^{th}$ rejects})=1-\left[\frac{j}{n} \times\left(\frac{n-1}{n}\right)^{j-1} \times \prod_{y=2}^{x-1}\sum_{k=0}^{j-1} (\binom{j-1}{k}(\frac{n-1}{n})^{j-1-k}(\frac{1}{n})^k(1-\prod_{i=0}^{Y-1} P(i^{th} rejects))) \times \\ \sum_{k=0}^{j-1}(\binom{j-1}{k}(\frac{n-1}{n})^{j-1-k}(\frac{1}{n})^k \sum_{i=0}^{k}(\frac{1}{1+i}\binom{k}{i}(1-\prod_{i=0}^{x-1} P(i^{th} rejects))^{k-i}(\prod_{i=0}^{Y-1} P(i^{th} rejects))^i)))\right]$$The code calculates \$E(\sigma)\$ for all values of \$j\$ for \$1\le j\le n\$ for all values of \$n\$ for \$3\le n\le a\$. The code finds the value of \$j\$ which produces the highest value of \$E(\sigma)\$ for each \$n\$ and expresses this as a percentage of \$n\$. Then it plots the graph of this percentage versus \$n\$The thing is it requires many loops and the computation for \$a=100\$ seems like it's taking forever. Is there a way to reduce the computation time? If not, if I leave the computer running for a few hours will it eventually compute successfully? I'm really new to programmming so please excuse the messy code. \$P(\text{$1^{st}$ rejects})\$ and \$P(\text{$2^{nd}$ rejects})\$ are not calculated in the same loop as the rest of \$P(\text{$x^{th}$ rejects})\$ but rather separately as prob1r and prob2r respectively as putting them in the same loop creates some problems (such as for y in 2:1).Edit: I looked around and discovered the magic of parallelization. By simply changing my for loop to a paralllelized foreach loop, time taken decreased from Time difference of 58.93389 secs to Time difference of 16.01991 secs for \$n=40\$. Still, help will be appreciated if further reduction can be achieved as I want to calculate to as high a value of \$n\$ as possible.a=as.numeric(a)registerDoParallel(cores=3)z=matrix(nrow=1, ncol=a)y=(foreach(n=3:a, .combine='c') %dopar%{  q=matrix(nrow=1,ncol=n)  for(j in 1:n)  {    m=matrix(nrow=1,ncol=n)    prob1r=0    for (k in 0:(j-1))    {      prob1r=prob1r+(((dbinom(k,(j-1),(1/n))))*(1/(k+1)))    }    m[1,1]=1-((j/n)*prob1r)    prob2r=0    for(k in 0:(j-1))    {      probir=0      for (i in 0:k)      {        probir=probir+((dbinom(i,k,m[1,1])))*(1/(1+i))      }      prob2r=prob2r+((dbinom(k,j-1,(1/n))))*probir    }    m[1,2]=1-((j/n)*(((n-1)/n)^(j-1))*prob2r)    for(h in 3:n)    {      prob1=1      for(y in 2:(h-1))      {        prob=0        proby1=1        for(k in 1:(y-1))        {          proby1=proby1*m[1,k]        }        for(k in 0:(j-1))        {            prob=prob+((dbinom(k,j-1,(1/n))))*((1-proby1)^k)        }        prob1=prob1*prob      }      prob2=0      for(k in 0:(j-1))      {        probrej=1        for (g in 1:(h-1))        {          probrej=probrej*(m[1,g])        }        probi=0        for (i in 0:k)        {          probi=probi+((dbinom(i,k,probrej)))*(1/(1+i))        }        prob2=prob2+((dbinom(k,j-1,(1/n))))*probi      }      m[1,h]=1-((j/n)*(((n-1)/n)^(j-1))*prob1*prob2)    }    probno=1    for(p in 1:n)    {      probno=probno*(m[1,p])    }    q[1,j]=(1-probno)*((((n*(n+1))/2)-(((n-j)*(n-j+1))/2))/j)  }  z[1,n]=((which.max(q))/n)*100})x=3:alm.r=lm(y~x)plot(x,y)abline(lm.r)
Calculate and plot expected value involving summations and products
beginner;time limit exceeded;combinatorics;r