packageを使わないでやってみた(2次元ニュートンラフソン)

#2dim-newton raphson
library(dplyr)
f<-deriv(~x^2/16+y^2/9-1,c("x","y"),function(x,y){})
g<-deriv(~x^2-y,c("x","y"),function(x,y){})
A=t(t(rep(9,2)))
times=100
val_f=val_g=df=dg=0
for(j in 1:times){
  a=A[1];b=A[2]
  df=attr(f(a,b),"gradient")
  dg=attr(g(a,b),"gradient")
  H=matrix(c(df,dg),ncol=2,nrow=2)
  A=A-((det(H)^(-1))*matrix(c(dg[1,2],-dg[1,1],-df[1,2],df[1,1]),ncol=2,nrow=2))%*%c(f(a,b)[1],g(a,b)[1])
}

A

#多次元勾配上昇法

library(dplyr)

d=3

f=function(x){

return(exp(-(x[1]-2)^2)+exp(-(x[2]-1)^2)+exp(-(x[3]-3)^2))

}


ite=100000

value=rep(0.1,d)

eta=0.01;h=0.01

for(j in 1:ite){

for(i in 1:d){

vec=rep(0,d);vec[i]=vec[i]+h  

value[i]=value[i]+eta*(f(value+vec)-f(value))/h  

}  

#print(value)  

}

print(value) 

#Levenberg-Marquardt method

library(dplyr)

data=data.frame(num=1:20,y=c(0.20, 0.10, 0.49, 0.26, 0.92, 0.95, 0.18, 0.19, 0.44, 0.79, 0.61, 0.41, 0.49, 0.34, 0.62, 0.99, 0.38, 0.22,0.71,0.40),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56),x3=c(24.5,24,23,22.5,22,24,25,26,25,24,23,21,26,25.5,24,23,24,24,24,22))

X=as.matrix(data[,colnames(data) %in% c("x1","x2")])/1000

Y=as.numeric(data$y)

n=nrow(data)

l1=function(x,y,a){

return((y-sum(exp(x*a))))  

}

l2=function(x,y,a){

return(sqrt(abs(y-sum(x*a))))  

} 

a=rep(1,nrow(X))





ite=10000

h=0.001;mu=10000

for(j in 1:ite){



a_pre=a  

df_vec=c();dl_vec=array(0,dim=c(nrow(X),ncol(X)))

ls=c()

for(i in 1:length(Y)){
df=l1(X[i,],Y[i],a[i])*(l1(X[i,]+h,Y[i],a[i])-l1(X[i,],Y[i],a[i]))/h+l2(X[i,],Y[i],a[i])*(l2(X[i,]+h,Y[i],a[i])-l2(X[i,],Y[i],a[i]))/h  

dl1=(l1(X[i,]+h,Y[i],a[i])-l1(X[i,],Y[i],a[i]))/h
dl2=(l2(X[i,]+h,Y[i],a[i])-l2(X[i,],Y[i],a[i]))/h

ls=c(ls,l1(X[i,],Y[i],a[i])+l2(X[i,],Y[i],a[i]))

dl_vec[i,]=c(dl1,dl2)

df_vec=c(df_vec,df)

}

J=dl_vec

ds=solve(J%*%t(J)+diag(mu,n))%*%(-df_vec)

a=a+ds

#print(sum(abs(a-a_pre)))

print(sum(ls^2))

}