2016-06-01 9 views
0

Я использовал для использования SigmaPlot, чтобы соответствовать модифицированную версию уравнения Никольского-Эйзенман формыR фитинга уравнение (уравнение Никольского-Эйзенман) данным

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) 

enter image description here

График и отчет, сгенерированный 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 

enter image description here

ответ

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)) 

enter image description here

Заметьте, что «журнал» является натуральный логарифм в R, если вы хотите войти в базу 10, а затем использовать log10 - это немного изменяет P2 примерно до -62 вместо -27, как указано выше ...

С y наши новые данные и с «log10» вместо «войти» в выражении формулы I получают:

> 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 
+0

спасибо - завтра будет подробный взгляд, хотя первый взгляд подсказывает, что это не совсем что мне нужно, но это дает мне возможность начать.График Sigma имеет гладкую кривую, и параметры сильно отличаются от того, что я ожидаю от этих данных. –

+0

Я также отредактировал вопрос, чтобы данные соответствовали анализу Sigmaplot, чтобы облегчить прямое сравнение. –

+0

Я подозреваю, что «журнал» - это база данных 10 здесь ... – Spacedman

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") 

screenshot

Смежные вопросы