R中的回归(对数正态):为特定的y值(结果)寻找x值(预测值)

Question

这是一个类似于R中回归(逻辑(中发布的问题：为特定y值(结果(寻找x值(预测值(。我正试图找到已知y值的x值(在这种情况下为0.000001(，该值是通过使用遗传算法拟合对数正态曲线来获得的，该曲线拟合到与母树相距一定距离的树苗密度。该算法给出了最佳拟合对数正态曲线的a和b参数。

通过使用以下代码使用uniroot，我已经获得了其他曲线(如负指数(y=0.00001的x值(适用于这些曲线(：

##calculate x value at y=0.000001 (predicted near-maximum recruitment distance)
aparam=a
bparam=b
testfn <- function (y, aparam, bparam) {
## find value of x that satisfies y = a + bx
fn <- function(x) (a * exp(-b * x)) - y
uniroot(fn, lower=0, upper= 100000000)$root
}
testfn(0.000001)

不幸的是，使用对数正规公式的相同代码不起作用。我试着通过将下边界设置为零以上来使用uniroot。但是得到一个错误代码：

Error in uniroot(fn, lower = 1e-16, upper = 1e+18) : 
f() values at end points not of opposite sign

我的代码和数据(在代码下面给出(是：

file="TR maire 1mbin.txt"
xydata <- read.table(file,header=TRUE,col.names=c('x','y'))

####assign best parameter values
a = 1.35577
b = 0.8941521

#####Plot model against data
par(mar=c(5,5,2,2))
xvals=seq(1,max(xydata$x),1)
plot(jitter(xydata$x), jitter(xydata$y),pch=1,xlab="distance from NCA (m)",
ylab=quote(recruit ~ density ~ (individuals ~ m^{2~~~ -1})))
col2="light grey"

plotmodel <-  a* exp(-(b) * xvals)
lines(xvals,plotmodel,col=col2)


####ATTEMPT 1
##calculate x value at y=0.000001 (predicted near-maximum recruitment distance)
aparam=a
bparam=b
testfn <- function (y, aparam, bparam) {
fn <- function(x) ((exp(-(((log(x/b)) * (log(x/b)))/(2*a*a))))/(a * x * sqrt(2*pi))) - y
uniroot(fn, lower=0.0000000000000001, upper= 1000000000000000000)$root
}
testfn(0.000001)

数据为：

xydata
1    1 0.318309886
2    2 0.106103295
3    2 0.106103295
4    2 0.106103295
5    3 0.063661977
6    4 0.045472841
7    5 0.035367765
8    5 0.035367765
9    7 0.048970752
10   8 0.021220659
11   8 0.021220659
12   8 0.042441318
13   9 0.018724111
14  10 0.016753152
15  10 0.016753152
16  12 0.013839560
17  13 0.025464791
18  16 0.010268061
19  17 0.009645754
20  24 0.013545102
21  25 0.032480601
22  26 0.043689592
23  27 0.006005847
24  28 0.011574905
25  31 0.062618338
26  32 0.005052538
27  42 0.003835059
28  42 0.003835059
29  44 0.003658734
30  46 0.003497911
31  48 0.006701261
32  50 0.003215251
33  50 0.006430503
34  51 0.006303166
35  58 0.002767912
36  79 0.002027452
37 129 0.003715680
38 131 0.001219578
39 132 0.001210304
40 133 0.001201169
41 144 0.001109094
42 181 0.000881745
43 279 0.001142944
44 326 0.000488955

或者还有其他方法可以解决这个问题吗？我是一名生态学家，有时R就是没有意义！

Answer 1

我的r代码中似乎有一些错误，但主要问题是我的下限太低，对数正态曲线没有扩展到那个值(我的解释(。对我有效的解决方案是：

### define the formula parameter values
a = 1.35577
b = 0.8941521

### define your formula (in this instance a log normal) in the {}
fn <- function(x,a,b,y) { ((exp(-(((log(x/b)) * (log(x/b)))/(2*a*a))))/(a * x * sqrt(2*pi))) - y}

###then use uniroot()$root calling the known parameter values and defining the value of y that is of interest (in this case 0.000001)
uniroot(fn,c(1,200000),a=a,b=b,y=0.000001)$root