Olá Walmes!

 

Segue abaixo a rotina que estou utilizando:

# Modelo Linear Segmentado com Resposta Platô:

cv=c(x1,x2,x3,x4,x6,x8,x9,x12,x16,x18,x24,x32,x36,x48,x64,x72,x96,x144,x192,x288,x384,x576);cv
x=c(1,2,3,4,6,8,9,12,16,18,24,32,36,48,64,72,96,144,192,288,384,576);x

plot(x,cv)

# rm(list=ls())###limpa de dados anteriores
# chutes

chute=lm(cv~1+x);chute
vetor=as.vector(chute$coefficients);vetor
chuteb0=vetor[1];chuteb0
chuteb1=vetor[2];chuteb1
chuteb2=5
grp.fit=nls(cv~(b0+b1*x)*(x<=X0)+(b0+b1*X0)*(x>X0),
start=list(b0=chuteb0,b1=chuteb1,X0=chuteb2),trace=F);grp.fit
grp.coef=coef(grp.fit);grp.coef
b0=grp.coef[1]
b1=grp.coef[2]
X0=grp.coef[3]
p=b0+X0*b1;p

#  lrp.fit<-nls(cv ~(b0 + b1*x)*(x<=x0p)
#      +(b0+b1*x0p)*(x>x0p),
#      start=list(b0=, b1=-1, x0p=5),
#      trace=F)
#  lrp.fit

summary(grp.fit)
#b0<-27.5625749
#b1<-(-0.6635418)
#X0<26.7758314
cvp<-(b0+b1*X0);cvp
xs1<-seq(0,X0,length=22);xs1
xs2<-seq(X0,576,length=22);xs2
cvsplato1<-(b0 + b1*xs1);cvsplato1
cvsplato2<-rep(b0 + b1*X0,22);cvsplato2
par(mfrow=c(1,2))
plot(x,cv, xlab='tamanho de parcela X (UEB)',ylab='CV (%)'
,xlim=c(0.5,600), ylim=c(0,40),xaxt="n",yaxt="n",pch=5)###xaxt="n",yaxt="n") apaga os eixos para logo mexer nele a vontade
title("MLSRP")
axis(1,at=seq(0,600,by=50))#para mexer nos intervalos da lenda grafica eixo x
axis(2,at=seq(0,40,by=5))#para mexer nos intervalos da lenda grafica eixo x
lines(xs1,cvsplato1,lty=1,col="red")
lines(xs2,cvsplato2,lty=1,col="red")
points(X0,P,pch=7)

#segments(26.7758314,0,26.7758314,9.795692,lty=3)##XO,0,XO,X0(cvp)
#text(110,7.0,"(26,77;9,78)PL")

sqr<-deviance(grp.fit);sqr# soma de quadrados do resíduos
sqrc<-deviance(lm(cv~1));sqrc#é a soma de qudrados corrigida para a média
r2<-1-sqr/sqrc;r2###0.82
r2a<-1-((1-r2)*21/20);r2a#0.81
AIC(grp.fit)#119.56

###################

summary(grp.fit)
### Intervalos de confiança
#X0
alpha <- 0.05

IC.b0 <- c(coef(summary(grp.fit))[1,1] -
    qt(1 - 0.5*alpha,summary(grp.fit)$df[2])*coef(summary(grp.fit))[1,2],
    coef(summary(grp.fit))[1,1] +
    qt(1 - 0.5*alpha,summary(grp.fit)$df[2])*coef(summary(grp.fit))[1,2]
    )

IC.b1 <- c(coef(summary(grp.fit))[2,1] -
    qt(1 - 0.5*alpha,summary(grp.fit)$df[2])*coef(summary(grp.fit))[2,2],
    coef(summary(grp.fit))[2,1] +
    qt(1 - 0.5*alpha,summary(grp.fit)$df[2])*coef(summary(grp.fit))[2,2]
    )

IC.X0 <- c(coef(summary(grp.fit))[3,1] -
    qt(1 - 0.5*alpha,summary(grp.fit)$df[2])*coef(summary(grp.fit))[3,2],
    coef(summary(grp.fit))[3,1] +
    qt(1 - 0.5*alpha,summary(grp.fit)$df[2])*coef(summary(grp.fit))[3,2]
    )

IC <- data.frame("b0" = IC.b0, "b1" = IC.b1, "X0" = IC. X0)
IC
t(IC)

# Modelo Quadrático segmentado com Resposta Platô:

y=c(x1,x2,x3,x4,x6,x8,x9,x12,x16,x18,x24,x32,x36,x48,x64,x72,x96,x144,x192,x288,x384,x576);cv
x=c(1,2,3,4,6,8,9,12,16,18,24,32,36,48,64,72,96,144,192,288,384,576);x
#plot(x,y)

