通常用いられる最小２乗法は，計算を単純にするため「縦軸のデータだけに誤差がある場合で，横軸のほうの誤差は全くないかあるいは無視できるものと見なすことができる」という前提をおいて，目的関数を

　　ｓ＝Σｗi（ｙi−ｆ（ｘi））^2　　　（ｗi：重み）

あるいは

　　ｓ＝Σｗi（ｚi−ｆ（ｘi，ｙi））^2 　　　（ｗi：重み）

とし，誤差はあくまでもｙのみに含まれるものとしています．それとは逆に，ｘのみに誤差がある場合は従属変数と独立変数を入れ替えるだけで，この解析法がそのまま適用できます．

　

　プログラム・パッケージに入っている最小２乗法は，通常，ｙのみに誤差があることを前提とした最小２乗法ですから，ｘのみに誤差がある場合やｘ，ｙ双方に誤差があり，その比があまり１から離れない範囲の場合はそのまま用いると誤った結論に陥ってしまうことになりかねません．実際問題としてｘ，ｙのどちら側にも誤差があるほうが当り前ですから，今回のコラムではｘ，ｙ双方に誤差があるときの最小２乗法の取り扱いについて考えてみます．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

　ｘ側，ｙ側どちらにも誤差がある場合は独立変数，従属変数の区別がなく対等と考えられますから，ｆ（ｘ，ｙ）＝０のように陰関数表現しておく方が楕円のような多価関数となるデータに対しても曲線のあてはめができるので何かと都合がよくなります．そのような方法に，青山学院大学理工学部・物理を定年退官された粟屋隆先生による２次元フィッティングの方法があります．

　

　粟屋の方法は，品質管理・計量管理で有名なデミング（Deming）のクラシカルな方法を改良したものですが，デミングの方法が反復計算にそぐわない１回限りの計算法であるのに対し，粟屋の方法は反復計算によって精度を逐次高めていくことが可能になっています．

　

　この２次元フィッティングは実に念の入った膨大なものであり，未知数が多くなって計算はかなり複雑になります．この解法については

　　粟屋隆著「データ解析」（学会出版センター）

に詳しいのですが，特殊な問題にはこのような考察が必要になり，拡張した取り扱いを然るべく行なえば解が得られることをアウトラインだけでも知っておいたほうがよいと思われます．

　

　そのすべてをここに記すことは量的に無理であり，また小生の能力を超えたところも少なくないので，詳細は省略し理論の要点のみを記しておきます．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

［１］ラグランジュの未定係数法

　

　粟屋の方法ではラグランジュの未定係数法が使われていますが，ラグランジュの未定係数法は，たとえば，「ｘ^3−３ｘｙ＋ｙ^3＝０の条件のもとでｘ^2＋ｙ^2の極値を求めよ」といった条件つき極値問題や制約条件付きの最小２乗法などの解を得るために導入された方法です．

　

　制約条件付きの最小２乗法は１９世紀の測量士たちのおこなった三角測量でもっともよく用いられましたが，以下にその簡単な例を示します．

　「三角形の三つの内角を等精度で測定し，α＝５４°０５′，β＝５０°０１′，γ＝７６°０６′を得た．内角の和は２直角になるべきであるが，測定誤差のため１８０°１２になった．」

　

　このような場合に内角の最確値ｘ，ｙ，ｚおよび確率誤差を求めるにはどうしたらよいかを考えてみましょう．ｘ＋ｙ＋ｚ＝πが要請されている条件です．また，測定精度は等しいので荷重をｗi ＝１とします．

　

　このような制約条件付きの問題では，変数λを新たに導入して，目的関数を

　　ｓ＝Σ（測定誤差）^2−λ（条件式）

　　　＝（ｘ−α）^2＋（ｙ−β）^2＋（ｚ−γ）^2−λ（ｘ＋ｙ＋ｚ−π）

としてｓを最小にするｘ，ｙ，ｚおよびλを求めます．極値の必要条件により，　　∂ｓ／∂ｘ＝２（ｘ−α）−λ＝０・・・・・（１）

　　∂ｓ／∂ｙ＝２（ｙ−β）−λ＝０・・・・・（２）

