# Ajuste de curvas por mínimos quadrados.

# Dados os vetores x e y, supomos que y_i=f(x_i)+erro_i, como erro_i um erro "pequeno".
x=seq(1930,2010,by=10)
x   # Mudança de escala
n=length(x)

y=c(0.63,0.76,0.92,1.1,1.2,1.3,1.5,1.7,1.9)

n=length(x)
plot(x,y,col="red",xlab="t",ylab="P")

yy=log(y)
plot(x,yy,col="red",ylab="ln(P)",xlab="t")


A=matrix(1,9,2);A[,1]=x;A
aa=solve(t(A)%*%A,t(A)%*%yy);aa

reta<-function(s){aa[1]*s+aa[2]}
curve(reta,x[1],x[n],ylab="L(t)",xlab="t")
points(x,yy,col="red")
ye=y; for ( i in 1:9){ye[i]=reta(x[i])}
points(x,ye,col="blue")


expreta<-function(s){exp(aa[1]*s+aa[2])}
curve(expreta,x[1],x[n],ylab="P(t)",xlab="t")
points(x,y,col="red")
ye=y; for ( i in 1:9){ye[i]=expreta(x[i])}
points(x,ye,col="blue")

#-------Linearizado

Z=1:8
xx=0*Z
for ( i in 1:8){Z[i]=(y[i+1]-y[i])/(10*y[i])
xx[i]=1
}

plot(y[1:8],Z,col="red",xlab="P_k",ylab="Z_k")

A=sum(Z)/8

points(y[1:8],A*xx,col="blue")
points(y[1:8],A*xx,col="black",'l')

yy=log(y)
plot(x,yy,col="red",ylab="ln(y)")


A=matrix(1,9,2);A[,1]=1:9;A
aa=solve(t(A)%*%A,t(A)%*%yy);aa
P0=exp(aa[2]);P0

reta<-function(s){aa[1]*s+aa[2]}
curve(reta,1,9,ylab="L(t)",ylim=c(-0.5,0.7),xlab="t")
points(1:9,yy,col="red")
ye=y; for ( i in 1:9){ye[i]=reta(i)}
points(1:9,ye,col="blue")

expreta<-function(s){exp(aa[1]*s+aa[2])}
curve(expreta,1,9,ylab="P(t)",xlab="t")
points(1:9,y,col="red")
ye=y; for ( i in 1:9){ye[i]=P0*exp(aa[1]*i)}
points(1:9,ye,col="blue")




# Ajuste no modelo 2
x=(x-1930)/10

# Suponha que y=f(x) (aproximado) e que f(s)= ace^(as)/(1+bce^(as)) e f'(s)=(a-bf(s))f(s). Descobrir aproximações para a, b e c.
# 
# Note que f(0) nos dá c, caso conhecermos a e b. 
# 
# Para determinar os valores dos parâmetros resolvemos o problema de minimização 
#
#        min g(u)
# 
# em que g(u)= sum (y_i-f(x_i))^2, sendo u=(a,b,c(a,b)).
# 
# Para fazer isso usamos o método de Euler para resolver a equaçõ diferencial 
#
#      u'(t)=-grad g(u(t))
#      u(0)=u_0

# Seguem os comando para as definições de funções necessárias.

g<-function(u){    # função g(u)
a=u[1];b=u[2]; c=u[3]
p=0
for (i in 1:n){
p=p+( y[i]-a/( b+c*exp(-a*x[i]) ) )^2
}
p
}

# Teste do código para g(u).
g(c(1,2,1)) 

gradng<-function(u){   # função gradiente de g(u).
h=0.00001
p=0*u
m=length(u)
for ( i in 1:m){v=u; v[i]=u[i]+h;p[i]=g(v)-g(u)}
p/h
}

# Teste da função gradg(u)
gradng(c(1,2,1))


#---- Início da rotin para o método de Euler.

ZeroEuler<-function(u0,t,m){
w=u0
h=t/m
for (i in 1:m){w=w-h*gradng(w)}
w
}

u0=c(1,1,1)
g(u0)
u=ZeroEuler(u0,40,400000)
#---------------------------------------------------------------------------------------

print("Soulução"); u                                # Aproximação para u que minimisa |Au-v|^2.
a=u[1]; b=u[2]; c=u[3]

print(" Resíduo da aproximação"); g(u)
print("Gradiente no ponto aproximado"); gradng(u)
print("Parâmetros a, b e c encontrados"); a;b;c
fapmv<-function(s){a/(b+c*exp(-a*(s-1930)/10))}  # Aproximação da função linearizada.

print("Gráfico para ajuste")
curve(fapmv,1930,2010,xlab="t",ylab="P(t)")                           # Teste de ajuste
points(seq(1930,2010,by=10),y,col="red")
ye=y; for ( i in 1:9){ye[i]=fapmv(x[i])}
points(seq(1930,2010,by=10),ye,col="blue")


print("População estimada em 2021");fapmv(2021)

#-------Linearizado

Z=1:8
for ( i in 1:8){Z[i]=(y[i+1]-y[i])/(10*y[i])}

plot(y[1:8],Z,col="red")

A=matrix(1,8,2);A[,2]=-y[1:8];A
aa=solve(t(A)%*%A,t(A)%*%Z);aa

reta<-function(s){aa[1]-aa[2]*s}
curve(reta,0.6,1.6,ylim=c(0,0.023),,xlab="P_k",ylab="Z_k")
points(y[1:8],Z,col="red",ylim=c(0,0.023))
ye=y[1:8]; for ( i in 1:8){ye[i]=reta(y[i])}
points(y[1:8],ye,col="blue",ylim=c(0,0.023))


P=y
for (i in 1:8){
    P[i+1]=(1+aa[1]*10-aa[2]*10*P[i])*P[i]
    }
    plot(x,P,col="blue",xlab="t_k",ylab="P_k",ylim=c(0.5,1.9))
    points(x,y,col="red")
# Aproximação em 2021
P10=(1+aa[1]*11-aa[2]*11*y[9])*y[9];P10