我正在尝试计算几何平均值并为一组参数找到它的最大可能值。我正在对两个分布进行卷积（我认为，不确定英语中的术语），其中一个是我的方程的概率密度函数，并且想要找到几何平均值最高的参数集。我将尝试注释我的代码，以便更清楚我在做什么。我这样做是为了确定某个物种在某些情况下的最佳生态位宽度。适应度取决于物种性状与环境标准差之间的匹配以及sdH权衡函数。
但是，如果我更改其中一个参数，它会更改完整的输出，即使我不希望它这样做。也许我忽略了一些错误，或者这确实是可以预料的，我的期望是错误的，但我无能为力。也许有人知道如何做到这一点？
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
这是我在以下代码中使用的第一组参数：
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
}
这会产生我期望的结果：
> 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
但是，当我只更改范围sdH=seq(1,10,0.02)并保留其他所有内容时，结果会发生变化。我突然在第二列中有值，我之前有预期的零：
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
我希望第 2 列中的值保持不变，但不知何故他们没有。这更进一步，当我把sdH= seq(1,15, 0.02)我什至在第一列中得到非零值时：
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
有谁知道，是什么原因造成的？我实际上想上升，sdH= seq(1,30, 0.02)但当我不知道为什么会发生这种情况时，我不能相信我的结果。