　　∂ｓ／∂ｚ＝２（ｚ−γ）−λ＝０・・・・・（３）

　　∂ｓ／∂λ＝−（ｘ＋ｙ＋ｚ−π）＝０・・・（４）

前の３つの式を加え合わせて，

　　２（ｘ＋ｙ＋ｚ−α−β−γ）−３λ＝０

が得られます．（４）式より，ｘ＋ｙ＋ｚ＝πを代入すると

　　λ＝２／３（π−α−β−γ）

が得られます．このλを（１），（２），（３）式に代入して

　　ｘ＝α＋１／３（π−α−β−γ）＝５４°０１′

　　ｙ＝β＋１／３（π−α−β−γ）＝４９°５７′

　　ｚ＝γ＋１／３（π−α−β−γ）＝７６°０２′

が求める解となります．すなわち，等精度で三角形の内角を測定したときの各補正量は測定値の和と１８０°との差を３等分したものです．

　

　このように，制約条件のある最小２乗問題を解くには未定係数λを導入して，未知のパラメータと未定係数λについて方程式を解けばよいことになり，制約条件のない場合の手法を用いて容易に解くことができることが理解されるでしょう．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

［２］粟屋の方法

　

　２次元フィッティングでは，取り扱う陰関数をパラメータも含めｆ（ｘ，ｙ，ａ）＝０と書くことにします．ここで，データ（ｘi，ｙi）の真の値を（ｘi0，ｙi0），真の係数を示すベクトルａ0＝（ａ10，・・・，ａm0）とすると，

　　ｆ（ｘi0，ｙi0，ａ0）＝０

が恒等的に成り立ちます．したがって，制約条件式にはｆ（ｘi0，ｙi0，ａ0）＝０を用いることになります．

　

　目的関数ｓは

　　ｓ＝Σ（ｘi−ｘi0）^2／２σxi^2＋Σ（ｙi−ｙi0）^2／２σyi^2

　　　＝ΣΣ（ｘi−ｘi0）^2／２σxi^2＋Σ（ｙi−ｙi0）^2／２σyi^2−Σλiｆ（ｘi0，ｙi0，ａ0）　　　　　　（λi ：ラグランジュの未定係数）

と表わすことができますから，

　　∂ｓ／∂ｘi0＝−（ｘi−ｘi0）／σxi^2−λi∂ｆ（ｘi0，ｙi0，ａ0）／∂ｘi0＝０

　　∂ｓ／∂ｙi0＝−（ｙi−ｙi0）／σyi^2−λi∂ｆ（ｘi0，ｙi0，ａ0）／∂ｙi0＝０

　　∂ｓ／∂ａj0＝−Σλi∂ｆ（ｘi0，ｙi0，ａ0）／∂ａj0＝０

　　　（ｉ＝１〜ｎ，ｊ＝１〜ｍ）

を解けばよいことになります．

　

　ここで解くべき連立方程式の未知数は，ｘi0，ｙi0，λiがｎ個ずつ，ａj0がｍ個で合計３ｎ＋ｍ個となり，未知数が多すぎて計算は大変複雑となります．そこで，まずλiを消去することを考えます．以下，表記を簡単にするため，

　　∂ｆ（ｘi0，ｙi0，ａ0）／∂ｘi0≡〜　ｆxi

　　∂ｆ（ｘi0，ｙi0，ａ0）／∂ｙi0≡〜　ｆyi

　　∂ｆ（ｘi0，ｙi0，ａ0）／∂ａj0≡〜　ｆji

と表わすことにします．

　

　アプロ−チの途中を省略しますが，誤差伝播の法則より，ｚ＝ｆ（ｘ，ｙ）の母分散σz^2は，分布形にかかわらず，

　　σz^2＝（∂ｆ／∂ｘ）^2σx^2＋（∂ｆ／∂ｙ）^2σy^2

で近似することができることより，重み係数に相当するのが

　　ｗi＝｛（σxiｆxi）^2＋（σyiｆyi）^2｝^(-1)

また，２変数関数のテイラー展開より

