使用 nls 再现 PROC NLIN 输出-6ren

使用 nls 再现 PROC NLIN 输出

转载作者：行者123 更新时间：2023-12-04 12:16:37

我想在 R 中实现以下 SAS 代码。

$$N_e = N_o{1-exp[\frac{(d+bN_o)(T_h N_e - T)}{(1+c N_o)}]}$$

其中 $b>0$、$c\geq 0$、$T_h>0$ 和 $T=72$。

SAS 代码是

 DATA NOTONECT;
  INPUT N0  REP FATE    NE  N02 N03 PROPEAT;
  DATALINES;
5   1   0   0   25  125 0
5   2   0   1   25  125 0.2
5   3   0   1   25  125 0.2
5   4   0   2   25  125 0.4
5   5   0   2   25  125 0.4
5   6   0   2   25  125 0.4
5   7   0   2   25  125 0.4
5   8   0   3   25  125 0.6
7   1   0   0   49  343 0
7   2   0   0   49  343 0
7   3   0   1   49  343 0.14286
7   4   0   1   49  343 0.14286
7   5   0   2   49  343 0.28571
7   6   0   2   49  343 0.28571
7   7   0   2   49  343 0.28571
7   8   0   3   49  343 0.42857
10  1   0   1   100 1000    0.1
10  2   0   1   100 1000    0.1
10  3   0   2   100 1000    0.2
10  4   0   2   100 1000    0.2
10  5   0   3   100 1000    0.3
10  6   0   3   100 1000    0.3
10  7   0   3   100 1000    0.3
10  8   0   4   100 1000    0.4
10  9   0   7   100 1000    0.7
15  1   0   1   225 3375    0.06667
15  2   0   1   225 3375    0.06667
15  3   0   3   225 3375    0.2
15  4   0   3   225 3375    0.2
15  5   0   4   225 3375    0.26667
15  6   0   5   225 3375    0.33333
15  7   0   5   225 3375    0.33333
15  8   0   5   225 3375    0.33333
20  1   0   3   400 8000    0.15
20  2   0   4   400 8000    0.2
20  3   0   7   400 8000    0.35
20  4   0   7   400 8000    0.35
20  5   0   8   400 8000    0.4
20  6   0   8   400 8000    0.4
20  7   0   9   400 8000    0.45
20  8   0   11  400 8000    0.55
25  1   0   4   625 15625   0.16
25  2   0   5   625 15625   0.2
25  3   0   6   625 15625   0.24
25  4   0   7   625 15625   0.28
25  5   0   9   625 15625   0.36
25  6   0   9   625 15625   0.36
25  7   0   13  625 15625   0.52
25  8   0   14  625 15625   0.56
30  1   0   5   900 27000   0.16667
30  2   0   8   900 27000   0.26667
30  3   0   10  900 27000   0.33333
30  4   0   11  900 27000   0.36667
30  5   0   11  900 27000   0.36667
30  6   0   12  900 27000   0.4
30  7   0   14  900 27000   0.46667
30  8   0   20  900 27000   0.66667
45  1   0   4   2025    91125   0.08889
45  2   0   7   2025    91125   0.15556
45  3   0   8   2025    91125   0.17778
45  4   0   10  2025    91125   0.22222
45  5   0   11  2025    91125   0.24444
45  6   0   14  2025    91125   0.31111
45  7   0   15  2025    91125   0.33333
45  8   0   19  2025    91125   0.42222
60  1   0   9   3600    216000  0.15
60  2   0   14  3600    216000  0.23333
60  3   0   14  3600    216000  0.23333
60  4   0   16  3600    216000  0.26667
60  5   0   18  3600    216000  0.3
60  6   0   21  3600    216000  0.35
60  7   0   24  3600    216000  0.4
60  8   0   26  3600    216000  0.43333
80  1   0   7   6400    512000  0.0875
80  2   0   11  6400    512000  0.1375
80  3   0   12  6400    512000  0.15
80  4   0   15  6400    512000  0.1875
80  5   0   17  6400    512000  0.2125
80  6   0   12  6400    512000  0.15
80  7   0   21  6400    512000  0.2625
80  8   0   23  6400    512000  0.2875
100 1   0   7   10000   1000000 0.07
100 2   0   8   10000   1000000 0.08
100 3   0   10  10000   1000000 0.1
100 4   0   11  10000   1000000 0.11
100 5   0   15  10000   1000000 0.15
100 6   0   24  10000   1000000 0.24
100 7   0   26  10000   1000000 0.26
100 8   0   33  10000   1000000 0.33
;

