拡散モデルーーデータ生成技術の数理　岡野原先生　岩波書店

```r
#p.8 積分値により行った場合
data("EuStockMarkets")

x=as.matrix(round(EuStockMarkets))

min_x=apply(x,2,min)
max_x=apply(x,2,max)


#f_z(x)=z1*x1+z2*x2+....

L=function(s){
val1=-mean(x%*%s)  
val2=0
for(i in 1:ncol(x)){
exp_val=function(r){
return(exp(-r*s[i]))  
}
val2=val2+log(integrate(exp_val,min_x[i],max_x[i])$value)  
}
val2=-val2
return(val1+val2)
}


ite=2*10^1
eta=10^(-4)
sita=rep(0.0003,ncol(x))

for(l in 1:ite){
val1=c();val2=c();
for(i in 1:length(sita)){
exp_val1=function(r){
return(exp(-r*sita[i]))  
}
exp_val2=function(r){
return(r*exp(-r*sita[i]))  
}
val1=c(val1,integrate(exp_val1,min_x[i],max_x[i])$value)
val2=c(val2,integrate(exp_val2,min_x[i],max_x[i])$value)  
}
sita=sita*eta*(-apply(x,2,mean)+val1/val2)

print(L(sita))  
}


```

```r

#p.8 
#data("EuStockMarkets")
data("LifeCycleSavings")

x=as.matrix((LifeCycleSavings))[,1:4]
y=as.matrix((LifeCycleSavings))[,5]

min_x=apply(x,2,min)
max_x=apply(x,2,max)


#f_z(x)=z1*x1+z2*x2+....

L=function(s){
val1=-mean(x%*%s)  
val2=-log(sum(exp(-x%*%s)))
return(val1+val2)
}


ite=10^4
eta=10^(-7)
sita=rep(0.3,ncol(x))

for(l in 1:ite){
sita=sita+eta*(-apply(x,2,mean)+t(x)%*%exp(-x%*%sita)/sum(exp(-x%*%sita)))
print(L(sita))  
}

#各標本の発生確率
exp(-x%*%sita)/sum(exp(-x%*%sita))


#ランジュバンMC(p(x)=exp(-x1*sita1+x2*sita2+...))
library(mvtnorm)

X=rmvnorm(1,rep(0,length(sita)),diag(1,length(sita)))
a=0.01;K=100
X_mat=array(0,dim=c(K+1,length(sita)))
X_mat[1,]=X

for(l in 1:K){
X=X-a*c(sita)+sqrt(2*a)*rmvnorm(1,rep(0,length(sita)),diag(1,length(sita)))
X_mat[(l+1),]=X
}


#p.16 スコアマッチング
#ここからは先程の確率のみしか知りえないとする
#この場合はスコアは不明である

#確率分布pとその標本データdf
p=exp(-x%*%sita)/sum(exp(-x%*%sita))
data("LifeCycleSavings")
x=as.matrix((LifeCycleSavings))[,1:4]
y=as.matrix((LifeCycleSavings))[,5]
prob_op=0

#exp(-(Y-X%*%sita)^2/(2*sig^2))をXで微分して、s(sita,sig,X)=sita*(Y-X%*%sita)/(sig^2)とする


S=function(s){
s1=s[1]
s2=s[-1]
return((y-x%*%s2)%*%s2/(2*s1^2))
}

trdS=function(s){
s1=s[1]
s2=s[-1]
return(rep(-sum(s2^2)/(2*s1^2),length(y)))  
}

JISM=function(s){#暗黙的スコアマッチングの値
if(prob_op==1){
return(sum(sqrt(apply(S(s)^2,1,sum))*p/2)+sum(trdS(s)*p))  
}else{  
return(mean(sqrt(apply(S(s)^2,1,sum))/2)+mean(trdS(s)))
}}

sita2=rep(0.01,5)
eta=0.0001
ite=1000
h=0.001

for(l in 1:ite){
print(JISM(sita2)) 
for(i in 1:length(sita2)){
sita2_sub=sita2;sita2_sub[i]=sita2_sub[i]+h
sita2[i]=sita2[i]-eta*(JISM(sita2_sub)-JISM(sita2))/h
}}


```

```r

#p.22　デノイジングスコアマッチング
#ほんとは過学習が起きそうな機械学習モデルなんかに用いると思う

#data("EuStockMarkets")
data("LifeCycleSavings")

x=as.matrix((LifeCycleSavings))[,1:4]
y=as.matrix((LifeCycleSavings))[,5]

min_x=apply(x,2,min)
max_x=apply(x,2,max)


#f_z(x)=z1*x1+z2*x2+....

L=function(s){
val1=-mean(x%*%s)  
val2=-log(sum(exp(-x%*%s)))
return(val1+val2)
}


ite=10^4
eta=10^(-7)
sita=rep(0.3,ncol(x))

for(l in 1:ite){
sita=sita+eta*(-apply(x,2,mean)+t(x)%*%exp(-x%*%sita)/sum(exp(-x%*%sita)))
print(L(sita))  
}

#各標本の発生確率
p=exp(-x%*%sita)/sum(exp(-x%*%sita))
sig=3

px=function(x1,x2){
return(exp(-sum((x1-x2)^2)/(2*sig^2))/(sqrt(2*pi*sig^2)^(ncol(x)/2)))
}

pxx=function(xval,l){
px_val=0
dpx_val=0
for(i in 1:nrow(x)){
px_val=px_val+px(xval,x[i,])*p[i] 
dpx_val=dpx_val+(x[i,]-xval)*px(xval,x[i,])*p[i]
}
if(l==1){
return(px_val)
}
if(l==2){
return(dpx_val/px_val)  
}
}

#学習における標本生成
library(mvtnorm)
N=100
X=array(0,dim=c(N,ncol(x)))
XX=array(0,dim=c(N,ncol(x)))
Y=c()
px_vec=c()
dpx_vec=array(0,dim=c(N,ncol(x)))
for(j in 1:N){
no=sample(1:nrow(x),1) 
Y=c(Y,y[no])
XX[j,]=x[no,]
X[j,]=x[no,]+rmvnorm(1,rep(0,ncol(X)),diag(sig^2,ncol(X)))
px_vec=c(px_vec,pxx(X[j,],1))
dpx_vec[j,]=pxx(X[j,],2)
}

S=function(s){
return((Y-X%*%s)%*%s/(2*sig^2))
}

JESM=function(s){
return(sum(apply((dpx_vec-S(s))^2,1,sum)*px_vec)/2)
}

#最急降下法
sita2=rep(0.01,ncol(X))
eta=0.001
ite=1000
h=0.001

for(l in 1:ite){
print(JESM(sita2)) 
for(i in 1:length(sita2)){
sita2_sub=sita2;sita2_sub[i]=sita2_sub[i]+h
sita2[i]=sita2[i]-eta*(JESM(sita2_sub)-JESM(sita2))/h
}}

```