title: R code produces different output for equation in loop every time
body: <p>I'm trying to calculate the geometric mean and find the maximum possible value for it for a certain set of parameters. I am convoluting two distributions (I think, not sure about the terms in english), one of which is the probability density function of my equation, and want to find the set of parameters where the geometric mean is the highest. I'll try to annotate my code, so that it gets clearer what I'm doing. I'm doing this to determine the optimal niche width for a species under certain cirumstances. The fitness depends on the match between the trait of the species and the standard deviation of the Environment <code>sdH</code> and a trade-off function.</p>
  
  <p>However, if I change one of the parameters, it changes the complete output, even though I would not expect it to do so. Maybe I am overlooking some mistake or this is acutally to be expected and my expectations are wrong, but I am at my wits end. Maybe someone has any idea how to do this? </p>
  
  <pre><code>deltax=0.5    #This is the step size for one of the distributions
  alpha=c(1, 1.5, 2, 3, 4)    #These are one set of parameters I want to study. They are each a distinct case.
  H=seq(-7, 7, deltax)    #This is the scope of the x-values for one distribution
  sdH=seq(1, 7, 0.02)    #This is the scope of the standard deviation I want to study
  optimumgvariance= matrix(0,ncol = length(alpha), nrow = length(sdH)) #a matrix to store my final values
  </code></pre>
  
  <p>This is the first set of parameters I'm using in the follwing code:</p>
  
  <pre><code>for (l in 1:length(sdH)){  #for each value of sdH I do the following:  
    g=seq(0.1, 9, 0.1)    # set up a vector of values, which determine the width of the other distribution
    gt=0    #initialize a helper variable
    tradeoff=numeric(length(g))     #initilaize a vector for the results of one of my equations
    expfitnessgeo=numeric(length(g))    #initiliaze a vector for the geometric mean
    optimumg=numeric(length(alpha))    #initialize a vector for the optimum value of g (the values from above) for the given set of parameters (depending on which loop we are in)
    for (k in 1:length(alpha)){    #for all values of alpha we do the following
      for (i in 1:length(g)){    #for all values of g
        gt=-(g[i]^2)    #I calculate the helper variabel gt (for clearer code futher below)
        tradeoff[i]=exp(gt/(2*alpha[k]^2))    #I calculate the results of the trade-off function for the combination of g and alpha in the current loop
        survival=numeric(length(H))    #I set up a vector to contain the results of the next equation
        for (j in 1:length(H)){    
          survival[j]=tradeoff[i]*exp(-((H[j])^2/(g[i]^2)))    #the survival
        }    
        expfitnessgeo[i]= exp(sum(log(survival)*dnorm(H, 0, sdH[l])*deltax)) # the maximum long term fitness can be identified by a maximum geometric mean of the probability density function    
      }    
      print(max(expfitnessgeo))    #just to control the output
      if(max(expfitnessgeo)&gt;=0.1){    #the geometric mean needs to be larger than 0.1 
        optimumg[k]=g[which(expfitnessgeo==max(expfitnessgeo))]  #the optimal value of g is the one where the geometric mean of the pdf is maximal.  
      } else optimumg[k]=0    #if the geometric mean of the pdf is smaller than 0.1 I set the optimum g to zero
  
    }    
    optimumgvariance[l,] =optimumg #store the results in the matrix    
  }    
  </code></pre>
  
  <p>This produces results as I would expect:</p>
  
  <pre><code>&gt; tail(optimumgvariance)
         [,1] [,2] [,3] [,4] [,5]
  [296,]    0    0  3.3  4.1  4.7
  [297,]    0    0  3.3  4.1  4.7
  [298,]    0    0  3.3  4.1  4.7
  [299,]    0    0  3.3  4.1  4.7
  [300,]    0    0  3.3  4.1  4.7
  [301,]    0    0  3.3  4.1  4.7
  &gt; max(optimumgvariance[,2])
  [1] 2.3
  &gt; max(optimumgvariance[,1])
  [1] 1.5
  </code></pre>
  
  <p>However, when I only change the scope of <code>sdH=seq(1,10,0.02)</code> and leave everthing else as it is, the results change. I suddenly have values in the second column, where I previously had the expected zeros:</p>
  
  <pre><code>deltax=0.5    
  alpha=c(1, 1.5, 2, 3, 4)    
  H=seq(-7, 7, deltax)    
  sdH=seq(1, 10, 0.02)    
  optimumgvariance= matrix(0,ncol = length(alpha), nrow = length(sdH))    
  
     &gt; tail(optimumgvariance)    
         [,1] [,2] [,3] [,4] [,5]    
  [446,]    0  2.9  3.4  4.1  4.8    
  [447,]    0  2.9  3.4  4.1  4.8    
  [448,]    0  2.9  3.4  4.1  4.8    
  [449,]    0  2.9  3.4  4.1  4.8    
  [450,]    0  2.9  3.4  4.1  4.8    
  [451,]    0  2.9  3.4  4.1  4.8    
  &gt; max(optimumgvariance[,1])    
  [1] 1.5    
  </code></pre>
  
  <p>I would expect that the values in column 2 stay as before, but somehow they don't. And this goes even further, when I put <code>sdH= seq(1,15, 0.02)</code> I even get values other than zero in the first column:</p>
  
  <pre><code>deltax=0.5
  alpha=c(1, 1.5, 2, 3, 4)
  H=seq(-7, 7, deltax)
  sdH=seq(1, 15, 0.02)
  optimumgvariance= matrix(0,ncol = length(alpha), nrow = length(sdH))
  
  &gt; tail(optimumgvariance)
         [,1] [,2] [,3] [,4] [,5]
  [696,]  2.4    3  3.4  4.2  4.8
  [697,]  2.4    3  3.4  4.2  4.8
  [698,]  2.4    3  3.4  4.2  4.8
  [699,]  2.4    3  3.4  4.2  4.8
  [700,]  2.4    3  3.4  4.2  4.8
  [701,]  2.4    3  3.4  4.2  4.8
  </code></pre>
  
  <p>Does anybody have any idea, what causes this? I actually wanted to go up to <code>sdH= seq(1,30, 0.02)</code> but I can't trust my results when I don't know why this happens.</p>