PROC NLIN DATA=NOTONEC2;
 PARMS BHAT= 0.001 0.01 0.1 /* initial parameter estimates */
 CHAT= 0.001 0.01 0.1 
 DHAT= 0 THHAT=3.0; 
 BOUNDS BHAT>0,CHAT>=0,
 THHAT>0; /* parameter bounds */
 T=72; /* experimental period in H */
 X=NE; /* initial predicted value */ 
 A=(DHAT+BHAT*N0)/(1+CHAT*N0); /* expression for A */
 /* define the implicit function */
 C1=EXP(-A*T); /* components of the implicit function */
 C2=A*THHAT;
 H=N0*C1*EXP(C2*X)+X-N0; /* the implicit function */
 ITER=0; /* iterations for Newton's method */
 /* Newton's method employed to find predicted number eaten */
 DO WHILE(ABS(H)>0.0001 AND ITER<50); /* stop criteria for Newton's method */
      X=X-H/(N0*C1*C2*EXP(C2*X)+1); /* new predicted value */
      H=N0*C1*EXP(C2*X)+X-N0; /* new value of implicit function */
      ITER=ITER+1; /* iteration counter */
 END;
MODEL NE=X; /* model for nonlinear least squares */
Run;

SAS 输出是

                                       The NLIN Procedure
                                  Dependent Variable NE

                                       Grid Search
                                                                  Sum of
                    BHAT        CHAT        DHAT       THHAT     Squares

                 0.00100     0.00100           0      3.0000      3714.9
                  0.0100     0.00100           0      3.0000     47602.6
                  0.1000     0.00100           0      3.0000      100896
                 0.00100      0.0100           0      3.0000      3118.3
                  0.0100      0.0100           0      3.0000     29395.6
                  0.1000      0.0100           0      3.0000     97841.3
                 0.00100      0.1000           0      3.0000      1804.1
                  0.0100      0.1000           0      3.0000      6625.5
                  0.1000      0.1000           0      3.0000     53940.3



                                   The NLIN Procedure
                                  Dependent Variable NE
                                  Method: Gauss-Newton

                                     Iterative Phase
                                                                       Sum of
             Iter        BHAT        CHAT        DHAT       THHAT     Squares

                0     0.00100      0.1000           0      3.0000      1804.1
                1    0.000396      0.0176     0.00190      3.2240      1774.5
                2    0.000327     0.00578     0.00206      3.5044      1762.4
                3    0.000307     0.00193     0.00208      3.6507      1754.8
                4    0.000302    0.000793     0.00208      3.7006      1752.3
                5    0.000299    0.000265     0.00208      3.7246      1751.1
                6    0.000299    0.000123     0.00208      3.7312      1750.7
                7    0.000298    0.000039     0.00208      3.7351      1750.6
                8    0.000298    0.000013     0.00208      3.7363      1750.5
                9    0.000298    5.224E-6     0.00208      3.7367      1750.5
               10    0.000298    2.547E-7     0.00208      3.7369      1750.5
               11    0.000298    6.092E-8     0.00208      3.7369      1750.5
               12    0.000298    3.59E-10     0.00208      3.7369      1750.5
               13    0.000298    6.37E-11     0.00208      3.7369      1750.5
               14    0.000298    -287E-13     0.00208      3.7369      1750.5
               15    0.000298    -864E-13     0.00208      3.7369      1750.5
               16    0.000298    -955E-13     0.00208      3.7369      1750.5
               17    0.000298    -983E-13     0.00208      3.7369      1750.5
               18    0.000298    -992E-13     0.00208      3.7369      1750.5
               19    0.000298    -997E-13     0.00208      3.7369      1750.5
               20    0.000298    -999E-13     0.00208      3.7369      1750.5
               21    0.000298      -1E-10     0.00208      3.7369      1750.5
               22    0.000298      -1E-10     0.00208      3.7369      1750.5
               23    0.000298      -1E-10     0.00208      3.7369      1750.5
               24    0.000298      -1E-10     0.00208      3.7369      1750.5
               25    0.000298      -1E-10     0.00208      3.7369      1750.5
               26    0.000298      -1E-10     0.00208      3.7369      1750.5
               27    0.000298      -1E-10     0.00208      3.7369      1750.5
               28    0.000407           0    0.000884      4.0885      1733.8
               29    0.000416           0    0.000825      4.0645      1733.2
               30    0.000415           0    0.000832      4.0640      1733.2
               31    0.000416           0    0.000832      4.0640      1733.2


            NOTE: Convergence criterion met.


                                   Estimation Summary

                          Method                   Gauss-Newton
                          Iterations                         31
                          Subiterations                     108
                          Average Subiterations        3.483871
                          R                             1.51E-6

                                   The NLIN Procedure

                                   Estimation Summary

                          PPC(DHAT)                    0.000057
                          RPC(DHAT)                    0.000663
                          Object                       2.84E-10
                          Objective                    1733.203
                          Observations Read                  89
                          Observations Used                  89
                          Observations Missing                0


                  NOTE: An intercept was not specified for this model.

                                          Sum of        Mean               Approx
        Source                    DF     Squares      Square    F Value    Pr > F

        Model                      3      9526.8      3175.6     157.57    <.0001
        Error                     86      1733.2     20.1535
        Uncorrected Total         89     11260.0


                                   Approx       Approximate 95%
     Parameter      Estimate    Std Error      Confidence Limits     Label

     BHAT           0.000416     0.000228    -0.00004    0.000869
     CHAT                  0            0           0           0
     DHAT           0.000832      0.00382    -0.00676     0.00842
     THHAT            4.0640       0.3575      3.3534      4.7747
     Bound0            107.9        315.6      -509.1       725.0    0 <= CHAT


                             Approximate Correlation Matrix
                           BHAT            CHAT            DHAT           THHAT

          BHAT        1.0000000               .      -0.8928900       0.7190236
          CHAT         .                      .        .               .
          DHAT       -0.8928900               .       1.0000000      -0.5474038
          THHAT       0.7190236               .      -0.5474038       1.0000000

