2018/03/18

Excel分析ツールを見直す(11)

　前回の続きです。今回はRとPythonでフーリエ解析を行います。データは前回のExcelのB、C列(t、y)をcsvファイルにしたものを使いました。

2. R
　フーリエ解析の基本的なアプローチは、Excelで行ったやり方と同じです。プログラムは以下の通りです。

x <- read.csv("FourierData.csv",header = FALSE) #Load data

datanum <- length(x[,1])              #data number
sampintv <- 0.01                      #sampling interval[sec]
funfreq <- 1/(datanum*sampintv)       #fundamantal frequency

anafreq <- funfreq * seq(0,datanum-1) #analysis frequency
res <- fft(x[,1])                     #FFT results
spec <- abs(fft(x[,1]))               #FFT spectra

#Graph output
plot(x[1:128,1],type="l",col=2)       #data plot(128data)
plot(anafreq,spec,type="l",col=2)     #spectra 
#File output
write.csv(res, file = "FFT_res.csv",quote = FALSE)
write.csv(cbind(anafreq,spec), file = "FFT_spec.csv",quote = FALSE)

FFTを実行しているのは8行目です。その結果(FFT_res.csv)は以下の通りで、Excelの結果と同じでした。

,x
1,21.3114464479788+0i
2,21.3198777521124-0.5839506913525i
3,21.3452101541248-1.1692561369683i
･･･
1022,21.3875595649728+1.7572814709631i
1023,21.3452101541249+1.1692561369683i
1024,21.3198777521126+0.5839506913526i

　また、9行目でスペクトルを求めますが、abs関数で虚数の計算も対応できています。Excelでは虚数をabs関数で処理するとエラーを吐き、imabs関数を使わなくてはならず面倒ですが、Rでは区別がなく使えるので便利です。結果(FFT_spec.csv)は以下の通りで、これもExcelと同じでした。

,anafreq,spec
1,0,21.3114464479788
2,0.09765625,21.3278734470867
3,0.1953125,21.3772111473314
･･･
1022,99.70703125,21.4596305260238
1023,99.8046875,21.3772111473315
1024,99.90234375,21.327873447087

波形グラフとスペクトルの結果も、当然ながらExcelと同じものでした。
3. Python
　Pythonのプログラムは、Rのプログラムをそのまま水平展開したものです。

import numpy as np
import matplotlib.pyplot as plt

x = np.loadtxt('FourierData.csv')         #Load data

datanum = len(x)                          #data number
sampintv = 0.01                           #sampling interval[sec]
funfreq = 1/(datanum * sampintv)          #fundamantal frequency

anafreq = funfreq * np.arange(0, datanum) #analysis frequency
res = np.fft.fft(x)                       #FFT results
spec = abs(res)                           #FFT spectra

#Graph output
plt.xlim(1, 128)                          #data plot(128data)
plt.plot(x)
plt.show()

plt.plot(anafreq, spec)                   #spectra 
plt.show()

#File output
np.savetxt("FFT_res.csv",res)
np.savetxt("FFT_spec.csv",np.c_[anafreq, spec])

11行目でFFTを実行し、12行目でスペクトルを求めています。Pythonもabs関数で虚数の計算も行ってくれます。便利ですね。FFTの結果(FFT_res.csv)は以下の通りで、Excel、Rと同じ結果です。

 (2.131144644797880972e+01+0.000000000000000000e+00j)
 (2.131987775211241853e+01+-5.839506913524710185e-01j)
 (2.134521015412475009e+01+-1.169256136968316984e+00j)
･･･
 (2.138755956497289290e+01+1.757281470963324432e+00j)
 (2.134521015412481759e+01+1.169256136968219950e+00j)
 (2.131987775211233327e+01+5.839506913524976639e-01j)

スペクトルの結果(FFT_spec.csv)も同様でした。

0.000000000000000000e+00 2.131144644797880972e+01
9.765625000000000000e-02 2.132787344708677679e+01
1.953125000000000000e-01 2.137721114733136574e+01
･･･
9.970703125000000000e+01 2.145963052602384735e+01
9.980468750000000000e+01 2.137721114733142613e+01
9.990234375000000000e+01 2.132787344708669153e+01

