<div>Olá Walmes!</div>
<div> </div>
<div>Segue abaixo a rotina que estou utilizando:</div>
<div><br /><strong># Modelo Linear Segmentado com Resposta Platô:</strong> <br /><br />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<br />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<br /><br />plot(x,cv)<br /><br /># rm(list=ls())###limpa de dados anteriores<br /># chutes</div>
<div><br />chute=lm(cv~1+x);chute<br />vetor=as.vector(chute$coefficients);vetor<br />chuteb0=vetor[1];chuteb0<br />chuteb1=vetor[2];chuteb1<br />chuteb2=5<br />grp.fit=nls(cv~(b0+b1*x)*(x<=X0)+(b0+b1*X0)*(x>X0),<br />start=list(b0=chuteb0,b1=chuteb1,X0=chuteb2),trace=F);grp.fit<br />grp.coef=coef(grp.fit);grp.coef<br />b0=grp.coef[1]<br />b1=grp.coef[2]<br />X0=grp.coef[3]<br />p=b0+X0*b1;p</div>
<div><br /># lrp.fit<-nls(cv ~(b0 + b1*x)*(x<=x0p)<br /># +(b0+b1*x0p)*(x>x0p),<br /># start=list(b0=, b1=-1, x0p=5),<br /># trace=F)<br /># lrp.fit</div>
<div><br />summary(grp.fit)<br />#b0<-27.5625749<br />#b1<-(-0.6635418)<br />#X0<26.7758314<br />cvp<-(b0+b1*X0);cvp<br />xs1<-seq(0,X0,length=22);xs1<br />xs2<-seq(X0,576,length=22);xs2<br />cvsplato1<-(b0 + b1*xs1);cvsplato1<br />cvsplato2<-rep(b0 + b1*X0,22);cvsplato2<br />par(mfrow=c(1,2))<br />plot(x,cv, xlab='tamanho de parcela X (UEB)',ylab='CV (%)'<br />,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<br />title("MLSRP")<br />axis(1,at=seq(0,600,by=50))#para mexer nos intervalos da lenda grafica eixo x<br />axis(2,at=seq(0,40,by=5))#para mexer nos intervalos da lenda grafica eixo x<br />lines(xs1,cvsplato1,lty=1,col="red")<br />lines(xs2,cvsplato2,lty=1,col="red")<br />points(X0,P,pch=7)<br /><br />#segments(26.7758314,0,26.7758314,9.795692,lty=3)##XO,0,XO,X0(cvp)<br />#text(110,7.0,"(26,77;9,78)PL")</div>
<div><br />sqr<-deviance(grp.fit);sqr# soma de quadrados do resíduos<br />sqrc<-deviance(lm(cv~1));sqrc#é a soma de qudrados corrigida para a média<br />r2<-1-sqr/sqrc;r2###0.82<br />r2a<-1-((1-r2)*21/20);r2a#0.81<br />AIC(grp.fit)#119.56</div>
<div><br />###################</div>
<div><br />summary(grp.fit)<br />### Intervalos de confiança<br />#X0<br />alpha <- 0.05<br /><br />IC.b0 <- c(coef(summary(grp.fit))[1,1] - <br /> qt(1 - 0.5*alpha,summary(grp.fit)$df[2])*coef(summary(grp.fit))[1,2],<br /> coef(summary(grp.fit))[1,1] + <br /> qt(1 - 0.5*alpha,summary(grp.fit)$df[2])*coef(summary(grp.fit))[1,2]<br /> )<br /><br />IC.b1 <- c(coef(summary(grp.fit))[2,1] - <br /> qt(1 - 0.5*alpha,summary(grp.fit)$df[2])*coef(summary(grp.fit))[2,2],<br /> coef(summary(grp.fit))[2,1] + <br /> qt(1 - 0.5*alpha,summary(grp.fit)$df[2])*coef(summary(grp.fit))[2,2]<br /> )<br /><br />IC.X0 <- c(coef(summary(grp.fit))[3,1] - <br /> qt(1 - 0.5*alpha,summary(grp.fit)$df[2])*coef(summary(grp.fit))[3,2],<br /> coef(summary(grp.fit))[3,1] + <br /> qt(1 - 0.5*alpha,summary(grp.fit)$df[2])*coef(summary(grp.fit))[3,2]<br /> )<br /><br />IC <- data.frame("b0" = IC.b0, "b1" = IC.b1, "X0" = IC.
X0)<br />IC<br />t(IC)<br /><br /><strong># Modelo Quadrático segmentado com Resposta Platô</strong>:<br /><br />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<br />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<br />#plot(x,y)<br /><br /># Obtendo chutes iniciais<br />chute=lm(y~1+x+I(x^2)) #considerando o modelo como se fosse linear (obter chutes)<br />vetor=as.vector(chute$coefficients)<br />chuteb0=vetor[1] # chute inicial para o b0<br />chuteb1=vetor[2] # chute inicial para o b1<br />chuteb2=vetor[3] # chute inicial para o b2<br /># Ajustando o modelo não linear<br />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),<br /> start=list(b0=chuteb0, b1=chuteb1, b2=chuteb2),trace=F)<br />qrp.fit<br />qrp.coef=coef(qrp.fit) # armazanando somente os coeficientes<br />qrp.coef<br /># Encontrando a abscissa do platô<br />b0=qrp.coef[1]<br />b
1=qrp.coef[2]<br />b2=qrp.coef[3]<br />X0=-b1/(2*b2);X0<br />#X0=-0.5*(b1/b2) #outro jeito, igual resultado<br />X0<br /># Obtendo o Platô<br />P=b0-(b1^2)/(4*b2)<br />P<br /># Observando em detalhes os resultados do ajuste<br />summary(qrp.fit)<br /># Intervalo de confiança<br />confint.default(qrp.fit)<br /># valores preditos<br />fitted(qrp.fit)<br /># R2<br />SQE=summary(qrp.fit)$sigma^2*summary(qrp.fit)$df[2]<br />SQE<br />SQT=var(y)*(length(y)-1) # dúvida se é 1, 3 ou 4 graus de liberdade?<br />SQT<br />R2=1-SQE/SQT<br />R2<br />AIC(qrp.fit)#115.60<br />AIC(qrp.fit,k=log(length(y)))#BIC 119.96<br /> <br /># Fazendo o gráfico</div>
<div><br />par(mfrow=c(1,2))<br />plot(x,y,ylim=c(0,40),<br /> xlab="Tamanho de Parcela X (UEB)",ylab="CV(%)",xaxt="n",yaxt="n",pch=5)<br />title("MQSRP")<br />axis(1,at=seq(0,600,by=50))#para mexer nos intervalos da lenda grafica eixo x<br />axis(2,at=seq(0,40,by=5))#para mexer nos intervalos da lenda grafica eixo x<br />curve((b0+b1*x+b2*x^2)*(x<=X0)+<br /> (b0-(b1^2)/(4*b2))*(x>X0),0,max(x),add=T,col="red")<br />points(X0,P,pch=7)<br />#segments(X0,-3,X0,p,lty=1,col="red")<br />#segments(-3,P,X0,P,lty=1,col="red")<br />###(xo,o,xo,cv)<br />#identify(x,cv)<br />#text(110,07,"(44.44;9.37)")<br />##(~xo,~x,(xo;cv)<br /><br /><strong># Método de Máxima Curvatura Modificada:</strong></div>
<div><br />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<br />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<br />par(mfrow=c(1,2))<br />#plot(x,cv)<br />modelo<-nls(cv ~ A/x^B,start=c(A=1, B=0.5),control=c(tol=1e-6, maxiter=100), trace=T)</div>
<div>#Estimativas do modelo<br />#A=34.0878<br />#B=0.2682<br />##Grafico<br /><br />xs<-seq(-2,576, length=200);xs<br />cvs<-A/xs^B;cvs<br /><br />plot(x,cv,xlab='tamanho de parcela X (UEB)',ylab='CV (%)', xlim=c(0,600), ylim=c(0,40),xaxt="n",pch=5)<br />axis(1,at=seq(0,600,by=50),title("MMCM"))<br />lines(xs,cvs,type="l",col="red")<br />x0<-((A^2*B^2*(2*B+1))/(B+2))^(1/(2+2*B)); x0<br />#CV do x0<br />cv.x0<-A/x0^B; cv.x0<br />#segments(4.91,0,4.91,22.24,lty=2,col="red")#linhas continua<br />#segments(4.91,0,4.91,22.24,ltx=2,col="red")#linhas continua<br /><br />#text(150.5,22,"(4.910354;22.24552)MC")<br />modelo<br />summary(modelo)<br />confint.default(modelo,level=0.95)##intervalo de confiança 95%<br />+R2(residuals(modelo),dados$cv)<br />#3sqr<-deviance(modelo);sqr###deviance é em regressão (sqr)#sqrc<-deviance(lm(cv~1));sqrc<br />#r2<-1-sqr/sqrc;r2<br /><br /> SQE <- summary(modelo)$sigma^2*summary(modelo)$df[2]<br /> SQE<br /><br /> SQT &
lt;- var(cv)*(length(cv)-1)# soma de quadrado total corrigida<br /> SQT<br /> R2 <- 1 - SQE/SQT<br /> R2#98.79588<br /> AIC(modelo)<br /> AIC(modelo,k=log(length(cv)))#BIC<br /><br /><br /><em>Att.</em></div>
<div><em>André</em></div>
<div><em> </em></div>
<hr style="border-top: 1px solid #ccc;" />
<div>Em 24/10/2012 22:19, <strong>Walmes Zeviani < walmeszeviani@gmail.com ></strong> escreveu:</div>
<div><span style="font-family: trebuchet ms,sans-serif;">Bem, você não passou CMR portanto eu não sei como anda sua programação, mas o linear platô seria algo assim</span></div>
<div><span style="font-family: trebuchet ms,sans-serif;"> </span></div>
<div>
<div><span style="font-family: courier new,monospace;">ym </span></div>
<div><span style="font-family: courier new,monospace;">xm </span></div>
<div><span style="font-family: courier new,monospace;">cm </span></div>
<div><span style="font-family: courier new,monospace;">curve(ym+cm*(x-xm)^1*(x<=xm), 0, 10, main="linear platô")</span></div>
<div><span style="font-family: courier new,monospace;"> </span></div>
<div><span style="font-family: courier new,monospace;">x </span></div>
<div><span style="font-family: courier new,monospace;">y </span></div>
<div><span style="font-family: courier new,monospace;">plot(y~x)</span></div>
<div><span style="font-family: courier new,monospace;"> </span></div>
<div><span style="font-family: courier new,monospace;">n0 </span></div>
<div><span style="font-family: courier new,monospace;">summary(n0)</span></div>
<div><span style="font-family: courier new,monospace;">confint(n0)</span></div>
<div><span style="font-family: courier new,monospace;">confint.default(n0)</span></div>
<div style="font-family: 'trebuchet ms',sans-serif;"> </div>
<div style="font-family: 'trebuchet ms',sans-serif;">À disposição.</div>
<div style="font-family: 'trebuchet ms',sans-serif;">Walmes.</div>
</div>
<div><span style="font-family: trebuchet ms,sans-serif;"><br /></span><span style="font-family: trebuchet ms,sans-serif;">==========================================================================</span><br /> <span style="font-family: trebuchet ms,sans-serif;">Walmes Marques Zeviani</span><br /><span style="font-family: trebuchet ms,sans-serif;">LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)</span><br /> <span style="font-family: trebuchet ms,sans-serif;">Departamento de Estatística - Universidade Federal do Paraná</span><br /><span style="font-family: trebuchet ms,sans-serif;">fone: (+55) 41 3361 3573</span><br /> <span style="font-family: trebuchet ms,sans-serif;">VoIP: (3361 3600) 1053 1173</span><br /><span style="font-family: trebuchet ms,sans-serif;">e-mail: <a href="../../../undefined//compose?to=walmes@ufpr.br" target="_blank">walmes@ufpr.br</a><br /> skype: walmeszeviani<br /></span><span style="font-family: trebuchet ms,sans-serif;
">twitter: @walmeszeviani</span><br /><span style="font-family: trebuchet ms,sans-serif;">homepage: <a href="http://www.leg.ufpr.br/%7Ewalmes" target="_blank">http://www.leg.ufpr.br/~walmes</a></span><br /> <span style="font-family: trebuchet ms,sans-serif;">linux user number: 531218</span><br /><span style="font-family: trebuchet ms,sans-serif;">==========================================================================</span><br /> </div>