配列のループ処理の出力がおかしい

Question

現在、c言語で行ったプログラムをpythonに変換してプログラムを書こうとしています。
目的としている処理はfor文の中で y, s, faiの値を更新していくプログラムです。

1. 初期値としてｘはすべて既知のデータでｗの初期値のみが与えられています。
2. ｘとｗの初期値をもとにｙを求める。
3. ｙをアークタンジェントとしたものをfaiとします。
4. yとfaiをプログラムに書かれている式に挿入し次ｗを求める。
5. 2~4を繰り返すことでｗを更新してゆく

図で示しているのは二行に列ｗの要素の一つを取り出したものです。
しかし、結果は図のようになり結果が更新されないうえに途中から0に収束しています。
手順4においてeta以降の数式を定数倍などに変更すると値が上昇することから、配列のいい関係などは間違っていないと思うのですが、解決方法がわかりません。
まだ、pythonに関しては初心者であるため読みにくいプログラムですが、教えていただけると幸いです。

#include<iostream>
#include<fstream>
#include<vector>
#include<cmath>
#include<stdio.h>

using namespace std;

#define N 80000 //データ点数
#define eta 0.0001 //ステップサイズ
#define signal1 "origin_data/a04.txt"
#define signal2 "origin_data/a27.txt"

void Load(vector<double>&s,string signal) //ファイルの読み込み
{
    double temp;
    ifstream ist(signal);
    for(int i=0;i<N;++i)
    {
        ist>>temp;
        s.push_back(temp);
    }
}

void Normalization(vector<double> &s) //原信号の正規化
{
    double sum = 0, avg, stdev;
    for(int i=0;i<N;++i)
    {
        sum += s[i];
    }
    avg = sum/N;
    for(int i=0;i<N;++i)
    {
        sum += pow(s[i] - avg,2);
    }
    stdev = sqrt(sum / N); //標準偏差を求める
    for(int i=0;i<N;++i)
    {
        s[i] = (s[i] - avg) / stdev; //正規化を行う
    }
}

void Print(vector<double>x1, vector<double>x2) //ファイルの出力
{
    ofstream ost1("output/a.txt");
    ofstream ost2("output/b.txt");
    for (int i = 0; i < N; i++)
    {
        ost1 << x1[i] << endl;
        ost2 << x2[i] << endl;
    }
}


int main()
{
    vector<double> s1,s2;
    double a1[2] = {1.0, 0.5}; //インパルス応答
    double a2[2] = {0.3, 2.0};

    Load(s1, signal1);
    Normalization(s1);

    Load(s2, signal2);
    Normalization(s2);

//混合信号の作成
    vector<double> x1, x2;
    for(int i=0;i<N;++i)
    {
        x1.push_back(a1[0]*s1[i] + a1[1]*s2[i]);
        x2.push_back(a2[0]*s1[i] + a2[1]*s2[i]);
    }

//分離更新
    vector<double>y1,y2; //分離信号
    vector<double>w11,w12,w21,w22; //分離行列
    vector<double>fai1,fai2; //スコア関数
    w11.push_back(1);
    w12.push_back(0);
    w21.push_back(0);
    w22.push_back(1);
    for (int i = 0; i < N; i++)
    {
        y1.push_back(w11[i]*x1[i] + w12[i]*x2[i]);
        y2.push_back(w21[i]*x1[i] + w22[i]*x2[i]);
        fai1.push_back(atan(y1[i]));
        fai2.push_back(atan(y2[i]));
        w11.push_back(w11[i] + eta * ((1-fai1[i]*y1[i])*w11[i] - fai1[i]*y2[i]*w21[i]));
        w12.push_back(w12[i] + eta * ((1-fai1[i]*y1[i])*w12[i] - fai1[i]*y2[i]*w22[i]));
        w21.push_back(w21[i] + eta * ((1-fai2[i]*y2[i])*w21[i] - fai2[i]*y1[i]*w11[i]));
        w22.push_back(w22[i] + eta * ((1-fai2[i]*y2[i])*w22[i] - fai2[i]*y1[i]*w12[i]));
    }

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as sp
import math

N, eta = 80000, 0.0001
A = np.array([[1.0, 0.5], [0.3, 2.0]])
s1 = np.loadtxt('origin_data/a04.txt')
s2 = np.loadtxt('origin_data/a27.txt')

s1 = sp.zscore(s1)
s2 = sp.zscore(s2)

x = np.array([[0 for i in range(N)] for j in range(2)])
for i in range (N):
    for j in range(2):
        x[j][i] = A[j][0]*s1[i] + A[j][1]*s2[i]

y = np.array([[0 for i in range(N)] for j in range(2)])
w = np.array([[[0 for i in range(N)] for j in range(2)] for k in range(2)])
fai = np.array([[0 for i in range(N)] for j in range(2)])
w[0][0][0], w[0][1][0], w[1][0][0], w[1][1][0]  = 1, 0, 0, 1

for j in range(2):
    y[j][0] = w[j][0][0]*x[0][0] + w[j][1][0]*x[1][0]
    fai[j][0] = math.atan(y[j][0])

for i in range(N-1):
    for j in range(2):
        y[j][i+1] = w[j][0][i]*x[0][i] + w[j][1][i]*x[1][i]
        fai[j][i+1] = math.atan(y[j][i])
        w[0][j][i+1] = w[0][j][i] + eta * ((1-fai[0][i]*y[0][i])*w[0][j][i] - fai[0][i]*y[1][i]*w[1][j][i])
        w[1][j][i+1] = w[1][j][i] + eta * ((1-fai[1][i]*y[1][i])*w[1][j][i] - fai[1][i]*y[0][i]*w[0][j][i])

Answer 1

ちゃんと見ていないのでこれが原因かはちょっと自信ないのですが、
Pythonに変更する際、yとfaiとwを、y1, y2という別変数から配列にしているようです。
yとfaiはいいのですが、元々はwは1つ(実際には4つ)であるものを2つ(実際には8つ)に増やされています。
また、wxxの値の計算にはy1とy2両方が必要なはずですが、jのループ内でwを計算しているため、少なくともjが0の時には計算に使用する値がおかしいものと思われます。

wは2次元配列に戻し、jのループではなくiのループで計算したら、またちょっと変わるかもしれません。

蛇足ですが、入力値が必要なプログラムの質問の時は、現象が再現するテストデータと、
それを実行した結果があると確認がしやすいです。

Comments

[itsukidesu]

回答ありがとうございます。wの数が増えていることに気づきませんでした。配列に対するイメージが足りていなかったようなので指摘の通りループ数を減らして、配列も次元を落とすことで求めていた結果が得られました。
質問に関するアドバイスもありがとうございます。回答が得られないこともあるのはそのような記述が足りていないためだお思うので今後追加していきたいと思います。