　　ｆ（ｘ＋Δｘ，ｙ＋Δｙ）＝ｆ（ｘ，ｙ）＋Δｘｆx＋Δｙｆy＋１／２（Δｘ^2ｆxx＋２ΔｘΔｙｆxy＋Δｙ^2ｆyy）＋・・・・・

　

　ニュートン法の漸化式ではｆ（ｘ＋Δｘ，ｙ＋Δｙ）＝０とおいてΔｘ，Δｙを求めるわけですから，残差の近似値に相当するのがｆi0＋（ｘi−ｘi0）ｆxi＋（ｙi−ｙi0）ｆyiになり，結局，

　　Ｆj≡Σｗi｛ｆi0＋（ｘi−ｘi0）ｆxi＋（ｙi−ｙi0）ｆyi｝ｆji＝０　　　（ｊ＝１〜ｍ）

を解けばよいことになります．

　

　これをガウス・ニュートン法と同様の方法で整理することにします．真の値ｘi0，ｙi0，ａ10，・・・，ａm0の近くの点（ｘi，ｙi，ａ1，・・・，ａm）で関数ｆ（ｘi0，ｙi0，ａ0）＝０をでテイラー展開し，１次までの項をとると

　　ｆi0＝ｆi＋ｆxiΔｘi＋ｆyiΔｙi＋ΣｆkiΔａk＝０

となり，また，微分係数は点（ｘi，ｙi，ａ1，・・・，ａm）における値で置き換えると，

　　Σｗi｛ｆi＋（ｘi−ｘi0）ｆxi＋（ｙi−ｙi0）ｆyi｝ｆji＝ΣｗiｆaiΣｆkiΔａk

が得られ，Δａkが計算されることになるのです．

　

　多価関数に対する取り扱い（２次元フィッティング）において，パラメータを計算するには最初に未知のパラメータの初期値を与えなくてはいけないことはガウス・ニュートン法と同じですが，粟屋の方法ではさらに測定値の真の値も与えなければなりません．測定値の真の値の初期値には測定値自身を用いることにし，パラメータと測定値の真の値の第１近似値を与えてやって逐次近似解法により漸近連立方程式を解きます．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

［３］デミングの方法による円近似

　

　粟屋の方法の欠点は計算に時間がかかるということです．未知数が多すぎることがその原因です．そこで，以下には粟屋の方法の原型となったデミングの方法による円近似プログラムを掲げます．

　

　デミングの方法については，

　　本間仁，春日屋伸昌「次元解析・最小２乗法と実験式」コロナ社

に譲りますが，（その４）に掲げた非線形最小２乗法プログラムの粟屋の方法の簡易版と比べて変わっているのは，重みをつけたこと（2290行〜2370行）と正規方程式の右辺の式（2860行）など数カ所を手直ししただけであって，変更点はごくわずかです．

　

　残差２乗和は，測定点（ｘ，ｙ）からあてはめるべき円に下ろした垂線の長さ（データから円までの最短距離）

　　ｓ＝√｛（ｘ−ｘ0）^2 ＋（ｙ−ｙ0）^2｝−ｒ

の平方和Σｓ^2が最小になるように定めていますが，この点は変わっていません（3600行）．

　

　垂線の足を（ｘ1，ｙ1）とすると，

　　ｘ1＝｛ｘｒ＋ａ（√｛（ｘ−ｘ0）^2 ＋（ｙ−ｙ0）^2｝−ｒ）｝／√｛（ｘ−ｘ0）^2 ＋（ｙ−ｙ0）^2｝

　　ｙ1＝｛ｙｒ＋ａ（√｛（ｘ−ｘ0）^2 ＋（ｙ−ｙ0）^2｝−ｒ）｝／√｛（ｘ−ｘ0）^2 ＋（ｙ−ｙ0）^2｝

より，直截的に

　　（ｘ−ｘ1）^2 ＋（ｙ−ｙ1）^2＝ｓ^2

が確かめられますが，円の場合，垂線の足の座標を求めることが楕円に比べてかなり簡単というわけです．

　

　ところが，この計算を実行すると最初は残差２乗和が減少するものの，その後，増加に転じてしまいます．そのため，プログラムでは反復計算を残差２乗和の減少が増加に転じる直前で打ち切っています（2540行〜2570行）．もともとデミングの方法は反復計算法ではなく１回限りの計算法なのですが，このプログラムでは無理矢理，反復計算にもちこんでいるのです．

　

