There is too much here for me to parse it right now. Did you walk
through the code line by line looking at the output anc checking what it
did at each step?
If you don't get a satisfactory reply from someone else and you can't
figure it out from a line-by-line analysis of your code, please simplify
the question: Make up a toy problem that illustrates in only a few
lines what you don't understand. You are more likely to get the answer
you want if it takes seconds rather than minutes to understand your
question.
hth. spencer graves
Mei-Chu Lo wrote:
Dear R users,
I am new in R and I want to use the nls package to analyze some
experimental data. The data is in the attached file data. It is the
response Sav measured at different C0. Basically, the C0 is a
function of C1, K2, and r, and the Sav is a function of C0, C1, K2,
and r. The math equations are shown in the attached fileequations.
The parameters K2 and r are the physical properties I want to get from
the non-linear regression. The R codes I wrote is in the attached
fileRcode. Basically, I wrote two functions. First , I calculated
the C1 for different C0 at the estimated K2 and r using the binary
search method and implemented it in the function fn1. Then for the
calculated C1 and estimated K2 and r, I used the function fn2 to get
the estimated Sav for different C0. nls was used to minimize the
differences between the calculated Sav and observed Sav. When I run my
R script, I got the error message step factor reduced below minFactor
. If I changed the minFactor to zero, the nls continued but did not
converge (exceed the maxiteration). I changed the tolerance to higher
value, nls finished but the fitting is bad. From the published results,
the best fitted values for K2 and r are 0.2237 and 1.296*10^7,
respectively. I can't get these numbers using my R script.
I know there are a lot math and R experts on the R mailing-list. I
will appreciate it if anyone can tell me what is wrong in my R script or
in the methods I used to get these parameters.
Mei-Chu Lo
C0 Sav
0.6766 6.0875
1.6165 6.4249
1.8796 6.5374
2.4436 7.025
4.8872 7.5125
5.3759 7.625
5.7518 7.8499
7.218 8.0749
9.2105 8.1125
12.7067 9.4624
12.5939 10.025
16.203 11.7125
17.2932 12.0124
18.1578 12.5749
21.5413 12.6875
21.5413 13.0625
fn1-function(C0,k2,r){
m-matrix(nrow=26,ncol=length(C0))
Ct-vector(numeric)
C1-vector(numeric)
dC-vector(numeric)
j=1
repeat{
C1[j]-C0[j]/2;dC[j]-C0[j]/4
repeat{
m[1,j]=C1[j];m[2,j]=k2*C1[j]^2
m[26,j]=26*(k2/2)^25*r*C1[j]^26
i=3
repeat{
m[i,j]-i*(k2/2)^(i-1)*C1[j]^i
i-i+1
if(i=26) break}
Ct[j]-sum(m[,j])
if(abs(Ct[j]-C0[j])10^-6*C0[j]) break
if(Ct[j]C0[j]) {
C1[j]=C1[j]-dC[j];dC[j]-dC[j]/2}
else{
C1[j]=C1[j]+dC[j];dC[j]-dC[j]/2}
}
j=j+1
if (jlength(C0)) break }
C1
}
fn2-function(C1,C0,k2,r){
m-matrix(nrow=26,ncol=length(C0))
m[1,]-5.8*(1-0.018*C0)*C1;
m[2,]-2^(2/3)*5.8*(1-0.018*C0)*k2*C1^2
m[26,]-42*(1-0.019*C0)*26*(k2/2)^25*r*C1^26
j=3
repeat{
m[j,]-j^(5/3)*5.8*(1-0.018*C0)*(k2/2)^(j-1)*C1^j
j=j+1
if(j=26)break
}
p-vector(numeric)
k=1
repeat{
p[k]-sum(m[,k])/C0[k]
k=k+1
if (klength(C0))break
}
p
}
# load data
ass-read.table(data.txt,header=T)
plot(ass$C0,ass$Sav,cex=1.5, xlab=C0(mg/ml),ylab=Sav(S))
# give initial guess of k2 and r
k2=0.3;r=10^6
fit-nls(Sav~fn2(fn1(ass$C0,k2fit,rfit),ass$C0,k2fit,rfit),data=ass,start=list(k2fit=k2,rfit=r))
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