波形グラフとスペクトルの結果も、当然ながらExcel、Rと同じものでした。
　Excel、R、Pythonで「フーリエ解析」をやってみましたが、どれを使っても手間暇は変わらないように思いました。Excelは虚数の絶対値を求めるのに、abs関数ではなく、imabs関数を使わなくてはならない所が唯一面倒に思いました。

• コメント : 0

• カテゴリ : データ解析

2018/03/18

Excel分析ツールを見直す(10)

　今回はフーリエ解析です。時系列データを周波数成分(時間間隔)ごとに分離・抽出してくれるので、使い方によっては何かと便利だなと思うのですが、必要性に迫られなかったせいか、この分析ツールも今まで使ったことがありませんでした。学生時代にフーリエ解析は勉強した経験がありますが、あまり良く覚えていません。勉強って、必要に迫られて日頃から使っていないと、こんなもんですね・・。

1. Excel
　まずは、解析用に波形データを作成しました。赤字で示したような振幅Aと角速度ωの異なる3つのSinカーブを合成しました。振幅Aはそれぞれ、4、2、10、周波数fはそれぞれ40、20、5Hzとしました。角速度は2πfで計算できます。
　時間tはsampling周期0.01sec(100Hz)で設定し(ワークシートB列)、C列に合成波形の値を計算しました。プロットすると、ギザギザしたSinカーブのような波形になりました。この波形にフーリエ解析を実行します。ここで注意することは、データ数が2ⁿ個でないと解析ができない点です。今回はデータ数1024(2¹⁰)で実行しました。

　Excelのメニューから、「データ」－「データ分析」を選択すると、「データ分析ウィンドウ」が立ち上がりますので、「フーリエ解析」を選択し、OKボタンを押します。　入力範囲にC列のデータ(合成されたギザギザSinカーブ)、出力先をD列に指定し、OKボタンを押します。　解析結果が表示されました。一見して、解析結果に虚数を含んでいますので、一瞬面喰いますが、それぞれn番目の周波数成分を表しています。細かい説明はこのサイトやこのサイトを参考にさせていただき、最終形のスペクトルまで持って行きます。解析周波数(赤枠)に対するスペクトル(青枠)をプロットすると、下図の通りです。　初めに仕込んでおいた、40、20、5Hzの周波数成分を抽出できました。ここで、解析可能周波数は今回の場合はサンプリングの定理より50Hzなので、高周波数側のピークは無視しています。

一通り、Excelでもフーリエ解析ができるようですね。感心、関心です・・。

• コメント : 0

• カテゴリ : データ解析

2018/03/11

DLLファイルを作成し、C++、C#、Pythonで呼び出す(2)

　前回の続きです。今回は前回作成したDLLファイルの関数をC++、C#、Pythonで呼び出したいと思います。

２．C++でDLLファイルを呼び出し
　Visual Studio 2017を立ち上げ、「新しいプロジェクト - Visual C++ - Windowsデスクトップ」の中の「Windowsコンソールアプリケーション(Visual C++)」を選択し、プロジェクト名を入れて、OKボタンを押します。プロジェクト名は「DllRead」としました。ソリューションエクスプローラには、2つのソースファイル(DllRead.cpp、stdafx.cpp)と2つのヘッダファイル(stdafx.h、targetver.h)が作成されていましたが、追記・修正するのは「DllRead.cpp」のみです。以下の通り、記載して「x64」でビルドしました。

// DllRead.cpp : アプリケーションのエントリ ポイントを定義します。
//

#include "stdafx.h"
#include <windows.h>
#include <iostream>

typedef double(*RFunc)(double x1, double x2);

int main(void)
{
	//Load DllTest.dll
	HINSTANCE handle_Dll = LoadLibrary(TEXT("DllTest.dll"));
	if (!handle_Dll) {
		std::cout << "DLL load error";
		return 1;
	}
	//std::cout << "DLL loaded \n";
	
	//DllAdd Function
	RFunc Func_Add = (RFunc)GetProcAddress(handle_Dll, "DllAdd");
	std::cout <<  Func_Add(3.5, 7.2) << std::endl;

	//DllSub Function
	RFunc Func_Sub = (RFunc)GetProcAddress(handle_Dll, "DllSub");
	std::cout << Func_Sub(2.2, 7.2) << std::endl;

	//DllMul Function
	RFunc Func_Mul = (RFunc)GetProcAddress(handle_Dll, "DllMul");
	std::cout << Func_Mul(5.0, 3.0) << std::endl;

	//DllDiv Function
	RFunc Func_Div = (RFunc)GetProcAddress(handle_Dll, "DllDiv");
	std::cout << Func_Div(8.0, 2.0) << std::endl;
	std::cout << Func_Div(8.0, 0.0) << std::endl;
	
	//DLLを解放
	FreeLibrary(handle_Dll);
	return 0;
}

