ベイズ回帰

まずは普通の回帰分析†

• 根本さんとかの振動解析モデルで、端部(\( x_{1j} \))とか中央支柱接合部(\( x_{2j} \))、、、と、 ヤング率を低下させる箇所を\( x_{ij} \)の\( i \)で表す。
• 箇所\( i \)ごとに、ヤング率を変化させた具体的な数値を\( x_{11}=100 \)MPa, \( x_{12}=200 \)MPa, \( x_{13}=300 \)MPa....みたいに、 \( j=1,2,3,... \)に対して与えていく。
• 全ての欠損箇所(\( i \))に\( j=1 \)のヤング率を与えて、振動解析を行って求まった固有振動数を\( y_{1} \)に代入する。
• \( j=1,2,3,... \)に対してこれを行えば\( y_{j} \)が求まる。
• \( y_{j} \)を目的変数、\( x_{ij} \)を説明変数として重回帰分析を行う。
• \( y=b_{1}x_{1}+b_{2}x_{2}+...+b_{0} \)みたいな重回帰式が求まる。
• この\( x_{1}, x_{2}, ... \)の説明変数のうち、どれが\( y \)に対する影響が強いのか弱いのかを、偏回帰係数の検定で調べるとかもできるか？
  • 重回帰分析のツールでF値とかが計算されるやつを利用できるかも。まずは、下の感度解析をやってみる。
• 感度解析のやり方
  • 根本さんとかの振動解析モデルで、端部(\( x_{1j} \))とか中央支柱接合部(\( x_{2j} \))、、、と、 ヤング率を低下させる箇所を\( x_{ij} \)の\( i \)で表す。
  • ある1箇所の\( i \)に対して、ヤング率を変化させた具体的な数値を\( x_{i1}=100 \)MPa, \( x_{i2}=200 \)MPa, \( x_{i3}=300 \)MPa....みたいに、 \( j=1,2,3,... \)に対して与えていく。
  • ある1箇所の\( i \)に対して、\( j=1 \)のヤング率を与えて、振動解析を行って求まった固有振動数を\( y_{1} \)に代入する。
  • \( j=1,2,3,... \)に対してこれを行えば\( y_{j} \)が求まる。
  • \( y_{j} \)を目的変数、ある特定の\( i \)の\( x_{ij} \)を説明変数として単回帰分析を行う。
  • \( y=bx+b_{0} \)みたいな回帰式が求まる。
  • 上記のように箇所\( i \)だけのヤング率を変えて単回帰を行った場合の決定係数\( R^{2} \)を縦軸に、\( i \)を横軸にプロットすると、どの箇所\( i \)が最も影響があるかが推定できる。

固有振動数の変化を説明変数、腐朽箇所のヤング率の変化を目的変数とした回帰†

• 青山さんや青野さんがやっているのは、特定の腐朽箇所のヤング率の変化を説明変数として、 特定の振動モードの固有振動数の変化を目的変数とする単回帰で、その回帰直線の傾きから感度を求めている。
  • 青野さんは、複数箇所が同時に腐朽した場合の腐朽箇所のヤング率の変化を説明変数とした場合についても同様の単回帰を行っている。
  • 今あるデータですぐにできそうなこととして、複数の腐朽箇所のヤング率の変化を複数の説明変数として、1つの振動モードの固有振動数の変化を目的変数として重回帰分析を行ってみるとうことはやってもいいかもしれない。
  • そうすると、どことどこの腐朽箇所の組み合わせが、ある特定の振動モードに対して感度が大きいかということはわかるかもしれない。
