LSTMをPythonで

最近、RNNについて勉強しています。今ちょうど、LSTMについて勉強しているところなのですが、tensorflowやkerasの使い方はなんとなく分かるけれど応用できるほどの技量はないので、イメージを膨らますためにもとりあえずPythonのnumpyのみでコードを書いてみようと思いました。一応、入出力ゲート、忘却ゲート、CEC、覗き穴結合などを搭載した自作LSTMを作ってみたのですが、うまく動きません。うまく動かないというのは、出力がnanになってしまいます。多分、重みを更新する際にすでに重み自体がnanになっているのではないかと考えられます。そこで、以下のコードで怪しい点などありますか？（個人の見解としては誤差逆伝播あたりが原因な気がしてます。）
また、そもそもtensorflowなどを使わないのは無謀ですか？そこらへんまだよくわからないので教えていただきたいです。

---コードについて---
alpha、gammaは学習率、equation()はy=xの線形変換をしてるだけなので意味はないです。（見やすかったので...）tanhとSigmoidの微分は頭にdをつけてます。また、各行列の宣言は行っていますが、長くなるので省きました。重みの初期値は平均0、標準偏差0.1の正規分布です。バイアスの初期値は0行列です。

コード

        #input gate
        i_[:, t] = np.dot(w[0, t, :, :], x[:, t]) + np.dot(u[0, t, :, :], h[:, t-1]) + np.dot(v[0, t, :, :], c[:, t-1]) + b[0, :, t]
        i[:, t] = Sigmoid(i_[:, t])
        #forget gate
        f_[:, t] = np.dot(w[2, t, :, :] , x[:, t]) + np.dot(u[2, t, :, :], h[:, t-1]) + np.dot(v[2, t, :, :], c[:, t-1]) + b[2, :, t]
        f[:, t] = Sigmoid(f_[:, t])
        #activated
        a_[:, t] = np.dot(w[3, t, :, :] , x[:, t]) + np.dot(u[3, t, :, :], h[:, t-1]) + b[3, :, t]
        a[:, t] = tanh(a_[:, t])
        #CEC
        c[:, t] = i[:, t]*a[:, t] + f[:, t]*c[:, t-1]
        #output gate
        o_[:, t] = np.dot(w[1, t, :, :] , x[:, t]) + np.dot(u[1, t, :, :], h[:, t-1]) + np.dot(v[1, t, :, :], c[:, t]) + b[1, :, t]
        o[:, t] = Sigmoid(o_[:, t])
        #hidden 
        h[:, t] = o[:, t] * tanh(c[:, t])
        #output layer
        q[:, t] = np.dot(vv[t, :, :], h[:, t]) + d1[:, t]
        #output
        y[:, t] = equation(q[:, t])

        e_h = 1
        e_c[:, t] = e_h * o[:, t] * dtanh(c[:, t])
        e_c[:, t-1] = e_c[:, t] * f[:, t]
        e_o[:, t] = y[:, t] - T[:, t]

        e_i_[:, t] = e_c[:, t] * a[:, t] * dSigmoid(i_[:, t])
        e_o_[:, t] = e_h * tanh(c[:, t]) * dSigmoid(o_[:, t])
        e_f_[:, t] = e_c[:, t] * c[:, t-1] * dSigmoid(f_[:, t])
        e_a_[:, t] = e_c[:, t] * i[:, t] * dtanh(a_[:, t])

        dw[0, t, :, :] = np.dot(e_i_[:, t], x[:, t])
        du[0, t, :, :] = np.dot(e_i_[:, t], h[:, t-1])
        dv[0, t, :, :] = np.dot(e_i_[:, t], c[:, t-1])
        db[0, :, t] = e_i_[:, t]

        dw[1, t, :, :] = np.dot(e_o_[:, t], x[:, t])
        du[1, t, :, :] = np.dot(e_o_[:, t], h[:, t-1])
        dv[1, t, :, :] = np.dot(e_o_[:, t], c[:, t])
        db[1, :, t] = e_o_[:, t]

        dw[2, t, :, :] = np.dot(e_f_[:, t], x[:, t])
        du[2, t, :, :] = np.dot(e_f_[:, t], h[:, t-1])
        dv[2, t, :, :] = np.dot(e_f_[:, t], c[:, t-1])
        db[2, :, t] = e_f_[:, t]

        dw[3, t, :, :] = np.dot(e_a_[:, t], x[:, t])
        du[3, t, :, :] = np.dot(e_a_[:, t], h[:, t-1])
        db[3, :, t] = e_a_[:, t]

        dvv[t, :, :] = np.dot(e_o[:, t], h[:, t])
        dd1[:, t] = e_o[:, t]

        w[:, t+1, :, :] = w[:, t, :, :] - alpha * dw[:, t, :, :]
        u[:, t+1, :, :] = u[:, t, :, :] - alpha * du[:, t, :, :]
        v[:, t+1, :, :] = v[:, t, :, :] - alpha * dv[:, t, :, :]
        b[:, :, t+1] = b[:, :, t] - alpha * db[:, :, t]

        vv[t+1, :, :] = vv[t, :, :] - gamma * dvv[t, :, :]
        d1[:, t+1] = d1[:, t] - gamma * dd1[:, t]