実行結果は以下の通りです。問題なく関数が呼び出され、実行されていることを確認できました。

３．C#でDLLファイルを呼び出し
　C#では「.NET Framework」と「.NET Core」の2種類を確認しました。

３．１．.NET Framework
　まず、プロジェクトの作成です。「新しいプロジェクト - Visual C#」の中の「コンソールアプリ(.NET Framework)」を選択し、プロジェクト名を入れて、OKボタンを押します。プロジェクト名は「DllReadCS」としました。Program.csを以下の通り、追記・修正し、ビルドしました。

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Runtime.InteropServices; //DLL import

namespace DllReadCS
{
    class Program
    {
        [DllImport("DllTest.dll", EntryPoint = "DllAdd")]
        public static extern double DllAdd(double x1, double x2);

        [DllImport("DllTest.dll", EntryPoint = "DllSub")]
        public static extern double DllSub(double x1, double x2);

        [DllImport("DllTest.dll", EntryPoint = "DllMul")]
        public static extern double DllMul(double x1, double x2);

        [DllImport("DllTest.dll", EntryPoint = "DllDiv")]
        public static extern double DllDiv(double x1, double x2);

        static void Main(string[] args)
        {
            Console.WriteLine(DllAdd(3.5, 7.2));
            Console.WriteLine(DllSub(2.2, 7.2));
            Console.WriteLine(DllMul(5.0, 3.0));
            Console.WriteLine(DllDiv(8.0, 2.0));
            Console.WriteLine(DllDiv(8.0, 0.0));
        }
    }
}

この実行結果も問題ないことを確認できました。
３．２．.NET Core
　「新しいプロジェクト - Visual C#」の中の「コンソールアプリ(.NET Core)」を選択し、プロジェクト名を入れて、OKボタンを押します。プロジェクト名は「DllReadCS2」としました。Program.csを以下の通り、追記・修正し、ビルドしました。プログラムの中身は３．１．と全く同じです。

using System;
using System.Runtime.InteropServices; //DLL import

namespace DllReadCS2
{
    class Program
    {
        [DllImport("DllTest.dll", EntryPoint = "DllAdd")]
        public static extern double DllAdd(double x1, double x2);

        [DllImport("DllTest.dll", EntryPoint = "DllSub")]
        public static extern double DllSub(double x1, double x2);

        [DllImport("DllTest.dll", EntryPoint = "DllMul")]
        public static extern double DllMul(double x1, double x2);

        [DllImport("DllTest.dll", EntryPoint = "DllDiv")]
        public static extern double DllDiv(double x1, double x2);

        static void Main(string[] args)
        {
            Console.WriteLine(DllAdd(3.5, 7.2));
            Console.WriteLine(DllSub(2.2, 7.2));
            Console.WriteLine(DllMul(5.0, 3.0));
            Console.WriteLine(DllDiv(8.0, 2.0));
            Console.WriteLine(DllDiv(8.0, 0.0));
        }
    }
}

dotnetコマンドでプログラムを実行させた結果は以下の通りです。これも問題なしです。.NET Coreはマルチプラットフォームに対応しているので、念のため、Macで確認してみました。結果はNG。Windows環境で作成したDLLファイルの読み込みは当然利用できない結果でした。

４．PythonでDLLファイルを呼び出し
　PythonでDLLファイルを呼び出すときは、ctypesライブラリを利用します。ポイントは各関数について、argtypesで引数、restypeで戻り値を設定する所です。プログラムは以下の通りです。

import ctypes
from ctypes import cdll

mdll = cdll.LoadLibrary('DllTest.dll')

#DllAdd Function
Func_Add = mdll.DllAdd
Func_Add.argtypes = [ctypes.c_double, ctypes.c_double]
Func_Add.restype = ctypes.c_double