• でも、我々が最終的に知りたいのは、振動測定して得られた、複数のモードに対する固有振動数の変化（初期の健全状態は数値解析での予測値として）から、どこが腐朽している（確率が高い）かを予測することだ。
  • 今、青野さんが大量に計算したおかげで、複数の腐朽箇所を1箇所ずつ腐朽させていったときに、複数の振動モードのそれぞれが、 どれくらい変化するかという大量のデータがある。
  • では、回帰分析の説明変数と目的変数を逆にして、つまり、複数のモード（逆対称1次とか、水平対称1次とか、ねじり1次とか、それぞれの2次とか）の固有振動数の変化を説明変数にして、特定の腐朽箇所のヤング率の変化を目的変数にした重回帰分析をしたらどうだろうか。
    • \( y_{腐朽箇所1のE}=a_{11}x_{水平1次の\omega}+a_{12}x_{鉛直逆対称の\omega}+a_{13}x_{ねじれの\omega}+a_{14}x_{鉛直対称の\omega} \)
    • \( y_{腐朽箇所2のE}=a_{21}x_{水平1次の\omega}+a_{22}x_{鉛直逆対称の\omega}+a_{23}x_{ねじれの\omega}+a_{24}x_{鉛直対称の\omega} \)
    • ....
    • \( y_{腐朽箇所nのE}=a_{n1}x_{水平1次の\omega}+a_{n2}x_{鉛直逆対称の\omega}+a_{n3}x_{ねじれの\omega}+a_{n4}x_{鉛直対称の\omega} \)

みたいな重回帰式が、腐朽箇所の個数だけできて、測定値とか、テストケースの固有振動数の変化を上記の\( x \)に代入してみて、 一番変化の大きい\( y_{腐朽箇所iのE} \)が、腐朽している可能性が高いみたいな推定はできないだろうか。

まずは簡単なモデルで試してみるとして、長方形断面の単純梁に100箇所の腐朽箇所を作って、上記のことをやってみるとか。 上の例は、目的変数が1つの重回帰分析を\( n \)箇所の腐朽箇所に対してやるということだが、 複数の目的変数に対する重回帰分析も調査すべきか。多目的最適化問題も調査すべきか。

• 少なくともCalculixならDakotaを使える

LibreOfficeで対数回帰†

以下の方法で簡単に対数回帰の決定係数も求められそう。但し、LibreOfficeの決定係数の定義がどれを使っているのかとうのは、調べて確認しておく必要がある。

• x,yデータを2列に書いて、マウスで選択
• 挿入→グラフ→散布図→完了でひとまずグラフを描く
• プロットの1つをクリックして、プロットが水色に変わったら右クリック
• 近似曲線を挿入→タイプで線形とか対数とか関数を選択、決定係数にチェック→OK

LibreOfficeの決定係数†

• LibreOffice7.3

重回帰分析†

• 石村貞夫「すぐわかる多変量解析」(東京図書)にのってる例題

  年	医療費の割合(%)\( x_{1} \)	タンパク質摂取量(gf)\( x_{2} \)	平均寿命\( y \)(年)
  1955	3.27	69.7	65.7
  1960	3.06	69.7	67.8
  1965	4.22	71.3	70.3
  1970	4.10	77.6	72.0
  1975	5.26	81.0	74.3
  1980	6.18	78.7	76.2

• \( x_{1}, x_{2} \)を説明変量、\( y \)を目的変量として重回帰分析する。
• 以下のように並べたデータ zyumyou.d (k2のprogに入れてある)を作る。

3
6
3.27, 69.7, 65.7
3.06, 69.7, 67.8
4.22, 71.3, 70.3
4.10, 77.6, 72.0
5.26, 81.0, 74.3
6.18, 78.7, 76.2

• 1行目の3は、説明変量2こと説明変量1この合計
• 2行目の6はデータ数(1955〜1980年まで6こ)
• 3行目は、各行に説明変量、右端に目的変量を並べる

• FORTRAN77で書かれたzyuukait.f (k2のprogに入れてある)を

  gfortran -o zyuukait zyuukait.f

  でコンパイルする。たくさんウォーニングが出るが、実行ファイル zyuukait ができる。

  ./zyuukait<zyumyou.d

  を実行。