我不知道如何在 nls 中使用牛顿法.如果有人帮我解决这个问题，我将不胜感激。提前致谢。

最佳答案

我没有和你完全一样，但我认为答案应该没问题。

设置启动条件

s0 <- list(b=0.1,c=0.1,d=0,Th=3)
X <- read.table("rogersdat.txt",header=TRUE)

计算预期食用数量的函数:

predfun <- function(b,c,d,Th,Te,N0,debug=FALSE) {
  a <- (d+b*N0)/(1+c*N0)
  r <- N0 - (1/(a*Th))*lambertW(a*Th*N0*exp(a*(Th*N0-Te)))
  if (debug) cat(mean(a),b,c,d,Th,mean(r),"\n")
  r
}

绘制数据(只是为了确保):

with(X,plot(Ne~N0))
## check starting value
lines(1:100,with(c(s0,X),predfun(b,c,d,Th,Te=72,N0=1:100)))

合身:

library(emdbook)
n1 <- nls(Ne~predfun(b,c,d,Th,Te=72,N0),data=X,
    lower=c(1e-6,1e-6,-Inf,1e-6),algorithm="port",
    start=list(b=0.1,c=0.1,d=0.002,Th=3))
summary(n1)

Formula: Ne ~ predfun(b, c, d, Th, Te = 72, N0)

Parameters:
    Estimate Std. Error t value Pr(>|t|)   
b  0.0004155  0.0008745   0.475  0.63591   
c  0.0000010  0.0657237   0.000  0.99999   
d  0.0008318  0.0067374   0.123  0.90203   
Th 4.0639937  1.4665102   2.771  0.00686 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 4.516 on 85 degrees of freedom

Algorithm "port", convergence message: relative convergence (4) 

confint.default(n1)

          2.5 %      97.5 %
b  -0.001298523 0.002129547
c  -0.128815083 0.128817083
d  -0.012373153 0.014036799
Th  1.189686504 6.938300956

估计匹配得很好，置信区间不匹配(无论如何，这些类型的置信区间在边界上有点冒险......)

我实际上建议使用二项式误差的最大似然估计会好一点。

pdat <- data.frame(N0=1:100,Ne=predict(n1,newdata=data.frame(N0=1:100)))
with(pdat,lines(N0,Ne,col=2))
library(ggplot2)
ggplot(X,aes(x=N0,y=Ne))+stat_sum(aes(size=factor(..n..)),alpha=0.5)+theme_bw()+
  geom_line(data=pdat,colour="red")

关于使用 nls 再现 PROC NLIN 输出，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/6701051/

文章推荐： r - 如何将列表分配为数据表中的观察值？

文章推荐： polymer-2.x - 聚合物应用程序路由 : intercept refresh page

文章推荐： jetty - Apache Zeppelin 0.7.3-浏览器中的HTTP错误503