#DllSub Function
Func_Sub = mdll.DllSub
Func_Sub.argtypes = [ctypes.c_double, ctypes.c_double]
Func_Sub.restype = ctypes.c_double

#DllMul Function
Func_Mul = mdll.DllMul
Func_Mul.argtypes = [ctypes.c_double, ctypes.c_double]
Func_Mul.restype = ctypes.c_double

#DllDiv Function
Func_Div = mdll.DllDiv
Func_Div.argtypes = [ctypes.c_double, ctypes.c_double]
Func_Div.restype = ctypes.c_double

print(Func_Add(3.5, 7.2))
print(Func_Sub(2.2, 7.2))
print(Func_Mul(5.0, 3.0))
print(Func_Div(8.0, 2.0))
print(Func_Div(8.0, 0.0))

　実行結果は以下の通りです。これまた問題なしです。
　C++、C#はVisual Studioでビルドするのに手間がかかりますが、Pythonならばすぐに確認できそうで便利ですね。今後、いろいろ使えそうです。

• コメント : 0

• カテゴリ : データ解析

2018/03/11

DLLファイルを作成し、C++、C#、Pythonで呼び出す(1)

　前回、PythonからR.dllの関数を呼び出すお話をしました。ただ、DLLファイルのことをあまり理解できていなかったので、今回、一から勉強し直しました。その備忘録です。
　まず、実際に自分でDLLファイルを作成する所から始めました。DLLファイルはC++で作成し、その後に、C++、C#、Pythonの3パターンでDLL内の関数を呼び出すテストをしました。それでは、一つずつ見て行きましょう。

１．C++でDLLファイルを作成
　C++でDLLファイルを作成する方法をネットで調査していると、いろいろ情報・やり方があって、どれを選択すればよいのか分からないですね。今回は「__declspec(dllexport)」をコードに記載した「オーソドックスな(昔ながらの)方法?」と「defファイルを使う方法」の2パターンを試してみました。

１．１．オーソドックスな方法(__declspec(dllexport)を利用)
　Visual Studio 2017を立ち上げ、「新しいプロジェクト - Visual C++ - Windowsデスクトップ」の中の「ダイナミックリンクライブラリ(DLL)」を選択し、プロジェクト名を入れて、OKボタンを押します。プロジェクト名は「DllTest」としました。ソリューションエクスプローラには、3つのソースファイル(dllmain.cpp、DllTest.cpp、stdafx.cpp)と2つのヘッダファイル(stdafx.h、targetver.h)が作成されていましたが、追記・修正するのは「DllTest.cpp」のみです。以下の通り、記載して「x64」でビルドしました。初め「x86」でビルドして、PythonでDLLを読み込もうとしていたのですが、エラーが出てしまいました。自分のPC環境に合わせてビルドする必要があります。

// DllTest.cpp : DLL アプリケーション用にエクスポートされる関数を定義します。
//

#include "stdafx.h"
#include <iostream>

extern "C" __declspec(dllexport) double DllAdd(double a, double b)
{
	return a + b;
}

extern "C" __declspec(dllexport) double DllSub(double a, double b)
{
	return a - b;
}

extern "C" __declspec(dllexport) double DllMul(double a, double b)
{
	return a * b;
}

extern "C" __declspec(dllexport) double DllDiv(double a, double b)
{
    if (b == 0)
    {
        printf("Divided by zero!");
        return 0;
    }
    return a / b;
}

作成した関数は、定番の?加減乗除の関数(DllAdd、DllSub、DllMul、DllDiv)です。main関数は記載せず、個別の関数について、「extern "C" __declspec(dllexport)」を付けて作成しました。作成したDLLファイルをdumpbinで確認すると以下の通りで、エクスポートされた関数が表示されました。問題なくDLLファイルが作成されています。

１．２．defファイルを使う方法
　次に、defファイルを使う方法ですが、その前に、「__declspec(dllexport)」を記載せずに以下のコードをビルドすると、ビルドは成功してDLLファイルは作成されます。

// DllTest.cpp : DLL アプリケーション用にエクスポートされる関数を定義します。
//

#include "stdafx.h"
#include <iostream>

double DllAdd(double a, double b)
{
	return a + b;
}

double DllSub(double a, double b)
{
	return a - b;
}

