pythonでフィッティングするときのためのめも。
コード
fitlin.py
""" 線形フィッティング """ import matplotlib.pyplot as plt import numpy as np from scipy.optimize import curve_fit # fitting関数 def fit(x, a, b): return a * x + b # データ入力 array_x = np.array([0, 1, 2, 3, 4, 5, 6]) array_y = np.array([0.2, 1.4, 2.2, 3.1, 4.5, 5.2, 5.8]) # paramはフィッティングパラメータ、covは共分散行列 param, cov = curve_fit(fit, array_x, array_y) # 共分散のルートでフィッティングパラメータの誤差 error = np.sqrt(np.diag(cov)) # fitting係数出力 print("a:{:.6f} b:{:.6f}".format(param[0], param[1])) print("sigma_a:{:.6f} sigma_b:{:.6f}".format(error[0], error[1])) # fitting関数からデータ作成 array_y_fit = fit(array_x, param[0], param[1]) # グラフ描画 plt.plot(array_x, array_y) plt.plot(array_x, array_y_fit) plt.show()
・入力データ
lognormal_dat.csv
0.01,0.0001 0.1,0.017 1,0.24 2,0.2 3,0.13 4,0.09 5,0.06 6,0.04 7,0.03 8,0.03 9,0.02 10,0.01
fitlognorm.py
""" 対数正規分布をフィッティング """ import os import sys import matplotlib.pyplot as plt import numpy as np import pandas as pd from scipy.optimize import curve_fit # 対数正規分布 def lognormal(x, a, sigma, mu, c): return a * (1 / (np.sqrt(2 * np.pi * sigma**2) * x) * np.exp(-((np.log(x) - mu) ** 2) / (2 * sigma**2))) + c if __name__ == "__main__": argvs = sys.argv argc = len(argvs) # エラー処理 if argc != 2: print("引数1：対数正規分布データ") input("キーを何か押すと終了します。") sys.exit() # ファイル保存設定 savedir = os.path.dirname(os.path.abspath(argvs[1])) filename, ext = os.path.splitext(os.path.basename(argvs[1])) # データ読み込み data = np.loadtxt(argvs[1], delimiter=",") pdf_x = data[:, 0] pdf_y = data[:, 1] # === 対数正規分布フィッティング === # paramはフィッティングパラメータ、covは共分散行列 param_ini = (1, 1, 1, 0) # 初期値 param, cov = curve_fit(lognormal, pdf_x, pdf_y, p0=param_ini) # 共分散のルートでフィッティングパラメータの不確かさ err = np.sqrt(np.diag(cov)) # fitting係数出力 print("・対数正規分布") print("a : {:.6f} +- {:.6f}".format(param[0], err[0])) print("sigma : {:.6f} +- {:.6f}".format(param[1], err[1])) print("mu : {:.6f} +- {:.6f}".format(param[2], err[2])) print("c : {:.6f} +- {:.6f}".format(param[3], err[3])) print("\n") # fitting関数からデータ作成 pdf_y_fit = lognormal(pdf_x, param[0], param[1], param[2], param[3]) # フィッティングデータ保存 # パラメータ name_pdf = np.array(["a", "sigma", "mu", "c", "a_sigma", "sigma_sigma", "mu_sigma", "c_sigma"]) data_pdf = np.array([param[0], param[1], param[2], param[3], err[0], err[1], err[2], err[3]]) data_pdf_df = pd.DataFrame([np.round(data_pdf, 6)]) data_pdf_df.to_csv(savedir + "\\" + filename + "_pdf_fitting_param.csv", header=name_pdf, index=False) # フィッティングy name_pdf = np.array(["x", "fit_y"]) data_pdf = np.vstack((pdf_x, pdf_y_fit)).T data_pdf_df = pd.DataFrame(np.round(data_pdf, 6)) data_pdf_df.to_csv(savedir + "\\" + filename + "_pdf_fitting.csv", header=name_pdf, index=False) # === グラフ描画 === plt.plot(pdf_x, pdf_y, "o", label="data") plt.plot(pdf_x, pdf_y_fit, label="fit") plt.title("log-normal PDF\n{0}".format(filename)) plt.xlabel("x") plt.ylabel("PDF") plt.grid(which="major", axis="x", color="gray", alpha=0.8, linestyle="--", linewidth=0.5) plt.grid(which="major", axis="y", color="gray", alpha=0.8, linestyle="--", linewidth=0.5) plt.legend() plt.savefig(savedir + "\\" + filename + "_pdf.png") plt.show()
fitlognormlhist.py
""" ヒストグラムをフィッティング 累積分布関数を生成してフィッティング """ import os import sys import matplotlib.pyplot as plt import numpy as np import pandas as pd from scipy.optimize import curve_fit from scipy.special import erf from scipy.stats import norm # 対数正規分布 def lognormal(x, a, sigma, mu, c): return a * (1 / (np.sqrt(2 * np.pi * sigma**2) * x) * np.exp(-((np.log(x) - mu) ** 2) / (2 * sigma**2))) + c # 対数正規分布 累積分布関数 def lognormal_erf(x, sigma, mu): return (1 / 2) * (1 + erf((np.log(x) - mu) / (np.sqrt(2 * sigma**2)))) if __name__ == "__main__": # === 対数正規分布作成 === np.random.seed(1) # 乱数を再現できるように mu = 3 sigma = 1 data_num = 1000 lognorm_dat = np.random.lognormal(mu, sigma, data_num) # dataframeに格納 label = "lognormal" df = pd.DataFrame(lognorm_dat, columns=[label]) # ヒストグラム作成 # ヒストグラムをplt.histで作成するのでここでfigureオブジェクト生成 plt.figure(1) plt.style.use("ggplot") freq, bins, _ = plt.hist( [df[label]], bins=100, alpha=0.5, # グラフの透過率 0が透明、1が不透明 未指定の場合は1 # color=['#f46d43'] #グラフの色 未指定の場合は自動 ec="navy", label=[label], # グラフのラベル # range=(x_range_min, x_range_max), # グラフの横軸の範囲 未指定の場合はデータの最大、最小値に合わせられる rwidth=1.0 # グラフの棒の幅 身指定の場合は1 # stacked=False #複数データを積み上げグラフにするか、横に並べるか 未指定の場合はFalse ) # 分布のx,ｙの値を設定 # plt.histは端の値が返ってくるのでbins+1のデータ数がある # フィッティングはビンの中心の値を使うのでビンの差分を2で割って足す pdf_x = bins[0:-1] + np.diff(bins) / 2 pdf_y = freq # 対数正規分布フィッティング # paramはフィッティングパラメータ、covは共分散行列 param_ini = (10000, sigma, mu, 0) # 初期値 param, cov = curve_fit(lognormal, pdf_x, pdf_y, p0=param_ini) # 共分散のルートでフィッティングパラメータの不確かさ err = np.sqrt(np.diag(cov)) # fitting関数からデータ作成 pdf_y_fit = lognormal(pdf_x, param[0], param[1], param[2], param[3]) # fitting係数出力 print("・対数正規分布") print("a : {:.6f} +- {:.6f}".format(param[0], err[0])) print("sigma : {:.6f} +- {:.6f}".format(param[1], err[1])) print("mu : {:.6f} +- {:.6f}".format(param[2], err[2])) print("c : {:.6f} +- {:.6f}".format(param[3], err[3])) print("\n") # フィッティングデータ保存 # パラメータ name_pdf = np.array(["a", "sigma", "mu", "c", "a_sigma", "sigma_sigma", "mu_sigma", "c_sigma"]) data_pdf = np.array([param[0], param[1], param[2], param[3], err[0], err[1], err[2], err[3]]) data_pdf_df = pd.DataFrame([np.round(data_pdf, 6)]) data_pdf_df.to_csv("pdf_fitting_param.csv", header=name_pdf, index=False) # フィッティングy name_pdf = np.array(["pdf_x", "pdf_y", "pdf_y_fit"]) data_pdf = np.vstack((pdf_x, pdf_y, pdf_y_fit)).T data_pdf_df = pd.DataFrame(np.round(data_pdf, 6)) data_pdf_df.to_csv("pdf_fitting.csv", header=name_pdf, index=False) # === 累積分布フィッティング === cdf_x = sorted(df[label]) cdf_y = [i / data_num for i in range(data_num)] param, cov = curve_fit(lognormal_erf, cdf_x, cdf_y) err = np.sqrt(np.diag(cov)) print("・対数正規分布 累積分布関数") print("sigma : {:.6f} +- {:.6f}".format(param[0], err[0])) print("mu : {:.6f} +- {:.6f}".format(param[1], err[1])) print("\n") cdf_y_fit = lognormal_erf(cdf_x, param[0], param[1]) # フィッティングデータ保存 # パラメータ name_cdf = np.array(["sigma", "mu", "sigma_sigma", "mu_sigma"]) data_cdf = np.array([param[0], param[1], err[0], err[1]]) data_cdf_df = pd.DataFrame([np.round(data_cdf, 6)]) data_cdf_df.to_csv("cdf_fitting_param.csv", header=name_cdf, index=False) # フィッティングy name_cdf = np.array(["cdf_x", "cdf_y", "cdf_y_fit"]) data_cdf = np.vstack((cdf_x, cdf_y, cdf_y_fit)).T data_cdf_df = pd.DataFrame(np.round(data_cdf, 6)) data_cdf_df.to_csv("cdf_fitting.csv", header=name_cdf, index=False) # === グラフ描画 === plt.plot(pdf_x, pdf_y_fit, label="fit") plt.title("log-normal") plt.xlabel("x") plt.ylabel("Frequency") plt.legend() # 凡例表示 plt.savefig("lognormal_hist.png") plt.figure(2) plt.plot(cdf_x, cdf_y, "o", label="data") plt.plot(cdf_x, cdf_y_fit, label="fit") plt.title("log-normal CDF") plt.xlabel("x") plt.ylabel("CDF") plt.legend() plt.savefig("lognormal_hist_cdf.png") plt.show()