# Obtendo chutes iniciais
chute=lm(y~1+x+I(x^2)) #considerando o modelo como se fosse linear (obter chutes)
vetor=as.vector(chute$coefficients)
chuteb0=vetor[1] # chute inicial para o b0
chuteb1=vetor[2] # chute inicial para o b1
chuteb2=vetor[3] # chute inicial para o b2
# Ajustando o modelo não linear
qrp.fit=nls(y~(b0+b1*x+b2*I(x^2))*(x<=-0.5*b1/b2)+(b0+I(-b1^2/(4*b2)))*(x>-0.5*b1/b2),
            start=list(b0=chuteb0, b1=chuteb1, b2=chuteb2),trace=F)
qrp.fit
qrp.coef=coef(qrp.fit) # armazanando somente os coeficientes
qrp.coef
# Encontrando a abscissa do platô
b0=qrp.coef[1]
b 1=qrp.coef[2]
b2=qrp.coef[3]
X0=-b1/(2*b2);X0
#X0=-0.5*(b1/b2) #outro jeito, igual resultado
X0
# Obtendo o Platô
P=b0-(b1^2)/(4*b2)
P
# Observando em detalhes os resultados do ajuste
summary(qrp.fit)
# Intervalo de confiança
confint.default(qrp.fit)
# valores preditos
fitted(qrp.fit)
# R2
SQE=summary(qrp.fit)$sigma^2*summary(qrp.fit)$df[2]
SQE
SQT=var(y)*(length(y)-1) # dúvida se é 1, 3 ou 4 graus de liberdade?
SQT
R2=1-SQE/SQT
R2
AIC(qrp.fit)#115.60
AIC(qrp.fit,k=log(length(y)))#BIC 119.96
 
# Fazendo o gráfico

par(mfrow=c(1,2))
plot(x,y,ylim=c(0,40),
     xlab="Tamanho de Parcela X (UEB)",ylab="CV(%)",xaxt="n",yaxt="n",pch=5)
title("MQSRP")
axis(1,at=seq(0,600,by=50))#para mexer nos intervalos da lenda grafica eixo x
axis(2,at=seq(0,40,by=5))#para mexer nos intervalos da lenda grafica eixo x
curve((b0+b1*x+b2*x^2)*(x<=X0)+
  (b0-(b1^2)/(4*b2))*(x>X0),0,max(x),add=T,col="red")
points(X0,P,pch=7)
#segments(X0,-3,X0,p,lty=1,col="red")
#segments(-3,P,X0,P,lty=1,col="red")
###(xo,o,xo,cv)
#identify(x,cv)
#text(110,07,"(44.44;9.37)")
##(~xo,~x,(xo;cv)

# Método de Máxima Curvatura Modificada:

cv=c(x1,x2,x3,x4,x6,x8,x9,x12,x16,x18,x24,x32,x36,x48,x64,x72,x96,x144,x192,x288,x384,x576);cv
x=c(1,2,3,4,6,8,9,12,16,18,24,32,36,48,64,72,96,144,192,288,384,576);x
par(mfrow=c(1,2))
#plot(x,cv)
modelo<-nls(cv ~ A/x^B,start=c(A=1, B=0.5),control=c(tol=1e-6, maxiter=100), trace=T)

#Estimativas do modelo
#A=34.0878
#B=0.2682
##Grafico

xs<-seq(-2,576, length=200);xs
cvs<-A/xs^B;cvs

plot(x,cv,xlab='tamanho de parcela X (UEB)',ylab='CV (%)', xlim=c(0,600), ylim=c(0,40),xaxt="n",pch=5)
axis(1,at=seq(0,600,by=50),title("MMCM"))
lines(xs,cvs,type="l",col="red")
x0<-((A^2*B^2*(2*B+1))/(B+2))^(1/(2+2*B)); x0
#CV do x0
cv.x0<-A/x0^B; cv.x0
#segments(4.91,0,4.91,22.24,lty=2,col="red")#linhas continua
#segments(4.91,0,4.91,22.24,ltx=2,col="red")#linhas continua

#text(150.5,22,"(4.910354;22.24552)MC")
modelo
summary(modelo)
confint.default(modelo,level=0.95)##intervalo de confiança 95%
+R2(residuals(modelo),dados$cv)
#3sqr<-deviance(modelo);sqr###deviance é em regressão (sqr)#sqrc<-deviance(lm(cv~1));sqrc
#r2<-1-sqr/sqrc;r2

 SQE <- summary(modelo)$sigma^2*summary(modelo)$df[2]
 SQE

 SQT & lt;- var(cv)*(length(cv)-1)# soma de quadrado total corrigida
 SQT
 R2 <- 1 - SQE/SQT
 R2#98.79588
 AIC(modelo)
 AIC(modelo,k=log(length(cv)))#BIC


Att.

Em 24/10/2012 22:19, Walmes Zeviani < walmeszeviani@gmail.com > escreveu:

==========================================================================
Walmes Marques Zeviani
LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)
Departamento de Estatística - Universidade Federal do Paraná
fone: (+55) 41 3361 3573
VoIP: (3361 3600) 1053 1173
e-mail: walmes@ufpr.br
skype: walmeszeviani
twitter: @walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
==========================================================================