#UNIVERSIDAD NACIONAL DEL ALTIPLANO
#INGENIERIA ESTADISTICA E INFORMATICA
#REGRESION AVANZADA
#REGRESION INTRISICAMENTE NO LINEAL

y = c(80.574 ,84.248 ,87.264 ,87.195 ,89.076 ,89.608 ,89.868 ,
      90.101 ,92.405 ,95.854 ,100.696 ,101.06 ,401.672 ,390.724 ,
      567.534 ,635.316 ,733.054 ,759.087 ,894.206 ,990.785 ,1090.109 ,
      1080.914 ,1122.643 ,1178.351 ,1260.531 ,1273.514 ,1288.339 ,
      1327.543 ,1353.863 ,1414.509 ,1425.208 ,1421.384 ,1442.962 ,
      1464.35 ,1468.705 ,1447.894 ,1457.628)

x = c(-3.067 ,-2.981 ,-2.921 ,-2.912 ,-2.84 ,-2.797 ,
      -2.702 ,-2.699 ,-2.633 ,-2.481 ,-2.363 ,-2.322 ,-1.501 ,-1.46 ,
      -1.274 ,-1.212 ,-1.1 ,-1.046 ,-0.915 ,-0.714 ,-0.566 ,-0.545 ,
      -0.4 ,-0.309 ,-0.109 ,-0.103 ,0.01 ,0.119 ,0.377 ,0.79 ,0.963 ,
      1.006 ,1.115 ,1.572 ,1.841 ,2.047 ,2.2)

#Representation de los datos
plot(y ~ x,xlab="Log de Densidad",
     ylab="Mobilidad de los electrones")

foo = function(x,b1,b2,b3,b4,b5,b6,b7){
  (b1 + b2*x + b3*x^2 + b4*x^3)/
    (1 + b5*x + b6*x^2 + b7*x^3)}

plot(y ~ x,xlab="Log de Densidad",
     ylab="Mobilidad de los electrones",main="Prueba 1")
curve(foo(x,1200,100,100,1,0.1,1,0.1),add=T,col="red")

plot(y ~ x,xlab="Log de Densidad",
     ylab="Mobilidad de los electrones",main="Prueba 2")
curve(foo(x,1200,100,100,1,-0.1,0.1,-0.1),add=T,col="blue")

m1start=list(b1=1200,b2=100,b3=100,b4=1,b5=-0.1,b6=0.1,b7=-0.1)

library(nls2)
## Warning: package 'nls2' was built under R version 4.0.5
## Loading required package: proto
## Warning: package 'proto' was built under R version 4.0.5
m1<- nls2(y ~ foo(x,b1,b2,b3,b4,b5,b6,b7),
          start=c(b1=1200,b2=100,b3=100,b4=1,b5=-0.1,b6=0.1,b7=-0.1),
          control = nls.control(warnOnly = TRUE))

B1=1000
B2=1000
B3=400
B4=40
B5=0.7
B6=0.3
B7=0.03

m.nls = nls(y ~ foo(x,b1,b2,b3,b4,b5,b6,b7),
            start = c(b1 = 1000, b2 = 1000, b3 = 400, b4 = 40,
                      b5 = 0.7, b6 = 0.3, b7 = 0.03),trace=T)
