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 sdH and a trade-off function.
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?
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
This is the first set of parameters I'm using in the follwing 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)>=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
}
This produces results as I would expect:
> 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
> max(optimumgvariance[,2])
[1] 2.3
> max(optimumgvariance[,1])
[1] 1.5
However, when I only change the scope of sdH=seq(1,10,0.02) 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:
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))
> 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
> max(optimumgvariance[,1])
[1] 1.5
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 sdH= seq(1,15, 0.02) I even get values other than zero in the first column:
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))
> 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
Does anybody have any idea, what causes this? I actually wanted to go up to sdH= seq(1,30, 0.02) but I can't trust my results when I don't know why this happens.
