0

我曾经使用SigmaPlot来拟合如下形式的 Nikolsky-Eisenman 方程的修改版本

y = P1 + P2 * log(10^(-x) + P3)

使用全局曲线拟合函数。参数的详细信息可以在下面的 Sigmaplot 报告中找到。我现在想在 R 中执行此操作。

一些数据：

pNO3 <- c(1.1203, 2.0410, 3.0155, 4.0048, 4.3045, 5.0, 6.0)
mV <- c(45.2, 100.9, 160.9, 215.7, 231.5, 244.5, 257.4)
data <- data.frame(pNO3, mV)
plot(data$pNO3, data$mV)

在此处输入图像描述

sigmaplot 为上述数据生成的图表和报告如下所示。谁能指出我如何在 R 中生成类似内容的正确方向？

NonLinear Regression - Global Curve Fitting     Wednesday, May 01, 2013, 13:04:55
Data Source: Data 1 in Notebook1
Equation: User-Defined, Nicolsky Eisenman
f=P1+P2*log(10^(-x)+P3)
Data Set Specifications:
Data Set    Independent Variable    Dependent Variable
1           Column 3        Column 7
Global Parameters:
A Global Parameter is shared across all data sets.
Global Goodness of Fit:
R   Rsqr    Adj Rsqr    Standard Error of Estimate
0.9997  0.9994  0.9991      2.4421  
Analysis of Variance: 
Analysis of Variance: 
    DF  SS  MS  
Regression    3 264242.5551 88080.8517
Residual    4   23.8549 5.9637  
Total   7   264266.4100 37752.3443  
Corrected for the mean of the observations:
    DF  SS  MS  F   P   
Regression  2   38844.3822  19422.1911  3256.7192   <0.0001 
Residual    4   23.8549 5.9637  
Total   6   38868.2371  6478.0395   
Statistical Tests:
Normality Test (Shapiro-Wilk)           Passed  (P = 0.4003)
W Statistic= 0.9106 Significance Level = 0.0500
Constant Variance Test      Passed  (P = 0.1209)
Number of Observations =  7 
Rsqr  = 0.9994 
Residual Sum of Squares = 23.8549 
Parameter Estimates:
    Coefficient Std. Error  t   P
P1  -24.3265      3.3330    -7.2987 0.0019
P2  -61.7088      1.2861    -47.9796    <0.0001
P3  2.8351E-005   4.6040E-006   6.1579  0.0035
Fit Equation Description:
[Variables]
f0_x = col(3)
f0_y = col(7)
[Parameters]
f0_P1 = 0 ' {{previous: -24.3265}}
f0_P2 = -5 ' {{previous: -61.7088}}
f0_P3 = 0 ' {{previous: 2.8351e-005}}
[Equation]
f0 = f0_P1+f0_P2*log(10^(-f0_x)+f0_P3)
fit f0 to f0_y
[Constraints]
[Options]
tolerance=0.000100
stepsize=100
iterations=100
Number of Iterations Performed = 4

在此处输入图像描述

4

2 回答 2

1

假设您的拟合标准是最小化平方和误差，您可以使用nls，但您确实需要一个公平的起始值。我不知道你的参数什么是合理的，所以我花了一段时间直到我从你的 sigmaplot 示例中复制了参数，我猜这对于可能与这个数据集相似的数据集是合理的。无论如何，如果您知道参数的含义，那么您可能会猜出合理的起始值。

> start=list(P1=-24,P2=-61,P3=2.8e-5)
> m = nls(formula= mV ~ P1 + P2 * log(10^(-pNO3) + P3),data=data,start=start)
> summary(m)

Formula: mV ~ P1 + P2 * log(10^(-pNO3) + P3)

Parameters:
     Estimate Std. Error t value Pr(>|t|)    
P1 -1.420e+01  4.642e+00  -3.059    0.055 .  
P2 -2.732e+01  9.257e-01 -29.514 8.54e-05 ***
P3  8.417e-05  1.818e-05   4.630    0.019 *

pNO3您可以通过创建一组新的度量来绘制数据和平滑曲线拟合：

 plot(data$pNO3,data$mV)
 newdata = data.frame(pNO3=seq(1,6,len=100))
 lines(newdata$pNO3,predict(m, newdata=newdata))

在此处输入图像描述

请注意，“log”是 R 中的自然对数，如果您想以 10 为底，请使用log10- 这会将 P2 稍微更改为大约 -62 而不是上面的 -27 ...

使用您的新数据并在公式表达式中使用“log10”而不是“log”，我得到：

> m10 = nls(formula= mV ~ P1 + P2 * log10(10^(-pNO3) + P3),data=data,start=start)
> summary(m10)

Formula: mV ~ P1 + P2 * log10(10^(-pNO3) + P3)

Parameters:
     Estimate Std. Error t value Pr(>|t|)    
P1 -2.433e+01  3.334e+00  -7.298  0.00187 ** 
P2 -6.171e+01  1.286e+00 -47.972 1.13e-06 ***
P3  2.835e-05  4.605e-06   6.157  0.00353 **

这看起来像您的 Sigmaplot 输出：

Parameter Estimates:
    Coefficient Std. Error  t   P
P1  -24.3265      3.3330    -7.2987 0.0019
P2  -61.7088      1.2861    -47.9796    <0.0001
P3  2.8351E-005   4.6040E-006   6.1579  0.0035
于 2016-06-01T18:35:45.783 回答
0

仅使用非线性输入的参数的"plinear"算法nls需要起始值。请注意，plinear需要更改公式以表示此模型：

fo <- mV ~ cbind(log10(10^(-pNO3) + P3), 1)
fm <- nls(fo, data, start = c(P3 = 0), alg = "plinear")
summary(fm)

给予：

Formula: mV ~ cbind(log10(10^(-pNO3) + P3), 1)

Parameters:
       Estimate Std. Error t value Pr(>|t|)    
P3     2.84e-05   4.60e-06    6.16   0.0035 ** 
.lin1 -6.17e+01   1.29e+00  -47.97  1.1e-06 ***
.lin2 -2.43e+01   3.33e+00   -7.30   0.0019 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.44 on 4 degrees of freedom

Number of iterations to convergence: 7 
Achieved convergence tolerance: 2.94e-06

我们可以像这样绘制它：

plot(data)
lines(fitted(fm) ~ pNO3, data, col = "red")

截屏

于 2016-06-02T11:06:27.193 回答