google-site-verification: google3bd66dd162ef54c7.html

Rコマンダーで気圧センサ(MPL115A2)の換算係数を算出

 前回の記事の続き、今回は温度/気圧ロガーの気圧センサ(MPL115A2)の圧力計算用の補正係数を統計手法を使って求めてみました。

◆まえがき
 まえがきと言うより、これからやろうとしていることのおさらいです。
 MPL115A2は、圧力と温度のDAコンバーターの値、PadcとTadcを使って以下の式で気圧を求めます。

 Pcomp = a0 + (b1 + c11*Padc + c12*Tadc)*Padc + (b2 + c22*Tadc)*Tadc --- (1)

 ここでa0などの係数はセンサーのメーカーのフリースケールからROMに書いた値で提供されています。
 これで結構正確な値が出ており、下のグラフのように気象庁の観測値と比べてもほとんど差はありません。

▼測定結果の比較
実測値と気象庁の値の時系列
 赤が(1)式で求めた測定値、黒が気象庁の観測値。
 気象庁のデータはここから入手。10分毎の観測値まであります。

 せっかくこういうデータがあるなら、測定結果から(1)式が逆算できるはずと考えました。

◆作戦は
 (1)式を一般化すると以下の形になります。この式は大量のデータを回帰分析(重回帰分析)することで係数のa~fを推定することができます。気象庁の値は10分毎に記録されており、一日で144件のデータが手に入ります。ということは数日分あれば1000件以上のデータが使えるので勝算がありそうです。

 気象庁の値 = a + b*Tadc + c*Padc + d*Tadc*Padc + e*Tadc^2 + f*Padc^2 ---(2)

◆ツールはRコマンダー
 回帰分析は統計の領域。エクセルでもプラグインで出来ますが、あまり洗練された仕上がりになっていません。そこで、今、流行っている「R」のライブラリで提供されている「Rコマンダー」を使うことにします。
 RについてはWebにいろいろ解説があるので、そちらを参照してください。例えばこちらとか
 というか私、Rについて語れるほどのスキルは持ち合わせていないです。

◆データーの分析
 以下、順に説明していきます

▼エクセルでデーターを準備
EXCELのデータ
・今回は9日分あるので、1296件のデーターがあります。
・分析精度を上げるためTadcとPadcの値は100回測定した平均値を使用。この値はデータロガーのプログラムで自動記録してます。
・Tadc2とPadc2は二乗した値をエクセルで計算(こんなことしなくてもRコマンダー側で対応可能かも)
・右端のkishouchouが気象庁のデータ。

 このデータをクリップボードに乗せて、Rコマンダーにインポートするとデータの準備完了。

▼Rコマンダーのメニュー選択
線形モデルの選択
 統計量、モデルへの適合、線形モデルを選択。

▼線形モデルの計算式を指定
線形モデル計算式の指定
 ここが山場。(2)式を画面のように指定します。前の説明で省略しましたが、Padc*Tadcは交互作用、Padc2やTadc2は二乗の成分になります。
 条件を設定してOKを押すと以下のような結果が出てきます。

Call:
lm(formula = Kishouchou ~ Padc + Tadc + Padc2 + Tadc2 + Padc *
Tadc, data = Dataset)

Residuals:
Min 1Q Median 3Q Max
-1.00512437 -0.19162741 0.00678466 0.19479549 1.05137462

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.221591e+03 6.193085e+02 5.20192 2.2923e-07 ***
Padc -5.323214e+00 2.138911e+00 -2.48875 0.01294476 *
Tadc -4.069904e+00 1.083385e+00 -3.75666 0.00017988 ***
Padc2 2.944575e-03 1.852555e-03 1.58947 0.11220025
Tadc2 2.175988e-03 4.873087e-04 4.46532 8.6944e-06 ***
Padc:Tadc 4.427681e-03 1.856895e-03 2.38445 0.01724843 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3076275 on 1289 degrees of freedom
Multiple R-squared: 0.9873965, Adjusted R-squared: 0.9873476
F-statistic: 20196.78 on 5 and 1289 DF, p-value: < 2.2204e-16

 この結果のEstimateの値が(2)式のa~fの値に相当します。
 Padc2の危険率が0.112と少し高めなのが気掛かりですが、他の因子はびしっと分離できてます。

 具体的には、気圧(hPa) =
 3221.591 - 4.069904*Tadc - 5.323214*Padc + 0.004427681*Tadc*Padc
 + 0.002175988*Tadc^2 + 0.002944575*Padc^2 ---(3)

 というのが知りたかった式になります。

 ここで注意は、表示の桁数を10桁くらい表示させておくこと。デフォルトだと7桁表示なので、桁落ちで充分な計算精度が出ません。やり方は、Rコンソールから options(digits=10)と入力。