1000 '

1010 ' **** circular approximation ****

1020 '　　　　　　　　　　　　　2003/01/15 (C) ｻﾄｳ ｲｸﾛｳ

1030 SCREEN 3,0:CONSOLE ,,0,1

1040 CLS 3:WIDTH 80,25:COLOR 6

1050 GOSUB *INITIALIZE

1060 GOSUB *READ.DATA

1070　GOSUB *INIT.VAL :L2$="LINEAR METHOD"

1080　GOSUB *EXECUTION:L2$="DEMING METHOD"

1090 GOSUB *SCATTER.DIAGRAM

1100 GOSUB *REPORT

1110 END

1120 '

1130 ' *** 初期設定 ***

1140 '

1150 *INITIALIZE:

1160 M=3:' No. of parameters

1170 DIM P(M),W(M)

1180 DIM A(M)

1190 DIM P0(M),P1(M),DP(M)

1200 DIM AZ(M,M),BZ(M,M)

1210 DIM HESSE(M,M)

1220 DIM SD(M),SE(M)

1230 DIM COV(M,M)

1240 '

1250 NV=2:' No. of independent variables

1260 DIM X(NV)

1270 DIM AVE(NV),DEV(NV),SEM(NV)

1280 DIM XMAX(NV),XMIN(NV)

1290 FOR J=0 TO NV

1300　 XMAX(J)=-1E+10:XMIN(J)=1E+10

1310 NEXT J

1320 '

1330 RX=.7:RY=.8:BOX=1:KIND=1:POWER=0

1340 IX=.7:IY=.8:MARK=1:CLR=2:DRAW.COLOR=4

1350 TXT.COLOR=6:GRP.COLOR=5:REF.COLOR=6

1360 SCALE.KIND=1:PLOT=1

1370 PI=3.14159

1380　JX=.709:JY=1

1390 'JX=1　 :JY=1

1400 RX=IX*JX:RY=IY*JY

1410 '

1420 'FORM$="((X(1)-P(1))/P(3))^2+((X(2)-P(2))/P(3))^2-1"

1430　FORM$="(X(1)-P(1))^2+(X(2)-P(2))^2-P(3)^2"

1440 RETURN

1450 '

1460 ' *** DEFINED EQUATION ***

1470 '

1480 *DEFINE.FORMULA:

1490 'DEF FNY(X)=((X(1)-P(1))/P(3))^2+((X(2)-P(2))/P(3))^2-1

1500　DEF FNY(X)=(X(1)-P(1))^2+(X(2)-P(2))^2-P(3)^2

1510 RETURN

1520 '

1530 *DEFINE.FORMULA2:

1540 ON J GOSUB *J1,*J2,*J3

1550 RETURN

1560 '

1570 '*J1:DEF FNJ(X)=-2*(X(1)-P(1))/P(3)^2:RETURN

1580 '*J2:DEF FNJ(X)=-2*(X(2)-P(2))/P(3)^2:RETURN

1590 '*J3:DEF FNJ(X)=-2*((X(1)-P(1))^2+(X(2)-P(2))^2)/P(3)^3:RETURN

1600　*J1:DEF FNJ(X)=-(X(1)-P(1))*2:RETURN

1610　*J2:DEF FNJ(X)=-(X(2)-P(2))*2:RETURN

1620　*J3:DEF FNJ(X)=-P(3)*2　　　 :RETURN

1630 '

1640 ' *** データ読み込み ***

1650 '

1660 *READ.DATA:

1670 CLS 3

1680 RESTORE *D

1690 READ N1

1700 'DFILE$="DATA01.TXT"

1710 'OPEN DFILE$ FOR INPUT AS #1

1720 'INPUT #1,N1

1730 '

1740 DIM XX1(NV,N1),YH(N1),WT(N1)

1750 DIM Y1(N1),E1(N1)

1760 DIM Z1(N1)

1770 DIM JACOBI(M,N1)

1780 FOR J=1 TO N1:WT(J)=1:NEXT J

1790 '

1800 FOR I=1 TO N1

