電子計算機を用いた実験式の計算法
岡崎重光
（昭和46年9月13日受理）Calculation of Empirical Formulas with an
Electronic Computer
ShigemituOKAZAKI
Synopsis The errors of experimental data， which have two variables， are studied． It is indic一 ated that the most suitable factors for all measured data can be calculated， if empirica1 formulas are expressed by the following equation． F（の＝α9ω十bh（りじ）十ck（め α， b， c： factors Next equation is one of judgment functions for the empirical formula． n 5一Σ（y，。p−Ycal。）2 r＝1 Yexp：measured data Ycal。：calculated value from empirical equation The value of S is approximated the following equation． n 5一Σ（F（y。xp）−F（y，。1。））29r 7＝1 1 9’ ”（F’（y，alc））・ If・9． i・assum・d t・be c・n・t・nt t・・， b and c， th・n・xt・equ・ti・n・an b…1・・d and the factors are obtained so as to minimize S．嘉一・，芸一α芸一・
1．緒
言 筆者は，四塩化ケイ素，トリクロルシランの蒸気圧 測定，その他の熱化学データ測定を行ない，電子計算 機を用いて，実験データの整理を行なった。蒸気圧の 測定は精度が良いので，電子計算機を用いることが必 要となる。 この方面の専門家による論文を読んでいない。数学 の専門書を見たが，ある関数を簡単な関数で近似する 方法が述べてあった。これは，真の値がわかっている 場合で，測定については少し事情が異なる。 化学系の本には1）2），反応速度と温度，蒸気圧と温度 のように，指数あるいは対数の含まれる実験式に対し て，最小二乗法が示されている。しかし，測定値と実 験式の値の差の二乗和が最小となるよう，計算する方 法ではない。測定値の対数と実験式の値の対数の差の 二乗和が最小となるように計算されている。データの 範囲が狭い場合，だいたい合う。しかし，データの精 度が良く，範囲が広い場合，全く無意味な計算を行な うことになる。このような場合について，最小二乗法 の計算方法を，5と6で述べる。 従来の最小二乗法は，yの誤差のみ考慮している。そこで，xの誤差とyの誤差の両方について検討する 方法を試みた。従来の最小二乗法と基本的な考えは同 じである。直線については同じ結果となる。一般の曲 線についても，傾きの変化が大きくない限り，ほとん ど差はない。 2． 誤差に対する一般的検討 この場合の誤差とは，バラツキであって，計器自身 のくるいなどは含まれない。 実験式 y一ゾ（x） （1） 測定点 （x。，Y。） （夕一1，……， n） xの測定誤差 ±α （α＞0） yの測定誤差 ±β （β＞0） α，βは，κおよびyのバラツキの大きさを示す。α， βが測定範囲により変化することも考えられる。ここ では，一定として取り扱った。そうすれぽ，以下に示 すようにα／βが問題となるのであって，α，βの値 はわからなくてもよい。 真の点を通る実験式を得るのが目的である。真の点 は不明であるから，測定点より真の点を推定する。そ の方法の一つとして，最小二乗法がある。得られる実 験式を推定点の集まりと考える。測定点と測定点に相 当する推定点の差を誤差とする。誤差は性質によって 取り扱いも異なるべきであるが，単なるバラツキと考 え，次の仮定で取り扱った。 数学的に厳密な誤差の取り扱いを知らない。以下の 取り扱いも筆者の独断が多いので数学的には不十分で あろう。 仮定：x，y，それぞれの誤差が±1となるような 平面上で，測定点と推定点の距離の二乗の和（次の式 のS）が最小となるように実験式を定める。 ns・一 Z d2r r＝＝L n：測定点の数 d，：夕番目の測定点と推定点の距離 （2） 測定点に対応する推定点を求める方法として，次の2 つを考えた。 （1） を取る。 これは， うした垂線の足となる。 x＝＝αX s Y＝βY この座標変換により，誤差の範囲をX， 推定点として，測定点に最も近い実験式上の点 図一1のように，測定点より実験式上にお （3） （4） γについて±1 とする。そしてXY平面で垂線の長さを求めれぽ，次 （x。，Y，） 1 Y・『 Gf（α Xr） 図一1XY平面上の垂線 Fig．1 A perpendicular on XY plane y，一 ii−f （・ Fr） 図一2XY平面上での誤差の配分 Fig．2 Distribution of error on XY PIane の式となる。 d…｛Y− 汲?iαx）｝・ 1 1 1＋ ｛ev／βft（αX）｝2 d・−奄堰oy−f（・・）｝・、＋｛α／ts，（x）｝・ （5） （6） （2）Xに対する誤差，Yの誤差が，1：1の比になる よう配分し，推定点を求める。 これを図一2に示す。実際は，誤差が等しい割合で寄 与しているとは限らないが，多数のデータに対する平 均として，この考えは正しいと思う。 ly一獅?iαx）1一醜／βft（αx）l d2＝2鳶2 式（3），（4），（7），（8）より，次の式が得られる。 d2一¥｛y一ゾ（・）｝2｛1＋α／毒（x）1｝・ （7） （8） （9） 式（6）と式（9）では，普通，大きな差はない。筆者は式⑨ を適当と考えるので，以後，式（9）を使う。式（2）に式（9） を代入すれば，次の式が得られる。
S’一w｛Y・−f（・r）｝・｛1＋α／tw1 ’（，。）1｝2 cm α，βは測定点に関係しない定数と考えているから， 次の式が最小になるよう計算すれぽよい。 S一m一f（Xr）｝・石＋α／βト㎞）1｝・⑪ 3． 実験式が直線の場合 f（x）−ax＋b （12） とおけば，次の式が得られる。 S−一
i1＋α1β1。D、塾一…−b）2 ⑬
Sは常に正または0であり，a， bに関して連続関数 であるから，最小点において，次の条件が満足されな けれぽならない。器一・ 器一・ ⑭
これを満足するα，bの組が1つならぽ，それが最小 点となる。式⑬を式⑭の条件で解くことはむずかし い。しかし，データの精度が十分に大きい場合，傾き はほぼ定まったものと考えられるから，1／（1＋α／β 1・D2は後の項に比べて変化カミ少ない・すなわち・近似 的に定数と見ることができる。定数と考えて，式⑬を 式⑭に代入すれば，次の式となる。 7t n ：1 Σ卿。一αΣ況。2−bΣκ。−o 聞 プ＝1 r＝＝1 7＝1 n n ， Σy。−aΣx．−b・n−0 （16） r＝・1 7＝1 n Σ（x。一め（Yr−7） r＝・1 a＝＝ 7t Σ（x。−hi）2 r＝1 b＝＝ynta・x 連立方程式⑬，⑯を解けぽ，a， bが求められる。こ れが直線の最小二乗法である。変形して，式聞の下の 式とすれば，一部，精度の低い計算で足りる。 次の4で述べるように，各データに対する重みg，を 考える。 n s一Σ9r（Y。一αXr−6）2 ⑰ rニ1 この値を最小とするα，bは，次の式より求められる。 n n n αΣ9。x，2＋δΣ9。Xr一Σ 9。x。y． r＝1 r＝1 r＝1 カ 2t n αΣ9。x。＋6Σ9。一Σ9。y． r＝1 t＝1 r＝1 4． 実験式が3個の係数を持つ多項式 ㈹ ⑲ 一般的に，次の式⑳で表わされる場合について，検 討する。式⑳は，その特別な場合である。 ゾ（x）−ax＋bx’＋cx’t （x，x’， x” e＊ xの任意の関数） ゾ（の一α沈2＋6κ＋c 式20）を式（11）に代入すれば，次の式となる。 s一Σ｛y。−ax，−bX。’−cxr”｝29r r＝・11
9「 ＝一｛1＋を／β1fl（脇5肝 grは一般に， a， b， cの関数である。 に関するものであるから， f’（x）が急に変化しなければ， ⑳ ⑳ 22） 23） しかし，傾き データの精度に比べて， 定数と見なしてよい。 すなわち，grは測定点それぞれに固有の値と見なす。 そうすれぽ，次の式の計算が容易となる。1／一・・器一・・票一・ ㈱
式⑭に，式22）を代入すれば，次の連立方程式となる。 n n タど αΣ9。x。2＋δΣ9。蹴／＋εΣ9。り6，バ r＝1 ，＝1 r＝1n
一Σ 9，X，Y， ㈲ r＝1 71 n 7t aΣ9。x。tりじ。＋bΣ9。κ。’2＋c Z 9，X，’：’：，” r＝1 ア＝1 7： 一Σ9。c。’yγ r＝1 Vl n 7＝1 rニ1 71 一Σ9。x，”Yr r＝1 これを解けば，α，b， r＝・1 aΣ9。バκ。＋b ： 9，Xr’W＋c X 9，πr’〃2n r＝1 26） ￠7） cが求められる。grはf（X） より，式囲に従って計算される。しかし，最初1・Xf（x） が決まっていない。そこで，grを適当な値，場合によ ってはg。＝1として，式2S）∼（2Zを計算し， a， b， cを求める。その値を用いてgrを計算し，再び連立 方程式を解けぽよい。普通の場合，二度計算すれば十 分である。f’（x）があまり変化しなけれぽ， g。−1で もよい。 なお，式2S）∼tl 7）の計算はもちろん，式⑬，（16）の計算 も，十分な精度で行なわなけれぽならない。測定値が 3桁の精度であれぽ，実験式．の係数は4桁程度求め る。式⑭また式⑳からわかるように，Sについて4桁 の精度を必要とし，測定値よりSを計算するとき7桁 以上の精度が必要となる。したがって，8桁以上の精 度でなけれぽ，信頼できない。以上は，指数，対数が 含まれない場合についての考察である。 5． 実験式が特別な場合（蒸気圧と温度の関係） 蒸気圧と温度の関係を表わす式はいくつかあるが，最も簡単な式は次の式である。
lnP−a／T−F 6 ⑳
（P：蒸気圧，T：絶対温度） この場合y＝lnP， x＝1／Tとおけぽ，3に示した最 小二乗法を用いることができると，一般の本には書か れている。しかし測定誤差について最小二乗法となっ ていない。特に広い温度範囲で測定したデータの場 合，全く役に立たない。 ここでは計算機を用いるので，さらに複雑な次の式に ついて考える。 ln1）＝＝a／T十1．751nT十bT−1−e ￠9） これをNernstの式という。この場合も，次のように 置けば，4で示した計算を行なうことができる。 y＝ln1）−1．751nT ｛ an−1／T （30） xt＝T x’t＝1 しかし，低温ではPが小さくなるため，1nPの誤差は 大きくなる。そこで，誤差の変化をおぎなうため，デ ータに対する重みとして蒸気圧を使った。ただし，式 倒の値は1とした。 lt S＝ΣPr（ln（1）r／760）一α／Tr−1．751nTr− bTr−c）2 r＝1 圃 式倒が最小となるよう計算した結果を，図一3のAに 示す。一方，24°C以上の高温側の実験データについ て，式鋤が最小となるよう計算した結果を，図一3の Bに示す。全データの誤差は，AよりBが小さい。す なわち，次の式の値を最小とするa，b， Cは，式⑪ について計算したのでは，不十分であることがわか る。 S＝Σ｛Pr−760exp（a／Tr十1．751nTγ十bT，十c）｝2 r＝1 働 最初，式（32）’vcついて，4に示したように直接計算す る方法は，わからなかった。そこで，計算機により多 数の方法を試みた。四塩化ケイ素の蒸気圧測定データ について，計算した結果を表1に示す。川∼（V）の場 合，最初，式β蜘こついて計算し，a， b， cを求めて おく。 川 a，b， cを少し変化させて，式㈱の値を計算 し，その値が最小となる値をさがす。変化量は最初大 きく，次第に小さくする。プログラムを工夫した。表 1に示すように，Sの減少は最初，比較的大きいが， ．終わりになると非常に小さくなる。 川 式閲をa，b， cの関数と見なし，式⑭を満足 表一1 各方法による計算結果 Table l Results of calculation by various methods method （1） （1） （ll） （川） （lv） （v） （Vi） computer time 7min 22 4 2 3 1 0 30sec 30 0 0 0 50 10 Sof equation （32） start end 597．1763 597．1763 597．1763 597．1763 597．1763 597．1763 440．771 440．4510 507．8584 585．2848 464．3581 586．8856 87． 2 するα，b， Cを次の方法により求める。 a・一・・−1−［㌃／箒1。＿＿、．。。＿．、㈹b，−bt−一［怨／莞］。＿＿。。＿（・4
c・一・・一・一際／莞］a＝＝a、，．、．b、，一。−ct−、倒 これらをα，b， cの順にくり返す。ニュートンの 反復法の応用である。 刷 次のようにおく。 α＝αz，b・＝・βt， c一γt β6） （α2十β2＋γ2−1） そして，次の式が最大となるよう，α，β，γを定 める。怨一器α＋器β＋怨・ （・力
次の計算をくり返して，α，b， cを求める。At一莞／箒 （・8
il≡il聞 ㈲
やはりニュートンの反復法の応用である。 （1v） （lli）では∂s／∂a，∂s／∂b，∂s／∂cの値に大きな差 があるため，十分でない。そこで，次の変換を行なっ た。 α一a。at， b−b。b’， C＝：coc’ ㈹ ao， bo， CO，は最初の値である。したがって， at， bt， ctは最初1である。α’， b’， ctについて，（ill）の 方法で計算する。 （v）（jV）と同じような方法であるが，4tを式㈱より 求めず，0、1程度から次第に小さくする。そして，式 働が最小となるa，b， Cをさがす。川と（iv）を混合し た方法である。以上5種の方法を試みた。同一の実験データについ て計算した結果を表一1に示す。使用した計算機は， FACOM 270−30である。 次に述べる（Vi）を除くと，良い方法は川または（Mであ る。しかし，どの方法も次第に収束は遅くなる。図一3 のAとBより，式B2）の値は100以下になりうることが わかっていた。プログラムの工夫，あるいは計算精度 が不足していたかもしれない。しかし，根本的に考え なければ，以上の方法の改良では期待できない。複雑 でも，くり返し演算でないことが望ましい。式助を式 佗4に代入して解く方法は考えられないので，データの 重みgrを利用する方法を再び考えた。今度は理論的に g，を決定した。 M）式？3）に示したgrを利用する方法
努一き 91）
であるから，∠1P≒P41nP ㈹
と書ける。実験値をP，。p，計算値をP，。1。で表わせば， 誤差の小さいとき，次の式が成り立つ。1：；塩慧㌦｝
式⑪はこの場合，次の式となる。 s一c（P・xp−Pca・c）2a＋α／β1‘P〃τ1）・ ⑬ ⑭ 式泌，⑬を用いて，式（44を書き変えれぽ，次の式が得 られる。 n S＝ｶ（1・P・x・−1・Pcal・）2’（、＋α／βldP／dTI）・ 1）2 ks） これならぽ，次のように置いて，式囲∼伽を計算す ることができる。 ・一i1＋α晶4TD・ 陶
最初，dP／dTは不明である。そこで，式⑰のよう に決め，式⑱が最小となるa，b， cを求める。その 値を用いて式鯛を計算する。9−P。xp2 ㈲
n S一ΣP。xp2（lnP，。r lnP，、1，）2 9S） 7＝1 式㈲として計算した結果を図一3トのCに，式陶にっい て計算した結果を，図一3のDに示す。 なお，根本的に同じであるが，別の考えで式㈲を求 めることができる。式㈹による変換を行なうと共に， そのときβが式⑫に従って変化すると考え，式⑩を計 算する。そうすれば，式㈲が得られる。 十3 一3 （mmHg） 十2 0 −tP’一”“”…．．．．」 一一Q 言 L ＋2 己 8 0 tL 9−2 口→2 ξ ‖。 烏 一2 一50℃ 0℃ 50℃ （C） @ 0℃ 1 一50℃ @ （D） @ 0℃ 50℃ @ 」 一50℃ @ Tcmperatllre 50ec1
（A）全データ（データ数123）について，式⑳の値が最小と なるように計算 Sof equation（31）is minimized on all data （123points） （B）24°C以上のデータ（データ数63）にっいて，式6Dの値 が最小となるよう計算 Sof equation（31）is minimized on the data higher than 24°C （63 points） （C）全データについて，式㈹の値が最小となるように計算 Sof equation（48）is minimized on all data． （D）全データについて，（C）の結果を用い，式㈲の値が最小 となるよう計算 Sof equation（45）is minimized on all data， using the result of（C） 全データに対する式㈹の値 S of equation （32） on all data （123 points） （A） 597．2 （B） 97．9 （C） 87．2 （D） 92．5 図一3各方法によるP。xp−P。al。の値 Eig．3 Pexp−Pca1。 calculated by above methods 6． 5に示した方法の一般化 式⑪が直接に，4の方法で解けない場合，その工夫 を5の⑭に示した。一般に，式（1）が次のような形式に 変換できる場合，以上の方法が可能である。 F（y）＝＝ag（x）十bh（x）十cle（x） ㈲ a，b， c：係数 c−1となっても，項の数がさらに多くても可能であ る。この場合，式（11）は5に示したと同様に次の式で近 似される。lt s一Σ｛．F（Yr）−F（プ（x，））｝29， 60） r＝1 9r− o、＋α／β1∫d｝z｛F，（f（Xr））｝・ 61） g，をa，b， cに対して一定と仮定すれば，4に示し た方法を用いることができる。測定精度が十分に良け れば，この仮定が可能である。 7． 結 言 以上，電子計算機を用いたデータ整理法について， 誤差の評価，最小二自法の改良について述べた。 以上のほか，多数の測定を行なうと，目盛の読み違 い，操作上の誤りなどが起こる。そして簡単にわから ないものもある。一定の基準を与えて，計算機により 選別することが必要である。これはデータによるが， 基準を考えるためには，2で述べたような誤差の評価 が必要である。 実験を行なえば，常に確率が含まれる。専門家から 見れぽ不十分な計算，自分の期待する解釈を行なうこ ともある。計算機を利用し，誤差についていろいろな 面から検討を加えれば，限られたデータから，さらに 良い結果を得ることができる。 以上は，東京大学生産技術研究所野崎研究室で研究 中に，主として考えたことである。 文 献 1）神野博：化学 第26巻（第3号）p．250（1971） 2） 清水祥一：化学と電子計算機 p．18 南江堂 （1971）