0SAMPLE SIZE             =    6
NUMBER OF VARIABLES     =    3      NUMBER OF VARIABLES DELETED =  0
SUM OF SQUARES ATTRIBUTABLE TO REGRESSION   =              73.33898320
SUM OF SQUARES OF RESIDUAL FROM REGRESSION  =               4.39601680
MULTIPLE CORRELATION COEFFICIENT          R =               0.97131286
SQUARE OF   M.C.C ( R)                   R2 =               0.94344868

• Rは重相関係数：「すぐわかる多変量解析」での計算値は 0.9714
• Rは決定係数\( R^{2} \): 「すぐわかる多変量解析」での計算値は 0.9437

VARIANCE OF ESTIMATE                        =               1.46533893
STANDARD ERROR OF ESTIMATE                  =               1.21051185

ziyuudo tyousei R2^=        0.905747804

• ziyuudo tyousei R2^は自由度調整済決定係数\( \hat{R}^{2} \): 「すぐわかる多変量解析」での計算値は 0.9062

0     ANALYSIS OF VARIANCE FOR MULTIPLE  LINEAR REGRESSION
0     SOURCE OF      DEGREES OF    SUM OF       MEAN OF         F
     VARIATION         FREEDOM   SQUARES       SQUARES         VALUE
0 DUE TO REGRESSION  .....    2  73.33898      36.66949      25.02458    
 RESIDUAL ABOUT REG. ....    3  4.396017      1.465339    
 TOTAL    ...............    5  77.73500    
0 VARIABLE NAME              MEAN      STANDARD DEVIATION       
 x1            4.34833      1.19054
 x2           74.66667      5.01305
 x3           71.05000      3.94297
 t(5%)=   3.1819999999999999     

• t値

 b0=   39.290425065219665     
hensuu           bi           gosa         t           t/t(5%)
 x1            2.07682      0.80467      2.58094      0.81111
 x2            0.30440      0.19110      1.59290      0.50060

• b0とbi: \( y=b_{1}x_{1}+b_{2}x_{2}+...+b_{0} \)の定数項\( b_{0} \)と係数\( b_{i} \)
• つまり、\( y=2.07682x_{1}+0.30440x_{2}+39.290425065219665 \)と回帰式が求まったということ。
• 「すぐわかる多変量解析」で求められている回帰式は\( y=2.083x_{1}+0.303x_{2}+39.368 \)なので、どちらか(または両方)の計算にそれなりの誤差があることになる。

1           TABLE OF RESIDUALS 
0  NO. VALUE(x3      )        ESTIMATE        RESIDUAL NORMAL RESIDUAL
    1        65.70000        67.29862        -1.59862        -1.32062
    2        67.80000        66.86249         0.93751         0.77447
    3        70.30000        69.75865         0.54135         0.44721
    4        72.00000        71.42718         0.57282         0.47321
    5        74.30000        74.87126        -0.57126        -0.47192
    6        76.20000        76.08180         0.11820         0.09764

• zyuukait.f は、後藤が大学院生ぐらいの時代の頃(つまり1990年代)の東北大の情報処理センターのFORTRAN77の数値計算ライブラリーのサブルーチンTMRAL(0054F ; COMPUTER CENTER TOHOKU UNIVERSITY)を使うためのプログラム。

         implicit real*8(a-h,o-z)
        CHARACTER*8 vname
        DIMENSION X(100,10),B(10),SB(10),DY(100),A(10,10),XMEAN(10),ANS(9)
       &   ,XSTD(10),vname(10)
        data vname /'x1','x2','x3','x4','x5','x6','x7','x8','x9','y'/
          read(*,*) n
          read(*,*) m
          m1=1
          mm=m
          do 10 i=1,m
           read(*,*) (x(i,j),j=1,n)
           write(*,*) (x(i,j),j=1,n)
  10       continue
          iprint=1
        call  TMRAL(M1,M,N,X,VNAME,MM,IPRINT,
       &               B0,B,SB,DY,ANS,XMEAN,XSTD,A,INDER)
         end

最新の100件