Run Code
|
API
|
Code Wall
|
Misc
|
Feedback
|
Login
|
Theme
|
Privacy
|
Patreon
16-09-2020ExVitorb
# 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=c(2,5,8,11,14,17,27,31,41,44) n=length(x) y=c(94.8,89.7,81.3,74.9,68.7,64.0,49.3,44.0,39.1,31.6) plot(x,y,col="red") # Suponha que y=f(x) (aproximado) e que f(s)= s/(a+bs). Descobrir aproximações para 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). # # 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] p=0 for (i in 1:n){ p=p+( y[i]-x[i]/(a+b*x[i]) )^2} p/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 } # #---- 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 } # Teste de aproximação de solução do problema de minimização. u0=c(-0.16253301, 0.03055026); t=20; m=200000 print("Erro na condição inicial") ; g(u0) # Teste de g(u0). u=ZeroEuler(u0,t,m) ; print("Solução"); u # Aproximadamente u(t) print("Erro na aproximação"); g(u) # Teste de minímo aproximado # Definição da aproximação de f(u). fap<-function(s){ s/(u[1]+u[2]*s)} curve(fap,2,44) #points(x,y,col="red") # Teste de ajuste. # Se quiser pode aumentar t e n ou escolher outra condição inicial # para se convencer do resultado obtido. # Ajuste linearizado # # Mudança de variável z=x/y=~ x/(f(x)) = a+bx. Determinar a e b. # Minimizar |Au-z|^2, em que A matriz, cuja linha i é o vetor c(1,x_i). # Resolver o sistema linear A^*(Au-z)=A^*Au-A^*z=0 pelo método A=QR. A=matrix(0,n,2); v=0*x; for ( i in 1:n){A[i,]=c(1,x[i]);v[i]=x[i]/(y[i])} # Agora é conhecido que o gradiente de g(u) é dado por A^*(Au-v). # # Como o problema é linear, vamos resolver pelo Método de decomposição A=QR. # Isso nos leva a resolver a equação Ru=t(Q)v. # ---------------Processo de ortogonalização de Gram-Schmdit pe<-function(u,v){ # Dados os vetores u e v, somo o produto de cada coordenada de u com a respectiva cooordenada de v. n=length(u) p=u[1]*v[1] for ( i in 2:n){p=p+u[i]*v[i]} p } mod<-function(u){sqrt(pe(u,u))} # módulo de u proj<-function(u,v){pe(u,v)*v} # Projeção de u sobre v (ortogonal), com norma de v igual 1 GramS<-function(A){ # Ortogonliza as colunas de A n=length(A[1,]) Q=0*A Q[,1]=A[,1]/mod(A[,1]) for (i in 2:n){ p=A[,i] for (j in 1:(i-1)){ p=p-proj(A[,i],Q[,j]) } Q[,i]=p/mod(p) } Q } Q=GramS(A); print("Matriz Q"); Q t(Q)%*%Q # Teste de ortogonalização R=t(Q)%*%A; print("Matriz R"); R #------ Resolução de sistema triangular Ru=t(Q)v. TSup<-function(A,b){ p=0*b; n=length(b) p[n]=b[n]/A[n,n] for (i in (n-1):1){ q=b[i] for (j in n:(i+1)){ q=q-A[i,j]*p[j] } p[i]=q/A[i,i] } p } v=t(Q)%*%(x/y);v u=TSup(R,v) #--------------------------------------------------------------------------------------- print("Soulução") ; u # Aproximação para u que minimisa |Au-v|^2. fapmv<-function(s){s/(u[1]+u[2]*s)} # Aproximação da função linearizada. curve(fapmv,2,44) # Teste de ajuste points(x,y,col="red")
run
|
edit
|
history
|
help
0
Linear regression wall
T4_EX01_9436451
ECO Stat
AjusteLinearizado
Linear Regression and log-log Transformation
R added
Ex11-03-21EdsonAmanda
Test
Hiring Discrimination
Predictive power of a test when False positive rate = False negative rate