CODEWITHSUNDEEP
rcurve-fittingdifferential-equations
Kenneth Hoadley
Kenneth Hoadley

Reputation: 11

Fitting raw data to difference equation using R

I am trying to fit a model to a data set using the following equation.

C[n] = C[n-1]+Ces*((1-C[n-1])/(1-C[n-1]*p))

The data set looks like this: saturating curve fit. The open circles are the data, and the red line is what I am looking to achieve using some sort of curve fitting routine. The first value from the data can be utilized to represent C[[1]] and I am looking to solve for C, Ces, and p. I have included an example of the raw data and some code to model it (resulting in the above figure). I have used the nls package in R for exponential fits, but this is a bit different and I am having trouble figuring out how best to run it.

C<-c(-0.003, 0.072, 0.149, 0.239, 0.320, 0.413, 0.456, 0.526, 0.554, 
0.611, 0.632, 0.648, 0.703, 0.714, 0.725, 0.771, 0.801, 0.770, 0.787,
 0.809, 0.836, 0.871, 0.866, 0.861, 0.873, 0.891, 0.887, 0.895, 
0.919, 0.930, 0.924, 0.927, 0.939, 0.929, 0.977, 0.929, 0.948, 0.944, 
0.947, 0.990, 0.981, 0.967, 0.973, 0.970, 1.005, 1.002, 0.978, 0.997, 
1.001, 1.008, 1.008, 1.016, 0.989, 1.037, 1.001, 1.030, 1.032, 0.999, 
1.031, 1.007, 1.015, 1.018)

Ces<-0.09
p<-0.01
Cd<-C[1]
Ce<-c(Cd)
I<-2
for (I in 2:length(C)) {
    Ct<-Cd+Ces*((1-Cd)/(1-Cd*p))
    Ce<-append(Ce, Ct)
    Cd<-Ct
}
summary(lm(C ~ Ce))$coefficients
plot(C, ylim=c(-0.5,1))
par(new=TRUE)
plot(Ce, type="l", col="red", ylim=c(-0.5,1))

Upvotes: 1

Views: 147

Answers (1)

dcarlson
dcarlson

Reputation: 11046

To use nls you need to make your code into a function, but first save your original fit:

Ce.orig <- Ce   # The fit from your code
satcurve <- function(C, Ces, p) {
    Cd<-C[1]
    Ce <- C[1]
    for (I in 2:length(C)) {
        Ct<-Cd+Ces*((1-Cd)/(1-Cd*p))
        Ce<-append(Ce, Ct)
        Cd<-Ct
    }
    return(Ce)
}

Now we can try to improve on your original values:

fit <- nls(C~satcurve(C, Ces, p), start=list(Ces= -0.09, p= 0.01))
summary(fit)
# 
# Formula: C ~ satcurve(C, Ces, p)
# 
# Parameters:
#      Estimate Std. Error t value Pr(>|t|)    
# Ces  0.100160   0.003492  28.682  < 2e-16 ***
# p   -0.251757   0.087765  -2.869  0.00568 ** 
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 0.01999 on 60 degrees of freedom
# 
# Number of iterations to convergence: 12 
# Achieved convergence tolerance: 3.149e-06
C.fit <- predict(fit)

The parameters changed slightly and reversed sign. Now plot the results. The difference is very small:

plot(C)
lines(Ce.orig, col="red")
lines(C.fit, col="blue")

Plot

Upvotes: 2

Related Questions