## 4528125 :  1e+03 1e+03 4e+02 4e+01 7e-01 3e-01 3e-02
## 2917852 :  1.017774e+03 1.154289e+03 5.125936e+02 7.798108e+01 8.450551e-01 3.683709e-01 6.079666e-02
## 2211486 :  1.051789e+03 1.190915e+03 5.099180e+02 7.455408e+01 8.587374e-01 3.623130e-01 5.733611e-02
## 1224759 :  1111.3074206 1269.7029496  520.9968386   71.9790313    0.8947483    0.3641991    0.0536446
## 382250.7 :  1.200294e+03 1.401290e+03 5.607441e+02 7.378599e+01 9.539234e-01 3.860543e-01 5.197037e-02
## 44573.78 :  1.288335e+03 1.490599e+03 5.806624e+02 7.350805e+01 9.668273e-01 3.970028e-01 4.801509e-02
## 7236.668 :  1.288056e+03 1.501817e+03 5.920922e+02 7.746637e+01 9.742614e-01 4.027477e-01 5.168883e-02
## 5672.79 :  1.288112e+03 1.489220e+03 5.822912e+02 7.528677e+01 9.653672e-01 3.978014e-01 4.920344e-02
## 5644.677 :  1.288113e+03 1.493501e+03 5.849135e+02 7.573829e+01 9.679714e-01 3.987010e-01 5.018898e-02
## 5643.374 :  1.288145e+03 1.489664e+03 5.822542e+02 7.522382e+01 9.653567e-01 3.975440e-01 4.943436e-02
## 5643.05 :  1.288133e+03 1.492128e+03 5.839721e+02 7.555999e+01 9.670080e-01 3.982977e-01 4.993472e-02
## 5642.85 :  1.288143e+03 1.490417e+03 5.827773e+02 7.532649e+01 9.658518e-01 3.977708e-01 4.959276e-02
## 5642.776 :  1.288137e+03 1.491544e+03 5.835633e+02 7.548015e+01 9.666098e-01 3.981163e-01 4.982007e-02
## 5642.738 :  1.288141e+03 1.490775e+03 5.830262e+02 7.537516e+01 9.660905e-01 3.978796e-01 4.966585e-02
## 5642.722 :  1.288138e+03 1.491288e+03 5.833842e+02 7.544515e+01 9.664361e-01 3.980372e-01 4.976913e-02
## 5642.714 :  1.288140e+03 1.490940e+03 5.831415e+02 7.539770e+01 9.662015e-01 3.979302e-01 4.969933e-02
## 5642.711 :  1.288139e+03 1.491174e+03 5.833042e+02 7.542951e+01 9.663587e-01 3.980019e-01 4.974622e-02
## 5642.71 :  1288.1400216 1491.0160857  583.1942944   75.4080283    0.9662525    0.3979535    0.0497146
## 5642.709 :  1.288139e+03 1.491122e+03 5.832682e+02 7.542247e+01 9.663238e-01 3.979860e-01 4.973587e-02
## 5642.709 :  1.288140e+03 1.491051e+03 5.832183e+02 7.541272e+01 9.662757e-01 3.979640e-01 4.972152e-02
## 5642.708 :  1288.1395737 1491.0986868  583.2519318   75.4192955    0.9663081    0.3979788    0.0497312
## 5642.708 :  1.288140e+03 1.491066e+03 5.832293e+02 7.541486e+01 9.662862e-01 3.979688e-01 4.972467e-02
## 5642.708 :  1.288140e+03 1.491088e+03 5.832445e+02 7.541784e+01 9.663009e-01 3.979756e-01 4.972906e-02
## 5642.708 :  1.288140e+03 1.491073e+03 5.832342e+02 7.541584e+01 9.662910e-01 3.979710e-01 4.972611e-02
## 5642.708 :  1288.1396581 1491.0832297  583.2411433   75.4171866    0.9662977    0.3979741    0.0497281
## 5642.708 :  1.288140e+03 1.491077e+03 5.832365e+02 7.541628e+01 9.662932e-01 3.979720e-01 4.972676e-02
## 5642.708 :  1.288140e+03 1.491081e+03 5.832396e+02 7.541689e+01 9.662963e-01 3.979734e-01 4.972766e-02
## 5642.708 :  1.288140e+03 1.491078e+03 5.832375e+02 7.541648e+01 9.662942e-01 3.979725e-01 4.972705e-02
## 5642.708 :  1.288140e+03 1.491080e+03 5.832389e+02 7.541675e+01 9.662956e-01 3.979731e-01 4.972746e-02
deviance(m.nls)
## [1] 5642.708
#Variance of just

m.lm <- lm(y~x)
anova(m.nls,m.lm)
## Analysis of Variance Table
## 
## Model 1: y ~ foo(x, b1, b2, b3, b4, b5, b6, b7)
## Model 2: y ~ x
##   Res.Df Res.Sum Sq Df  Sum Sq F value    Pr(>F)    
## 1     30       5643                                 
## 2     35     905627 -5 -899984  956.97 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Independence de los residues
plot(fitted(m.nls),residuals(m.nls),
     xlab="Valores ajustados",ylab="Residuos")
abline(a=0,b=0,col="blue")

#Test de normalized de los residues
qqnorm(residuals(m.nls))
qqline(residuals(m.nls))

shapiro.test(residuals(m.nls))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(m.nls)
## W = 0.93462, p-value = 0.0312
plot(fitted(m.nls),residuals(m.nls),
     xlab="Valores ajustados",ylab="Residuos")
abline(a=0,b=0,col="blue")