double DllMul(double a, double b)
{
	return a * b;
}

double DllDiv(double a, double b)
{
    if (b == 0)
    {
        printf("Divided by zero!");
        return 0;
    }
    return a / b;
}

しかし、dumpbinで確認すると以下の通りで、エクスポートされた関数が表示されません。当然、外部から関数を呼び出すことができませんでした。このコードで、関数をエクスポートするためには、エクスポート定義ファイル(defファイル)を作成する必要があります。ソリューションエクスプローラのプロジェクト(DllTest)を選択、右クリックし、「追加 - 新しい項目」を選択します。「Visual C++ - コード - モジュール定義ファイル(.def)」を選択し、ファイル名(DllTest.def)を入力します。　ソリューションエクスプローラのソースファイルにDllTest.defが追加されましたので、以下の内容を記入します。

LIBRARY DllTest

EXPORTS
	 DllAdd @1
	 DllSub @2
	 DllMul @3
	 DllDiv @4

ビルドする前に、プロジェクトのプロパティページの「リンカー - 入力」のモジュール定義ファイルの欄にこのファイル名(DllTest.def)を記入し、「x64」でビルドしました。作成したDLLファイルをdumpbinで確認すると以下の通りで、エクスポートされた関数が表示されました。「__declspec(dllexport)」でビルドしたDLLファイルの関数表記より関数の引数まで記載されており、大変分かりやすいですね。
　2パターンでDLLファイルを作成しましたが、2つとも問題なく機能しました。次回は、この作成したDLLファイルの関数をC++、C#、Pythonで呼び出したいと思います。

• コメント : 0

• カテゴリ : データ解析

2018/03/10

PythonからR.dllの関数を呼び出す

　Rをインストールしていない環境でも、R.dllの中の関数を呼び出すことができれば、何かと便利ではないかと考えていました。参考サイトにVC++2010から、R.dllの関数を呼び出す「見本」がありましたので、今回はPythonから呼び出してみようと思います。

　まず、R.dllの中身を調べてみます。使っているRはVer.3.4.2で、R.dllは30MB程のファイルです。dumpbinコマンドを実行するために、Visual Studio 2017をインストールした際にセットアップされている「Developer Command Prompt for VS 2017」を実行します。プロンプトが出てきたら、以下のコマンドを実行します。後で内容を確認するために、リダイレクトでファイル出力(d.txt)しました。

dumpbin /EXPORTS R.dll >d.txt

　実行結果の中身は以下の通りで、R.dllの中の関数を確認することができます。

Dump of file R.dll

File Type: DLL

  Section contains the following exports for R.dll

    00000000 characteristics
    59CCC296 time date stamp Thu Sep 28 18:36:22 2017
        0.00 version
           1 ordinal base
        1118 number of functions
        1118 number of names

    ordinal hint RVA      name

          1    0 00114130 ATTRIB
          2    1 01D4D848 AllDevicesKilled
          3    2 00115340 BODY
　　　　　･･･････
        554  229 00201C20 Rf_beta
　　　　　･･･････
        642  281 00201100 Rf_gammafn
           ･･･････
       1116  45B 0022EF00 xdr_double
       1117  45C 0022F280 xdrmem_create
       1118  45D 001472D0 xerbla_

  Summary

        1000 .CRT
      133000 .bss
       1F000 .data
        7000 .edata
        5000 .idata
       12000 .pdata
       F7000 .rdata
        A000 .reloc
     18B6000 .rodata
        1000 .rsrc
      356000 .text
        1000 .tls
       16000 .xdata

　今回は、見本のガンマ関数の他にベータ関数も呼び出しました。R.dllの中の関数はそれぞれRf_gammafn、Rf_betaです。

　事前にRで実行した結果は、以下の通りです。

> gamma(0.5)
[1] 1.772454
> beta(0.2,0.5)
[1] 6.268653

　Pythonで実行したプログラムは以下の通りです。

import ctypes
from ctypes import cdll

libc = cdll.LoadLibrary('R.dll')

#Rf_gammafn Function
gamma = libc.Rf_gammafn
gamma.argtypes = [ctypes.c_double]
gamma.restype = ctypes.c_double

#Rf_beta Function
beta = libc.Rf_beta
beta.argtypes = [ctypes.c_double, ctypes.c_double]
beta.restype = ctypes.c_double

