確率制御理論入門　コロナ社　K.J.アストロノーム著

```r

#確率制御理論入門　p.141
#システム変数の分散の計算式


a=c(1,0.7,0.5,-0.3)

b=c(1,0.3,0.2,0.1)

a0=a[1]

amat=array(0,dim=c(length(a),length(a)))

bmat=array(0,dim=c(length(b),length(b)))

alphavec=c()

betavec=c()

bvec=c()

for(j in 1:length(a)){
  
bvec=c(bvec,b[length(b)])  

amat[j,]=c(a,rep(0,(j-1)))

bmat[j,]=c(b,rep(0,(j-1)))

inva=a[length(a):1]

beta=b[length(b)]/a[1]

if(j<nrow(amat)){

alpha=inva[1]/a[1]

}else{
  
alpha=a  
  
}

betavec=c(betavec,beta)

alphavec=c(alphavec,alpha)

b=(b-inva*beta)[1:(length(b)-1)]

a=(a-inva*alpha)[1:(length(a)-1)]

}

I=sum(bvec^2/amat[,1])/a0

```

伝達関数の平滑化モデル(＃安定性の考慮が必要)：
$$ H(z)=\prod_{i = 1}^{K}\frac{\sum_{l=0}^{P}a_{l}z^{-l}}{\sum_{l=0}^{P}b_{l}z^{-l}} \space \space \space ただし、a_{0}=1、b_{0}=1$$

をデータ$H_d(e^{jw_{i}})$について二乗誤差と平滑化を加えた：
$$ Q(\theta)=\sum_{i=1}^{M}(|H(e^{jw_{i}})|-H_d(e^{jw_{i}}))^{2}+\zeta \sigma^{2} $$

を最小にするようにアルゴリズムを作成した。ただし、$\zeta$は平滑化パラメータであり、

$$\sigma^{2}=\int^{\pi}_{-\pi} H(z)H(z^{-1}) dz $$

である。$z=e^{jw}、\eta(w)=H(e^{jw})H(e^{-jw})=|H(e^{jw})|^{2}$はスペクトル密度関数と呼ばれるものであり、$\sigma^2$はこの変動の総計を表すものである。

$\theta$はパラメータで$\theta=(a_{1},...,a_{P},b_{1},...,b_{P})$である。

$\zeta=0.5での出力例:\sigma^{2}=0.05370661、学習率0.05$

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/230981/7443bf3b-4227-ed4b-cdb2-246eb5998576.png)

$\zeta=0.001での出力例:\sigma^{2}=18.28、学習率0.03$

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/230981/371489c9-0b54-45d6-ada8-793a1a026dcd.png)