1810　 FOR J=1 TO NV

1820　　 READ XX1(J,I)

1830　　'INPUT #1,XX1(J,I)

1840　 NEXT J

1850 NEXT I

1860 'CLOSE #1

1870 '

1880 PRINT "ただいまデータの取り込み中"

1890 FOR I=1 TO N1

1900 PRINT "*";

1910 Y1(I)=0

1920 X(0)=Y1(I)

1930 FOR J=1 TO NV

1940　 X(J)=XX1(J,I)

1950 NEXT J

1960 'Z1(I)=1

1970　Z1(I)=X(1)*X(1)+X(2)*X(2)

1980 GOSUB *CALCULATE

1990 NEXT I

2000 '

2010 FOR J=0 TO NV

2020　 SS=DEV(J)-AVE(J)*AVE(J)/N1

2030　 DEV(J)=SQR(SS/(N1-1))

2040　 SEM(J)=DEV(J)/SQR(N1)

2050　 AVE(J)=AVE(J)/N1

2060 NEXT J

2070 RETURN

2080 '

2090 ' *** CALCULATE ***

2100 '

2110 *CALCULATE:

2120 FOR J=0 TO NV

2130　 IF XMAX(J)<X(J) THEN XMAX(J)=X(J)

2140　 IF XMIN(J)>X(J) THEN XMIN(J)=X(J)

2150　 AVE(J)=AVE(J)+X(J)

2160　 DEV(J)=DEV(J)+X(J)*X(J)

2170 NEXT J

2180 RETURN

2190 '

2200 ' *** 非線形最小２乗法 ***

2210 '

2220 *EXECUTION:

2230 PRINT:PRINT "ただいま計算中"

2240 GOSUB *WEIGHT

2250 GOSUB *GAUSS.NEWTON

2260 GOSUB *FINAL

2270 GOSUB *GAUSS.NEWTON.ERROR

2280 RETURN

2290 '

2300 ' *** 重み ***

2310 '

2320 *WEIGHT:

2330 FOR I=1 TO N1

2340　 FOR V=1 TO NV:X(V)=XX1(V,I):NEXT V

2350　 WT(I)=1/((X(1)-P(1))^2*4+(X(2)-P(2))^2*4)

2360 NEXT I

2370 RETURN

2380 '

2390 ' *** GAUSS-NEWTON METHOD ***

2400 '

2410 *GAUSS.NEWTON:

2420　LMAX=50

2430 'LMAX=500

2440 LOOP=0

2450 GOSUB *RESIDUAL:SO=SS

2460 GOSUB *TEMPORARY

2470 '

2480 GOSUB *JACOBIAN

2490 GOSUB *HESSIAN

2500 FOR I=1 TO M:P1(I)=P(I):NEXT I

2510 GOSUB *MARQUARDT

2520 GOSUB *UPDATE

2530 GOSUB *TRANSFORM

2540 IF SN<SO THEN 2580

2550 SN=SO

2560 FOR I=1 TO M:P(I)=P1(I):NEXT I

2570 RETURN

2580 '

2590 IF SN<SO THEN 2610

2600 FOR I=1 TO M:DP(I)=DP(I)/2:NEXT I:GOSUB *TRANSFORM:' [DAMPING]

2610 GOSUB *CONVERGENCE:IF SW=0 THEN RETURN

2620 SO=SN:IF LOOP<LMAX THEN 2480

2630 RETURN

2640 '

2650 ' ** X **

2660 '

2670 *JACOBIAN:

2680 FOR I=1 TO N1

2690　 FOR V=1 TO NV:X(V)=XX1(V,I):NEXT V

2700　 FOR J=1 TO M

2710　　 GOSUB *DEFINE.FORMULA2

2720　　 JACOBI(J,I)=FNJ(X)

2730　 NEXT J

2740 NEXT I

2750 GOSUB *JACOBIAN.SUB

2760 RETURN

2770 '

2780 ' ** Y **

2790 '

2800 *JACOBIAN.SUB:

2810 FOR I=1 TO N1

2820　 FOR V=1 TO NV:X(V)=XX1(V,I):NEXT V

2830　 GOSUB *DEFINE.FORMULA:YH=FNY(X)