print(gamma(0.5))
print(beta(0.2, 0.5))

実行結果は、以下の通りで同じ結果が得られました。

1.772453850905516
6.268653124086035

　PythonでR.dllの関数を簡便に呼び出せると、わざわざC++でコンパイルする面倒もなくなりますね。また、R.dllの関数名と実際のRの関数名は似ていますが、異なるので、試行錯誤で関数を探す手間も含めてすごく楽になると思いました。
　今回、簡単にできたように見えますが、ここまで来るのにかなり難儀しました･･･。DLLファイルのことについてあまり良く理解していなかったので、その外堀から埋めていきました。プログラムの解説を含めて、次回にそのお話をしたいと思います。

• コメント : 0

• カテゴリ : データ解析

2018/03/08

Excel分析ツールを見直す(9)

　前回の続きです。今回は「F検定」をRとPythonで求めます。

2. R
　Rのプログラムは以下の通りです。

# F-test
va<-c(301, 311, 325, 291, 388, 412, 325, 361, 287)
vb<-c(197, 180, 247, 260, 247, 199, 179, 134, 163, 200)
var.test(va,vb)

#95 percent confidence interval
vr<-var(va)/var(vb)
vr/qf(0.025,8,9)
vr/qf(0.975,8,9)

#F-Distribution
x<-seq(0,10,0.1)
plot(x, df(x, 8, 9),ylim=c(0,1),type="l")
curve(pf(x, 8, 9),add=TRUE,type="l",col="red")
lower <- qf(0.05, 8, 9)
upper <- qf(0.95, 8, 9)
sprintf("Lower 0.05  : %f",lower)
sprintf("Upper 0.05  : %f",upper)
abline(v = lower,col="blue")
abline(v = upper,col="blue")

　F検定のみならば、2~4行目の3行のみで完了し、結果は以下の通りです。

	F test to compare two variances

data:  va and vb
F = 1.2001, num df = 8, denom df = 9, p-value = 0.7859
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.2925646 5.2290589
sample estimates:
ratio of variances 
          1.200087

　結果の4行目のFはF値、dfは自由度です。p-valueは前回のExcelで算出された片側のp値の2倍の値と同じです。5行目は「対立仮説：真の分散比は1に等しくない」との親切な覚書です。つまり、帰無仮説は「分散比は１に等しい(等分散である)」ということです。前回と同様に、p値から判断して帰無仮説は棄却できない、つまり「等分散である」ことを棄却できない(等分散である)ということになります。

　その次の「95 percent confidence interval:」でふと止まってしまいました。「95％信頼区間」って、いったい何だろうと思いながら、ネットで調べていると、プログラムの7~9行目に書いてあるように、分散比をF分布の上下限2.5%の確率密度(Rコマンドではqf)で割った値であることが分かりました。12行目以降は、確率密度関数と累積分布関数をRでプロットしたものです。15、16行目で上下限5%の確率密度を求めており、結果は以下の通りで、グラフ上に青線で示しています。

Lower 0.05  : 0.295148
Upper 0.05  : 3.229583

　Rでも、Excelと同様に簡単にF検定ができることが分かりました。

3. Python
　Pythonのプログラムは以下の通りです。

import numpy as np
from scipy import stats 
from matplotlib import pyplot as plt

#data
a = np.array([301, 311, 325, 291, 388, 412, 325, 361, 287])
b = np.array([197, 180, 247, 260, 247, 199, 179, 134, 163, 200])

#variance
var_a = np.var(a, ddof = 1)
var_b = np.var(b, ddof = 1)
print("Var a    : " + str(var_a))
print("Var b    : " + str(var_b))

#F-value
if var_a > var_b :
    Fvalue = var_a / var_b
else:
    Fvalue = var_b / var_a
print("F-value  : " + str(Fvalue))

#Degrees of Freedom
print("Freedom a: " + str(len(a)-1))
print("Freedom b: " + str(len(b)-1))

#F-border
upper = stats.f.ppf(0.95, len(a)-1, len(b)-1)
print("F-border  : " + str(upper))
lower = stats.f.ppf(0.05, len(a)-1, len(b)-1)
print("F-border  : " + str(lower))

