Gab1B-n, R - rextester

Language: Layout:

# 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=(x-1930)/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")

# 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, c e d.
# 
# 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=(a-y[1]*b)/y[1]
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))

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))

#---- 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)
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=(a-y[1]*b)/y[1]

g(u)
gradng(u)
print("Parâmetros a, b e c d 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)                           # Teste de ajuste
points(seq(1930,2010,by=10),y,col="red")

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

[ + ] Show input
Absolute running time: 7.28 sec, cpu time: 11.02 sec, memory peak: 55 Mb, absolute service time: 7,41 sec