2840　'E1(I)=Y1(I)-YH

2850　'E1(I)=SQR((X(1)-P(1))^2+(X(2)-P(2))^2)-ABS(P(3))

2860　 E1(I)=(X(1)-P(1))^2+(X(2)-P(2))^2-P(3)^2

2870　 YH(I)=YH

2880 NEXT I

2890 '

2900 ' ** XT*Y **

2910 '

2920 FOR J=1 TO M

2930　 S=0

2940　 FOR I=1 TO N1

2950　　 S=S+JACOBI(J,I)*E1(I)*WT(I)

2960　 NEXT I

2970　 W(J)=S

2980 NEXT J

2990 RETURN

3000 '

3010 ' ** XT*X **

3020 '

3030 *HESSIAN:

3040 FOR M1=1 TO M

3050　 FOR M2=M1 TO M

3060　　 S=0

3070　　 FOR I=1 TO N1

3080　　　 S=S+JACOBI(M1,I)*JACOBI(M2,I)*WT(I)

3090　　 NEXT I

3100　　 HESSE(M1,M2)=S:HESSE(M2,M1)=S

3110　 NEXT M2

3120 NEXT M1

3130 NZ=M

3140 RETURN

3150 '

3160 ' ** MARQUARDT **

3170 '

3180 *MARQUARDT:

3190 FOR I=1 TO M

3200　 FOR J=1 TO M:AZ(I,J)=HESSE(I,J):NEXT J

3210 NEXT I

3220 RETURN

3230 '

3240 ' ** UPDATE **

3250 '

3260 *UPDATE:

3270 GOSUB *INVERSE:' (XT*X)^(-1)

3280 GOSUB *UPDATE.SUB

3290 RETURN

3300 '

3310 ' ** (XT*X)^(-1)*XT*Y **

3320 '

3330 *UPDATE.SUB:

3340 FOR J=1 TO M

3350　 S=0

3360　 FOR K=1 TO M

3370　　 S=S+BZ(J,K)*W(K)

3380　 NEXT K

3390　 DP(J)=S

3400 NEXT J

3410 RETURN

3420 '

3430 ' ** TRANSFORM of PARAMETERS **

3440 '

3450 *TRANSFORM:

3460 FOR I=1 TO M

3470　 P(I)=P1(I)+DP(I)

3480 NEXT I

3490 GOSUB *RESIDUAL:SN=SS

3500 RETURN

3510 '

3520 ' ** SUM of RESIDUALS **

3530 '

3540 *RESIDUAL:

3550 SS=0

3560 FOR K=1 TO N1

3570　 FOR V=1 TO NV:X(V)=XX1(V,K):NEXT V

3580　 GOSUB *DEFINE.FORMULA:YH=FNY(X)

3590　'SS=SS+WT(K)*(Y1(K)-YH)^2

3600　 E1(K)=SQR((X(1)-P(1))^2+(X(2)-P(2))^2)-ABS(P(3))

3610　'SS=SS+WT(K)*(E1(K))^2

3620　 SS=SS+E1(K)^2

3630 NEXT K

3640 RETURN

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

［４］デミング法と粟屋法の比較

　

　計算速度の点でアドバンテージを期待したデミングの方法による計算結果は

　　ｘ0＝１０５．５７±０．４２

　　ｙ0＝　８０．３２±０．１０

　　ｒ ＝　８０．６７±０．４０

でした．

　

　最小２乗法ソフト「耕太郎」にサポートされている粟屋の方法の最終収束値

　　ｘ0＝１１８．７５±３．１３

　　ｙ0＝　８５．５８±０．３７

　　ｒ ＝　９１．９６±３．１１

と比較して，計算所要時間がはるかに短くて済むわりには，まあまあよい近似値と誤差が求められたといってもよいでしょう．

　

　もちろん，粟屋先生が開発された方法のほうが，収束安定性も高いことから安心して使うことができます．粟屋の手法の根拠は最尤法にあり，実にみごとな計算手法で特別な職人芸的技巧を必要とせず，しかも，made in Japan （純国産）の方法ですから，読者が粟屋の方法に興味をもち使って下さることを願ってここで宣伝しておきます．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