-1

私は次のコードを持っています*。これにより、変数「epsf1」を以前の値「epsf1 = eps [ii、1] -c1*epsf1」に変更します。その値を出力すると0になります。しかし、完全に異なるものになるはずです。なぜそれが起こっているのかわかりませんし、同じ問題でより短いコードを作成することもできません。「eps[ii]」の値は「epsf1=eps [ii、1] -c1*epsf1;」では使用されていないようです。声明。

*：

delka_experimentu = 1000;
lambda  = 1;
Q = diag(3);

e = rnorm(delka_experimentu);
u = rnorm(delka_experimentu);#rep(10,delka_experimentu)
y = rep(0,delka_experimentu);
y[1] = e[1];

for(ii in c(2:delka_experimentu)){
 y[ii] = 0.5*y[ii-1]+0.3*u[ii-1]+e[ii]-0.2*e[ii];
}

Res =  matrix(0,delka_experimentu,6);


P = 1000*diag(3);
R1 = 0.1*diag(3);
K = matrix(0,3,1);
eps = matrix(0,delka_experimentu,1)
epsf1 = y[1];
print(c(epsf1))
epsf2 = 0;
ef3 = 0;
yf1 = y[1];
yf2 = 0;
yf3 = 0;
uf1 = u[1];
uf2 = 0;
uf3 = 0;
eps[1,1] = y[1];

epst1 = 0;
pomocna1 = 0;
pomocna2 = 0;
theta = matrix(0,3,1);

for(ii in 2:30){ #delka_experimentu
a1 = theta[1];
b1 = theta[2];
c1 = theta[3];

pomocna1 = epsf1;
epsf1 = eps[ii,1]-c1*epsf1;
epsf2 = pomocna1;
print(c(epsf1))

pomocna1 = yf1;
yf1 = y[ii]-c1*yf1;
yf2 = pomocna1;

pomocna1 = uf1;
uf1 = u[ii]-c1*uf1;
uf2 = pomocna1;

eps[ii,1] = y[ii]+a1*y[ii-1]-b1*u[ii-1]-c1*eps[ii-1]

#if(eps[ii]*eps[ii]>100){paste("problem",ii)}

psi = matrix(c(yf2,uf2,epsf2),3,1);

b = t(psi)%*%P%*%psi;
P = (P - P%*%psi%*%solve(1+b[1,1])%*%t(psi)%*%P)/lambda #lambda*solve(Q)
K = P%*%psi%*%solve(1+b[1,1]); #lambda*solve(Q)
theta = theta + K%*%eps[ii];
Res[ii,] = theta;
}

plot(c(0,delka_experimentu),c(min(Res,-0.5),max(Res,1)),col="white")
lines(Res[,1],col="red")
#lines(Res[,2],col="red",type=o)
lines(Res[,3],col="blue")
#lines(Res[,4],col="blue")
lines(Res[,5])
#lines(Res[,6])

lines(c(0,delka_experimentu),c(-0.5,-0.5),col="red")
lines(c(0,delka_experimentu),c(0.3,0.3),col="blue")
lines(c(0,delka_experimentu),c(-0.2,-0.2))
4

1 に答える 1

1

問題は、0 を掛けていることです。2 行目eps(つまりeps[ii,1]when ii=2) は 0 であり、theta0 の行列であるため、明らかepsf1に 0 になります。問題はeps[ii,1]、ループで定義する場所が後であることです。あなたが割り当てepsf1ます。そのため、ループを変更してその行を移動すると、変化する変数が出力されます。編集したコードは次のとおりです。

rm(list=ls())
delka_experimentu = 1000;
lambda  = 1;
Q = diag(3);

e = rnorm(delka_experimentu);
u = rnorm(delka_experimentu);#rep(10,delka_experimentu)
y = rep(0,delka_experimentu);
y[1] = e[1];

for(ii in c(2:delka_experimentu)){
  y[ii] = 0.5*y[ii-1]+0.3*u[ii-1]+e[ii]-0.2*e[ii];
}

Res =  matrix(0,delka_experimentu,6);


P = 1000*diag(3);
R1 = 0.1*diag(3);
K = matrix(0,3,1);
eps = matrix(0,delka_experimentu,1)
epsf1 = y[1];
print(c(epsf1))
epsf2 = 0;
ef3 = 0;
yf1 = y[1];
yf2 = 0;
yf3 = 0;
uf1 = u[1];
uf2 = 0;
uf3 = 0;
eps[1,1] = y[1];

epst1 = 0;
pomocna1 = 0;
pomocna2 = 0;
theta = matrix(0,3,1);

for(ii in 2:30){ #delka_experimentu
  a1 = theta[1];
  b1 = theta[2];
  c1 = theta[3];
  eps[ii,1] = y[ii]+a1*y[ii-1]-b1*u[ii-1]-c1*eps[ii-1] #This is the line I moved
  pomocna1 = epsf1;
  epsf1 = eps[ii,1]-c1*epsf1;
  epsf2 = pomocna1;
  print(c(epsf1))
  pomocna1 = yf1;
  yf1 = y[ii]-c1*yf1;
  yf2 = pomocna1;

  pomocna1 = uf1;
  uf1 = u[ii]-c1*uf1;
  uf2 = pomocna1;



  #if(eps[ii]*eps[ii]>100){paste("problem",ii)}

  psi = matrix(c(yf2,uf2,epsf2),3,1);

  b = t(psi)%*%P%*%psi;
  P = (P - P%*%psi%*%solve(1+b[1,1])%*%t(psi)%*%P)/lambda #lambda*solve(Q)
  K = P%*%psi%*%solve(1+b[1,1]); #lambda*solve(Q)
  theta = theta + K%*%eps[ii];
  Res[ii,] = theta;
}  
}

plot(c(0,delka_experimentu),c(min(Res,-0.5),max(Res,1)),col="white")
lines(Res[,1],col="red")
#lines(Res[,2],col="red",type=o)
lines(Res[,3],col="blue")
#lines(Res[,4],col="blue")
lines(Res[,5])
#lines(Res[,6])

lines(c(0,delka_experimentu),c(-0.5,-0.5),col="red")
lines(c(0,delka_experimentu),c(0.3,0.3),col="blue")
lines(c(0,delka_experimentu),c(-0.2,-0.2))

次の出力が得られます。

[1] 0.003594182
[1] -0.9812494
[1] -5.430992
[1] -27.87607
[1] -143.5392
[1] -732.6409
[1] -3729.887
[1] -18962.87
[1] -96278.02
[1] -488131.6
[1] -2471325
[1] -12494038
[1] -63075221
[1] -317971881
[1] -1600659248
[1] -8043157617
[1] -40412374164
[1] -203058922762
[1] -1.020209e+12
[1] -5.125226e+12
[1] -2.308228e+13
[1] -1.039444e+14
[1] -4.680327e+14
[1] -2.107201e+15
[1] -9.47394e+15
[1] -4.259006e+16
[1] -1.914032e+17
[1] -8.600152e+17
[1] -3.905738e+18
于 2013-02-28T12:25:53.737 に答える