`
#p.26 修正中です
N=10
x=c(1:N)
u=sample(1:100,N,replace=T)
a_b=array(0,dim=c(2,N-1))
lam=1
integral=c()
for(i in 1:(N-1)){
u_vec=u[i:(i+1)];x_vec=x[i:(i+1)]
mat=cbind(rep(1,2),x_vec)
alpha=solve(mat)%*%u_vec
a_b[,i]=alpha
u_f=function(x){
N=c((x_vec[2]-x)/(x_vec[2]-x_vec[1]),(x-x_vec[1])/(x_vec[2]-x_vec[1]))
return(lam*sum(N*u_vec))
}
u_f2=function(x){
l_k=c(-1/(x_vec[2]-x_vec[1]),1/(x_vec[2]-x_vec[1]))
N=c((x_vec[2]-x)/(x_vec[2]-x_vec[1]),(x-x_vec[1])/(x_vec[2]-x_vec[1]))
return(as.numeric(l_k%*%t(t(u_vec)))*as.numeric(N%*%u_vec))
}
#monte carlo integral
Y1=0;Y2=0
n=100000
x1=sample(seq(x_vec[1],x_vec[2],0.001),n,replace=T);x2=sample(seq(x_vec[1],x_vec[2],0.001),n,replace=T)
for(j in 1:(n)){
Y1=Y1+u_f(x1[j]);
Y2=Y2+u_f2(x2[j])
}
#monte carlo積分値
print(Y1/n);print(Y2/n)
integral=c(integral,Y1/n+Y2/n)
}
u_dif=u[2:length(u)]-u[1:(length(u)-1)]
cor(integral,u_dif)
#p.48
x=c(0,1,2);y=c(1,0,2)
pthi=sample(1:10,3)/10
mat=cbind(rep(1,3),x,y)
a=solve(mat)%*%pthi
delta=det(mat)/2
A1=(x[2]*y[3]-x[3]*y[2])/(2*delta)
A2=(x[3]*y[1]-x[1]*y[3])/(2*delta)
A3=(x[1]*y[2]-x[2]*y[1])/(2*delta)
B1=(y[2]-y[3])/(2*delta)
B2=(y[3]-y[1])/(2*delta)
B3=(y[1]-y[2])/(2*delta)
C1=(x[3]-x[2])/(2*delta)
C2=(x[1]-x[3])/(2*delta)
C3=(x[2]-x[1])/(2*delta)
#N1==1
A1+B1*x[1]+C1*y[1]
#N2==1
A2+B2*x[2]+C2*y[2]
#N3==1
A3+B3*x[3]+C3*y[3]
B=c(B1,B2,B3);C=c(C1,C2,C3)
#(3.33)
G=(t(t(B))%*%t(B)+t(t(C))%*%t(C))*delta
#重みの最適化(pthi_e,pthi_k)でよいか？
Q_k=c(10,20)
#境界をc(0,1)->c(2,2)
L_k=sqrt((x[3]-x[2])^2+(y[3]-y[2])^2)
H=L_k*matrix(c(2,1,1,2),ncol=2)/6
pthi_e=rep(1,3);pthi_k=rep(1,2)
cost=function(pthi1,pthi2){
return((as.numeric(t(c(pthi1))%*%G%*%pthi)-as.numeric(t(c(pthi2))%*%H%*%Q_k))^2)
}
ite=100
r=0.001
for(l in 1:ite){
vec=pthi_e
for(j in 1:length(pthi_e)){
vec_sub=vec;vec_sub[j]=vec_sub[j]+0.01
pthi_e[j]=pthi_e[j]-r*(cost(vec_sub,pthi_k)-cost(vec,pthi_k))/0.01
}
vec=pthi_k
for(j in 1:length(pthi_k)){
vec_sub=vec;vec_sub[j]=vec_sub[j]+h
pthi_k[j]=pthi_k[j]-r*(cost(pthi_e,vec_sub)-cost(pthi_e,vec))/0.01
}
print(cost(pthi_e,pthi_k))
}