CODEWITHSUNDEEP
rstatistics
MYaseen208
MYaseen208

Reputation: 23938

Fitting Redlich–Kister equation in R

One of my friends asked for help to fit Redlich–Kister equation in R.

I took first three columns of Table 3 from this article (also shown in Figure below) and used the following R code. I'm getting different results than from the article output.

Data <-
structure(list(X2 = c(0, 0.0048, 0.0086, 0.0131, 0.0188, 0.0261,
0.0411, 0.0561, 0.0768, 0.0935, 0.1425, 0.1993, 0.2799, 0.3524,
0.4103, 0.4553, 0.5505, 0.6562, 0.7295, 0.8186, 0.8825, 0.9584,
1), d5 = c(0.99996, 0.9998, 0.99972, 0.99972, 0.9998, 1.00021,
1.00129, 1.00278, 1.00486, 1.00652, 1.00968, 1.00943, 1.00453,
0.99806, 0.99251, 0.98837, 0.97969, 0.97136, 0.96601, 0.96058,
0.95709, 0.95345, 0.95155), V5 = c(NA, -0.015, -0.028, -0.046,
-0.069, -0.105, -0.184, -0.274, -0.403, -0.51, -0.796, -1.037,
-1.224, -1.274, -1.257, -1.224, -1.081, -0.874, -0.691, -0.473,
-0.308, -0.114, NA)), .Names = c("X2", "d5", "V5"), class = "data.frame", row.names = c(NA,
-23L))

fm1 <-
  lm(
  formula=
     V5 ~ -1 +
          I(X2*(1-X2)) +
          I(X2*(1-X2)*(1-2*X2)) +
          I(X2*(1-X2)*(1-2*X2)^2) +
          I(X2*(1-X2)*(1-2*X2)^3) +
          I(X2*(1-X2)*(1-2*X2)^4) +
          I(X2*(1-X2)*(1-2*X2)^5)
   , data = Data)

summary(fm1)$coef

Output

                                   Estimate Std. Error     t value     Pr(>|t|)
I(X2 * (1 - X2))                  -4.636138 0.01554354 -298.267899 1.017962e-29
I(X2 * (1 - X2) * (1 - 2 * X2))   -2.729863 0.07286814  -37.463052 3.095331e-16
I(X2 * (1 - X2) * (1 - 2 * X2)^2) -1.695933 0.15189170  -11.165411 1.150270e-08
I(X2 * (1 - X2) * (1 - 2 * X2)^3) -1.658202 0.39261960   -4.223432 7.371920e-04
I(X2 * (1 - X2) * (1 - 2 * X2)^4)  3.010149 0.23867895   12.611709 2.184653e-09
I(X2 * (1 - X2) * (1 - 2 * X2)^5)  4.492079 0.45676290    9.834596 6.220233e-08

enter image description here enter image description here enter image description here enter image description here

Upvotes: 0

Views: 1710

Answers (1)

Roland
Roland

Reputation: 132969

I couldn't find your data in the table, so I typed it in myself:

Data <- data.frame(X2=c(0, 0.0059,0.0111,0.0209,0.0334,0.0462,0.06,0.0771,0.1145,0.1621,
                        0.2259,0.316,.378,.4361,.5214,.6594,.7835,.9103,.9545,1),
                   V5=c(NA,-.041,-.081,-.165,-.284,-.415,-.562,-.734,-1.055,-1.365,-1.618,-1.764,
                        -1.777,-1.726,-1.582,-1.232,-.882,-.414,-.217,NA))

I then had a closer look at the results table and noticed that they fit a polynomial of the seventh degree (for DMEA):

Data$V5x <- Data$V5/(Data$X2)/(1-Data$X2)
Data$X2a <- 1-2*Data$X2
fm1 <- lm(V5x ~ poly(X2a, 7, raw=TRUE), data = Data)

summary(fm1)$coef

#                            Estimate Std. Error     t value     Pr(>|t|)
#(Intercept)                -6.512876 0.03473349 -187.510013 4.573934e-19
#poly(X2a, 7, raw = TRUE)1  -4.066934 0.16822455  -24.175629 3.338491e-10
#poly(X2a, 7, raw = TRUE)2  -1.033662 0.43270844   -2.388819 3.803527e-02
#poly(X2a, 7, raw = TRUE)3   3.495208 1.15371167    3.029533 1.268724e-02
#poly(X2a, 7, raw = TRUE)4  -7.186365 1.24383354   -5.777594 1.783191e-04
#poly(X2a, 7, raw = TRUE)5 -10.215220 2.37360372   -4.303676 1.552439e-03
#poly(X2a, 7, raw = TRUE)6   8.982591 0.96372896    9.320661 3.017317e-06
#poly(X2a, 7, raw = TRUE)7  10.151813 1.50463933    6.747007 5.060875e-05

This is similar to their result (for DMEA, 5). The difference might be due to a typo in the data I used, rounding errors or them using a numeric optimizer.

Upvotes: 1

Related Questions