```r


#p.132　応用

w=c(-4.1	,
-4.099375	,
-4.09875	,
-4.098125	,
-4.0975	,
-4.096875	,
-4.09625	,
-4.095625	,
-4.095	,
-4.094375	,
-4.09375	,
-4.093125	,
-4.0925	,
-4.091875	,
-4.09125	,
-4.090625	,
-4.09	,
-4.089375	,
-4.08875	,
-4.088125	,
-4.0875	,
-4.086875	,
-4.08625	,
-4.085625	,
-4.085	,
-4.084375	,
-4.08375	,
-4.083125	,
-4.0825	,
-4.081875	,
-4.08125	,
-4.080625	,
-4.08	,
-4.079375	,
-4.07875	,
-4.078125	,
-4.0775	,
-4.076875	,
-4.07625	,
-4.075625	,
-4.075	,
-4.074375	,
-4.07375	,
-4.073125	,
-4.0725	,
-4.071875	,
-4.07125	,
-4.070625	,
-4.07	,
-4.069375	,
-4.06875	,
-4.068125	,
-4.0675	,
-4.066875	,
-4.06625	,
-4.065625	,
-4.065	,
-4.064375	,
-4.06375	,
-4.063125	,
-4.0625	,
-4.061875	,
-4.06125	,
-4.060625	,
-4.06	,
-4.059375	,
-4.05875	,
-4.058125	,
-4.0575	,
-4.056875	,
-4.05625	,
-4.055625	,
-4.055	,
-4.054375	,
-4.05375	,
-4.053125	,
-4.0525	,
-4.051875	,
-4.05125	,
-4.050625	,
-4.05	,
-4.049375	,
-4.04875	,
-4.048125	,
-4.0475	,
-4.046875	,
-4.04625	,
-4.045625	,
-4.045	,
-4.044375	,
-4.04375	,
-4.043125	,
-4.0425	,
-4.041875	,
-4.04125	,
-4.040625	,
-4.04	,
-4.039375	,
-4.03875	,
-4.038125	
)
w=max(w)-w
w=2*pi*w*10

Hd=c(0.015997742	,
0.015997742	,
0.016329822	,
0.015333584	,
0.06580965	,
-0.012893164	,
-0.023519705	,
0.057839745	,
0.002050408	,
0.012344869	,
0.058835983	,
0.002382487	,
0.016329822	,
0.00105417	,
0.010684472	,
0.023967648	,
0.042564093	,
-0.003262862	,
-0.051414373	,
0.038247061	,
0.020314774	,
0.036918744	,
0.088723128	,
0.016329822	,
0.121266908	,
0.056843507	,
-0.139415411	,
0.09337224	,
-0.077648645	,
-0.115505695	,
0.068466286	,
0.046216967	,
0.041235776	,
-0.036470801	,
0.084406096	,
0.088723128	,
-0.033814166	,
-0.003927021	,
0.030609235	,
-0.077980724	,
0.013341107	,
0.080089064	,
-0.012561085	,
-0.002930783	,
0.021643092	,
0.012344869	,
0.048209443	,
0.044888649	,
-0.040787833	,
0.043560331	,
0.036918744	,
-0.017542276	,
0.039243299	,
-0.021859308	,
-0.016546038	,
0.020978933	,
0.014337345	,
0.023303489	,
-0.010236529	,
0.047545284	,
0.091047684	,
-0.051082294	,
0.012676948	,
0.10399878	,
-0.021859308	,
-0.012229006	,
0.014669425	,
-0.077648645	,
0.052858554	,
0.21225666	,
0.048541522	,
-0.184578208	,
0.068134206	,
0.052858554	,
-0.10753579	,
-0.024183864	,
0.008691996	,
0.010020313	,
0.027288441	,
0.018654377	,
0.041899935	,
0.014669425	,
0.002714567	,
0.02263933	,
0.038247061	,
0.050866078	,
0.01168071	,
0.007363678	,
-0.013889403	,
0.003046646	,
0.061492618	,
-0.039791595	,
0.010020313	,
0.107651653	,
-0.013225244	,
0.015665663	,
0.030277156	,
-0.057059723	,
0.050866078	,
-0.006583656
)

w=w[1:60]
Hd=Hd[1:60]

M=length(w)
L=20

#Gauss法

f=function(x){
ak=x[1:(L)];x=x[-c(1:(L))]
bk=x
val1=1
val2=1
for(j in 1:L){
val1=val1+ak[j]*exp(w*1i)^(-j)
val2=val2+bk[j]*exp(w*1i)^(-j)  
}  
val=val1/val2
val=sqrt(Re(val)^2+Im(val)^2)
return(val)
}

#システム変数の分散
sig=function(x){
a=x[1:(L)];x=x[-c(1:(L))]
b=x
a0=a[1]
amat=array(0,dim=c(length(a),length(a)))
bmat=array(0,dim=c(length(b),length(b)))
alphavec=c()
betavec=c()
bvec=c()
for(j in 1:length(a)){
bvec=c(bvec,b[length(b)])  
amat[j,]=c(a,rep(0,(j-1)))
bmat[j,]=c(b,rep(0,(j-1)))
inva=a[length(a):1]
beta=b[length(b)]/a[1]
if(j<nrow(amat)){
alpha=inva[1]/a[1]
}else{
alpha=a  
}
betavec=c(betavec,beta)
alphavec=c(alphavec,alpha)
b=(b-inva*beta)[1:(length(b)-1)]
a=(a-inva*alpha)[1:(length(a)-1)]
}
I=sum(bvec^2/amat[,1])/a0
return(I)
}

Y=abs(Hd)
h=0.001
lam=1
zeta=0.001
ite=5*10^3

sita=sample(seq(1,100,1),2*L,replace=T)/100
eta=10^(-2)

for(l in 1:ite){
sita_pre=sita
Z=array(0,dim=c(length(Y),length(sita)))
vec=sita
for(j in 1:ncol(Z)){
vec_sub=vec;vec_sub[j]=vec_sub[j]+h
Z_sub=(f(vec_sub)-f(vec))/h
Z[,j]=Z_sub
}
beta=solve(t(Z)%*%Z+diag(lam,ncol(Z)))%*%t(Z)%*%(Y-f(sita))

dparam=c()
for(j in 1:length(sita)){
paramsub=sita
paramsub[j]=paramsub[j]+h
dparam=c(dparam,(sig(paramsub)-sig(sita))/h)
}
sita=c(sita+eta*beta-eta*zeta*dparam)
print(sum((Y-f(sita))^2)+zeta*sig(sita))
}


plot(w,f(sita),xlim=c(min(w),max(w)),ylim=c(min(c(f(sita),Y)),max(c(f(sita),Y))),xlab="w",ylab="H",col=2)

par(new=T)

plot(w,Y,xlim=c(min(w),max(w)),ylim=c(min(c(f(sita),Y)),max(c(f(sita),Y))),xlab="w",ylab="H",col=3)




```