◆結論(結果の比較)
相関(オリジナル係数)相関(新係数)

 どちらのグラフも横軸が気象庁の値で、縦軸が気圧の計算結果。左がオリジナルの係数、右が回帰分析で求めた係数を使って計算しています。

 あちゃー、大掛かりな解析計算やったのに両者に大差は無いです。もっと相関が良くなるかと思ってたのですが、_| ̄|○
 まあ、ともかく観察されたデータから変換式を割り出すことは出来たので、一応成功ということにします。

 大量のデータを使って変換式を割り出してもあまり結果が変わらなかった件、これはどういうことなんだろうと考えると、フリースケールが指定してくる係数はけっこう精度が高いということなんでしょう。
 このMPL115A2というセンサーは、同じ秋月で売っているSCP1000と比べるとバラツキが大きいとか精度が悪いとか言われてます。でも充分にアベレージングを行うと精度は悪くないような気がします。

◆おまけ
 せっかくRコマンダーを使ったので機能をちょっと紹介します。
クロス相関図
 散布図行列です。こういうのが一発で出てきます。

3D分布
 3Dグラフ(Padc/Tadc/kishouchou)。これ、実際の画面ではマウスでグリグリ動かせます。
 それにしてもかなりの急斜面で、ちょっとした値の違いで気圧が大きく変わることが判ります。

 この記事ではいきなり回帰分析を始めてます。でも、まずはこんなふうにデータの全体像の把握から始めないとヤバイです。

 Rの本体にはGUIは無いし統計の知識が要求されます。でもRコマンダーはGUI画面でとっつきやすいので、エンジニアなら知っておいて損は無いツールだと思います。何しろタダですから。
関連記事

tag : 気圧センサー 換算係数 重回帰分析 統計 気象庁 散布図

コメントの投稿

管理者にだけ表示を許可する

No title

すごいですね〜。やり方をそっくりそのまま、いただきます(..;)ぺこり

大学にいた頃は、統計計算のデータを端末のPC-98で作って、転送。SPSSで処理するのにメインフレームを数秒シェアさせてもらうために、書類を何枚も書きましたね。
いまじゃ、普通にパソコンで、3D表示で分布を視覚的に観察できる。
良い時代になりました。

Sunday Gamer さん、今晩は

コメントありがとうございます。これが何か応用のヒントになったら幸いです。

SPSSですか、値段が高いので名前を知ってるだけで使ったことがありません。
たしか今年になってSPSSからRの関数を呼び出す仕掛けが出来たような気がします。SPSSにとって自殺行為のような気がしますが、それだけ無視できない存在になっているんでしょうね。
カレンダー
05 | 2017/06 | 07
- - - - 1 2 3
4 5 6 7 8 9 10
11 12 13 14 15 16 17
18 19 20 21 22 23 24
25 26 27 28 29 30 -
プロフィール

ラジオペンチ

Author:ラジオペンチ
電子工作を中心としたブログです。たまに近所(東京都稲城市)の話題など。60過ぎて視力や器用さの衰えを感じつつ日々挑戦!
コメントを入れる時にメールアドレスの記入は不要です。なお、非公開コメントは受け付けていません。

記事が気に入ったらクリックを!
最新記事
カテゴリ
最新コメント
リンク
FC2カウンター
検索フォーム
月別アーカイブ
RSSリンクの表示
QRコード
QRコード