java - ORA-12705 : Cannot access NLS data files or invalid environment specified Japanese NLS
在 Eclipse 中运行我的程序时出现错误ORA-12705:无法访问 NLS 数据文件或指定的环境无效。到目前为止我尝试过的事情: 我已经尝试在注册表的\HKEY_LOCAL_MACHINE\S
r - nls() 的后续步骤？
在 R 中工作，我正在尝试使用 nls() 将一些数据拟合到以下模型: y ~ c - a * exp(-b * x) 我的数据: x <- c(8, 8, 10, 10, 10, 10, 12, 1
R nls 奇异梯度
我尝试搜索有关此主题的其他线程，但没有一个修复对我有用。我有一个自然实验的结果，我想显示符合指数分布的事件连续发生的次数。我的 R shell 粘贴在下面 f x [1] 1 2 3 4
r - nls - 收敛错误
对于这个数据集: dat = structure(list(x = c(5L, 5L, 5L, 5L, 10L, 10L, 10L, 10L, 15L, 15L, 15L, 15L, 17L, 17L
java - NLS 初始化消息传递
我正在使用 NLS.initializeMessages(BUNDLE_NAME, XYZ.class) 加载常量。当“BUNDLE_NAME”和“XYZ.class”存在于同一个项目中时，它可以工
R:nls() 公式中的多项式快捷符号
使用线性模型函数 lm() 多项式公式可以包含这样的快捷符号: m <- lm(y ~ poly(x,3)) 这是一个快捷方式，可以让用户不必创建 x^2 和 x^3 变量或将它们输入到公式中，例如
r - 将 nls 应用于列表中的元素
想象这两个列表如下: seq<-seq(1,16) L1<-list(A=seq, B=seq, C=seq, D=seq) v1=c(11,15,17,19,21,22,24,25,26,27,28
r - 用分类变量拟合 nls 模型
我想拟合一个线性平台 (nls) 模型，该模型将高度描述为年龄的函数，并且我想测试区域间模型的任何参数是否存在显着差异。这是我目前所拥有的: # Create data df1 head (df1
oracle - 给定周的星期一日期，独立于 NLS
我有一个输入日期(比如 sysdate)，我想获取上一个星期一的日期。我试过 select trunc(sysdate, 'D') from dual; 但它依赖于 NLS。此外，我不想通过名称检查结
r - 初始参数估计时 nls 奇异梯度矩阵的误差
我正在尝试使用 R 中的 nls 拟合矩形双曲线。 curve.nlslrc = nls(photolrc ~ (1/(2*theta))*(AQY*PARlrc+Am-sqrt((AQY*PARlr
R : catching errors in `nls`
我正在使用 nls 拟合一些指数数据. 我正在使用的代码是: fit <- nls(y ~ expFit(times, A, tau, C), start = c(A=100, tau=-3, C=0
r - nls 中的错误 - 迭代次数超过最大值
尝试拟合我的数据时遇到问题。这是数据: x = c(1, 1.071519305, 1.148153621, 1.230268771, 1.318256739, 1.412537545, 1.513
r - nls 奇异梯度矩阵 - 在积分上限中拟合参数
我正在尝试制作 nls适合一个有点复杂的表达式，其中包括两个积分，其中两个拟合参数在其上限中。我得到了错误 "Error in nlsModel(formula, mf, start, wts) :
r - 如何设置 nls 系数的界限？
有没有办法限制 NLS 系数在 R 中的取值范围？我知道我的数据应该存在的曲线形状；然而，NLS 无法通过产生 < 1 的功率系数来产生这样的曲线。从本质上讲，我正在尝试为一组幼树茎(树苗)数据生成
r - 对正弦模型进行非线性最小二乘 (nls) 拟合
我想将以下函数拟合到我的数据中: f(x) = Offset+Amplitudesin(FrequencyT+Phase)，或根据 Wikipedia : f(x) = C+alphasin(ome
r - 如何为 nls 函数找到好的起始值？
我不明白为什么我不能为这些数据使用 nls 函数。我尝试了很多不同的起始值，但总是出现相同的错误。这是我一直在做的事情: expFct2 = function (x, a, b,c) { a*
r - nls() 错误收敛(尽管起始值很好)
我一直在研究一个曲线拟合脚本，该脚本将 3 个指数修改的高斯 (EMG) 拟合到卷积曲线。我的基函数类似于高斯分布，但包括第三个参数(前两个是 mu 和 sigma )，它确定函数的指数分量的权重。
java - 如何在不设置区域设置的情况下设置 NLS 日期格式
我们的 java 类调用 PLSQL proc，它以由 NLS_DATE_FORMAT 定义的默认格式返回日期。我们的应用程序为国际化设置了自己的语言环境，但我希望日期格式保持为“DD-MON-RR”
R 脚本 - nls 函数
谁能给我一个很好的解释，说明参数“算法”在 R 中的 nls 函数中的作用？另外，公式是如何工作的？我知道它使用了 tilda，但我真的找不到对它的实际解释。另外，起始值有多重要？我是否需要尝试多
R:使用 nls 解析时出错
我正在尝试使用 nls 将高斯峰拟合到 R 中的密度图.当我使用以下等式时: fit :2:0: unexpected end of input 1: ~ ^ 谁能指出我哪里出错了？最佳答案试

行者123

个人简介

我是一名优秀的程序员,十分优秀！

作者热门文章

滴滴打车优惠券免费领取

全站热门文章

首页

博学

6Ren·AI

商城

使用 nls 再现 PROC NLIN 输出