#F-distribution
x = np.arange(0,10,0.1)
plt.xlim(0, 10)
plt.ylim(0, 1)
plt.plot(x, stats.f.pdf(x, len(a)-1, len(b)-1))
plt.plot(x, stats.f.cdf(x, len(a)-1, len(b)-1))
plt.vlines(lower,0,1,colors="blue",linestyles="dashed")
plt.vlines(upper,0,1,colors="blue",linestyles="dashed")
plt.show()

　ExcelやRのようにワンタッチで答えが出るような関数を見つけることができなかったので、一つずつ計算を進めていきます。10、11行目で分散、16~19行目でF値(分散比)、23、24行目で自由度、27、29行目で上下限5%の確率密度を求めました。結果は以下の通りで、当然ですが、ExcelやRと同じ結果でした。

Var a    : 1945.5277777777776
Var b    : 1621.1555555555558
F-value  : 1.2000870435354751
Freedom a: 8
Freedom b: 9
F-border  : 3.2295826126867744
F-border  : 0.29514804057607874

　33行目以降は、確率密度関数と累積分布関数をPythonでプロットしたものです。結果は見ての通りでExcel、Rと同じでした。

　Excel、R、Pythonと見てきて、ExcelとRは簡単に「F検定」できることが良く分かりました。今回も大変勉強になりました。

• コメント : 0

• カテゴリ : データ解析

2018/03/07

Excel分析ツールを見直す(8)

　前回からしばらく間が空きました。今回は「F検定」です。この機能も今まで使ったことがありませんでした。F検定やF分布は「等分散の検定」や「分散分析」で利用されるとのことです。早速使ってみましょう。

1.Excel
　Excelのメニューから、「データ」－「データ分析」を選択すると、「データ分析ウィンドウ」が立ち上がりますので、「F検定：2標本を使った分散の検定」を選択し、OKボタンを押します。データはA列の9個のデータとB列の10個のデータです。このA列とB列のデータの「分散が等しいか等しくないか」を検定します。変数1にA列のデータ、変数2にB列のデータを指定し、出力先を指定し、OKボタンを押します。　　
　結果が表示されました。・・が、これらが何を意味しているか、直感的に理解できないところが、統計ビギナーの悲しい所です。平均、分散、観測数、自由度は直感的にこのような計算をしているんだなと分かりますが・・・。これらの計算値はそれぞれExcel関数を持っているので再確認しました。「観測された分散比」はその名の通り、求められた2つの分散値の比を表しており、F値と呼ばれています。等分散ならばF値が1に近いということですね。
「P(F<=f)片側」はp値のことですね。
このF検定は、
　帰無仮説H0　：　データ間の分散は等しい。
　対立仮説H1　：　データ間の分散は等しくない。
に基づき、行われるので、p値(0.392929：両側ならば2をかけて0.785858)が0.05より大きいので帰無仮説H0は棄却されないことより、「分散は等しくない」とは言えない。つまり、「分散は等しい」ということになります。あぁ、しんど。帰無仮説と対立仮説のことを考えるといつも混乱しますね・・(私だけでしょうか??)。

「F境界値片側」はまた、何のことかピンときません。F分布とF検定を十分理解する必要がありそうです。このサイトを参考にF分布の式は次の形になることが分かりました。Γ(ガンマ)関数を使ったすごい形です。n1、n2は各データの自由度の値です。　(1)式を用いて、F分布の確率密度関数を求めてプロットしました。こんな複雑な計算をわざわざしなくても、Excelには便利な「F.DIST関数」がありました。
　　F.DIST(X,自由度1,自由度2,TRUE/FALSE)
で、4つ目の引数をTRUEにすると「累積分布関数の値」、FALSEにすると「確率密度関数の値」を返すとのこと。確かに、FALSEにした時の値は(1)式で計算した値と同じでした。

話を戻して、「F境界値片側」ですが、累積分布関数の値(グラフのオレンジ色の線)のyが0.95の時のxの値だったんですね。ようやく納得です。もう片側はxが0.3付近にあります。つまり、F値(分散比)は1.2で、0.3~3.23の間に入っていることからも分散は等しくないとは言えない(つまり、等しい)ということですね。

長くなりましたので、次回はR、Pythonで同じことをやってみます。

• コメント : 0

• カテゴリ : データ解析