日本包装学会誌 Vol.28 No.5 (2019)

一般論文～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～

両面段ボールに対する均質化解析とベイズ推定を用いた確率的妥当性評価技術の開発
島 村

雅 彦*，江 澤 良 孝**，田 村 善 昭***，高 清 水 聖*，佐 藤 大 亮****Development of Homogenization Analysis and Probabilistic
Validation Assessment Technology Using Bayesian Inference
for Single Wall Corrugated Fiberboards
Masahiko SHIMAMURA*, Yoshitaka EZAWA**, Yoshiaki TAMURA***,
Satoru TAKASHIMIZU* and Daisuke SATOU****
ボールに数値解析を適用する場合、その構造が周期性を有するため均質化解析は有効な手法である。しかし、均質化解析は 3 方向に周期性をもつ構造に対して理論が成り立つため、両面段ボールに均質化解析を適用した場
合、解析結果の変位が負荷試験の変位測定値と合わないことが知られている。本論文では測定値を利用して、両面段ボールに均質化解析が適用できる弾性係数を算出する方法を提案する。一方、提案による弾性係数は測定値のばらつきを含む。解析結果を有効に活用するには、このばらつきを評価する必要があるため、弾性係数の妥当性検証としてベイズ推定を利用した確率論的評価を行う。ここでは、推定精度向上を目的に情報量基準を用いてハイパーパラメータを最適化する方法を提案する。本論文では以上に述べた 2 つの提案から、両面段ボールに適用可能な弾性係数を算出し、得られた弾性係数のばらつきをベイズ推定により精度よく推定できることを示す。

Homogenization analysis is effective for the analysis of corrugated fiberboard which has periodic feature. As
homogenization analysis theoretically assumes periodicity in three directions, results of homogenization analysis
for single wall corrugated fiberboards are reported to have unignorable discrepancy with measurements. In thispaper, we propose a method to obtain elastic moduli for single wall corrugated fiberboards in homogenization
analysis, utilizing measurement data. Obtained elastic moduli include variations. In order to estimate thesevariations, probabilistic assessment method using Bayesian inference is applied to validate the elastic moduli. In
this paper, we propose a method to optimize hyper parameter of Bayesian inference using information criterion
for highly accurate estimation. In this paper, we show the two proposal leads to the elastic moduli applicable forsingle wall corrugated fiberboards and to the highly accurate estimation of their variations using Bayesianinference.
キーワード : 両面段ボール、均質化解析、妥当性検証、ベイズ推定

Keywords ：Single wall corrugated fiberboard, Homogenization analysis, Validation, Bayesian inference*

連絡者(Corresponding author), 東洋大学大学院工学研究科(〒350-8585 埼玉県川越市鯨井 2100)
2100 Kujirai, Kawagoe City, Saitama, 350-8585, Japan, Email: s46a01310013@toyo.jp**
東洋大学計算力学研究センター客員研究員***
東洋大学総合情報学部****
カネパッケージ(株)－ 321 －

両面段ボールに対する均質化解析とベイズ推定を用いた確率的妥当性評価技術の開発ている。均質化解析は 3 次元問題に適用する

1. 緒言
段ボールは優れた力学的構造の特徴をもち、にあたり、3 方向に周期性を有する構造体を段ボール箱や緩衝材として活用されている。

対象として理論が成立している。そのため、段ボールの既存研究は数多く存在し、例えば厚さ方向に周期性のない両面段ボールに均質1)の強度に関する基礎的研究がある。

化解析を適用すると、無視し得ない誤差が生

近年、設計段階において数値解析を活用する

じることが報告されている 6 ) 。しかし、解析ことの意義は大きくなっており、効率的に製の大幅な効率化を考慮すれば実用上有用な解品開発を行うことの一助となっている。段ボ析手法であるといえる。平板を対象とした均ールを対象とした数値解析の既存研究として、質化理論の適用は、連続体をモデル化したミ

中川ら 2 ) , 3 ) は、段ボールに対して有限要素解

ク ロ ス ケ ー ル の 問 題 に 対 し て Bernoulli-析が適用可能であるか動的解析を行なうこと

Euler の仮定などに基づいた平板の運動場をで検討を行っている。ここでは、段ボールシ導入する方法 6 ) や漸近展開に基づく定式化 7 )ートを一様の単一素材と捉えたシェル要素でなどが試みられている。

松島ら

簡易モデル化し、物性データは段ボールシー

本研究では、従来の均質化解析に対して数トの状態から測定した値を用いている。さら学的理論に基づいた解の漸近展開から始まるに、簡易モデルによる段ボールの緩衝材基本式展開や厚板の面内周期性に関する数学的処

モデルおよび木箱を用いた包装品モデルにつ理を行うなどすることなく、両面段ボールにいても検討を行っている。また、三上 4 ) らは、均質化解析を適用する方法を提案する。既存

垂直圧縮試験の有限要素解析を行なうことでの研究 2 ) , 3 ) , 4 ) では、単一のシェル要素による最適なフルート形状を検討している。この解モデル化や段ボールの形状対称性を利用した

析では、ライナと中芯の物性に対して等方性

詳細モデルの簡略化により節点数、要素数を弾塑性を仮定し、原紙の抄紙方向のリングク低減するに留まっており、均質化解析の活用

ラッシュ試験から得られた応力－ひずみ曲線

はなされていない。本研究では均質化解析をに、2 直線近似より得られた弾性率、降伏応

用いることで、効率的かつ高精度な解析を行

力、接線係数を適用した非線形解析を実施しうことを目的とする。

ている。
段ボールは、その形状が複雑であることか

2. 提案手法の概要

ら一般的な有限要素解析を用いた場合、節点本研究では両面段ボールを対象に以下に述数および要素数が多くなることで、メモリー

べる 2 つの手法を提案する。ここでは、そのの負荷やそれに伴い解析にかかる計算時間が特徴について述べる。

極めて高くなる。そのため、段ボールの周期性に着目した均質化解析 5 ) の適用が試みられ－ 322 －

日本包装学会誌 Vol.28 No.5 (2019)

2.1 両面段ボールの数値解析

そこで、本研究では事前分布の最適な確率モ

2.1.1 均質化解析適用の問題

デルを選択する方法論として情報量基準

均質化解析の数学的理論は、空間的に周期WAIC(Widely

Applicable

Information

変動する係数をもつ微分方程式において、そCriterion) 1 1 ) , 1 2 ) を 用 い る こ と で 、 事 前 分 布 のの解の漸近挙動を記述するものである。その上層部に属するハイパーパラメータを最適化

ため、面内にのみ周期的な非均質性を有するする方法を提案する。このハイパーパラメー

両面段ボールに直接適用すると解析結果の変

タは、事前分布に対する事前分布を意味する

位に無視し得ない誤差が生じる。近年、平板ものであり、ハイパーパラメータから事前分に対する均質化解析の研究が進められている布を推定することが可能である。このようにが、現状では解決に至っていない 8 ) 。

情報 量基 準 WAIC に よって 推定された 事前分布を用いることで調整済み弾性係数のばらつきを精度良く推定できる。

2.1.2 弾性係数の適用方法
前節で述べた問題を解決するため、本研究では曲げ剛性を定義する弾性係数に着目し、

3. 試験測定された変位に均質化解析の変位が合うよ

調整済み弾性係数は、試験と均質化解析か

う弾性係数を調整する方法を提案する。尚、

ら算出される値である。ここでは、そのため提案により算出された弾性係数を本論文ではの曲げ試験について述べる。

調整済み弾性係数と定義する。
3.1 曲げ試験
本 研 究 で は 表 裏 ラ イ ナ LB180 、 中 芯

2.2 ばらつきの確率的評価法
測定値は様々な要因でばらついており、そMC120 から成る B-flute を供紙材として曲げの測定値に合わせるように均質化解析を行な

試験 1 3 ) を行った。LB180 と MC120 は JIS P

うため、算出された調整済み弾性係数も同様

3902 、 JIS P 3904 に 準 拠 し た も の で あ る 。

にばらつきを含む。このようなばらつきを評

また、試験片は、温度 11℃、湿度 38%の条価する方法論としてベイズ推定9),10)がある。

件で前処理を施した。 Fig.1 に は曲げ試験 装

ベイズ推定は、設計者の主観を反映した事置を示すが、測定は試験片の両端を治具で完前確率密度関数(以後、事前分布)を有効に活全に固定し、圧子により圧縮速度 10mm/min

用することで、少ない観測データからでも比で鉛直方向の荷重を与え、段ボールシートか較的良い推定を行なうことが可能である。すら300

なわち、推定の対象となるパラメータに対し段が完全に潰れるまでの荷重－変位特性を測て確信度の高い事前分布を定義することが可定した。

能であれば、精度の高い推定が可能となる。

－ 323 －

300mmの 寸法に切 り抜いた試 験片の両面段ボールに対する均質化解析とベイズ推定を用いた確率的妥当性評価技術の開発

ぞれの平均として 2 方向の回帰直線を求め、Intender

この線上に存在する 20N の荷重を基準荷重

Jig

𝑁 とし、Fig.2 中に示した黒点、すなわち基Test piece

準荷重𝑁 おける変位を指標変位𝑇

、𝑇 と

する。

Sensor本研究では、Fig.2 に拡大図で示した線形範Fig.1 Bending test

囲の一部である 5N から 20N を調整済み弾性係数の妥当性を評価するための対象とした。

また、段ボールは中芯の方向で異なる特性をもつことから、抄紙方向と垂直方向に対して各方向に 10 回、総測定数 20 回の試験を𝑇実施した。尚、各試験につき新しい試験片を用いて、圧子の初期位置や試験片の固定具合など、測定の準備段階でばらつきの影響が極力小さくなるようにした。

3.2 試験結果
曲げ試験で測定された荷重－変位特性を(a) Machine direction

Fig.2 に示す。図中の各線は、それぞれの 測定値を示している。 Fig.2 に示した拡大図 よ𝑇り、測定値には、試験機の性能(センサの分
解能、サンプリングレートなど)といった試験機に依存するばらつきがあり、また、荷重負荷の小さい段階から測定値のばらつきが見受けられることから、試験準備や試験片の作成といった人的な要因で生じるばらつきも考えられる。
(b) Cross direction

また、これらの測定された変位は、曲げ剛

Fig.2 Load ‐displacement characters

性によって決定されるため、弾性係数を調整して均質化解析の変位を測定値に合わせる。
この調整済み弾性係数を算出するにあたり指4. 数値解析

標となる変位を以下のように定める。抄紙方

ここでは、均質化理論に基づく両面段ボー

向および垂直方向の 10 回の測定値からそれルのモデル化および境界条件について述べる。

－ 324 －

日本包装学会誌 Vol.28 No.5 (2019)

4.1 均質化法の理論 5 ) ， 1 4 )4.2.1 ミクロモデル(ユニットセル)
段ボールの 1 つの段をユニットセルとする均質化解析では周期性を伴うミクロ構造からなる構造体において、ミクロ構造の 1 周期

ミクロモデルを Fig.3 に示すように作成しにあたる構成単位をユニットセル、構造全体た。このミクロモデルから式(2)による均質をマクロ構造と定義する。ユニットセルにお化弾性係数を得ることができる。Fig.3 より

け る 変 位 𝒖は 、 ミ ク ロ ひ ず み 𝑬と 擾 乱 変 位 𝒖∗

𝑥軸方向は抄紙方向、𝑧軸方向は垂直方向を示

によって下記のように表される。

している。中芯は製造過程により依存し各段で異なった形状をしているため、Harrysson

𝒖

𝑬∙𝒙

∗

𝒖

(1)

and Ristinmaa 1 5 ) に倣い正弦波としてモデル化を行った。また、ライナと中芯は異なる特

ここで、𝒖∗ はユニットセルと同じ周期性を有性をもつため、それぞれ独立した弾性係数をする擾乱変位である。また、ユニットセル境

Material(1)および Material(2)として定義す

界 上 の 相 対 す る 2 点 を𝑥 、 𝑥 と す る と 弾 性る。

係 数 を 評 価 す る た め に 、 擾 乱 変 位 𝒖∗ 𝑥
𝒖∗ 𝑥、
Material (1): Liner

は次式を満足しなければならない。
𝒖∗ 𝑥

𝒖∗ 𝑥

Material (2):Corrugated medium

2

式 (2) を 解 く こ と で 、 ミ ク ロ 構 造 の 均 質 化 弾
性 係 数 を 得 る こ と が で き る 。 ま た 、 式 (1) 、3mm(2) よ り 、 周 期 境 界 条 件 は 次 式 の よ う に 表 さ

6mm

れる。

𝑦

6mm
𝒖 𝑥

𝒖 𝑥

𝑬∙ 𝑥

𝑥

3

𝑧

式 (3) の 周 期 境 界 条 件 を 用 い る こ と で 、 ミ ク

Fig.3 Composition and dimensions of
micro model(Unit cell)𝑥

ロ構造を考慮しながらマクロ構造を解くことができる。

ユニットセルの弾性係数に関して、複両面
段 ボ ー ル で は 算 出 式 16) な ど か ら 求 め た 値 を用いることで、測定値と概ね一致した結果を得 ら れ て い る 17) 。 し か し 、 両 面 段 ボ ー ル に

4.2 有限要素モデル
両面段ボールの均質化モデルを下記のように作成した。

算出式から求めた弾性係数を用いると解析結果と測定値とで合わないことから、算出式か－ 325 －

両面段ボールに対する均質化解析とベイズ推定を用いた確率的妥当性評価技術の開発

ら求まる弾性係数はここでは用いず、第 3 章

こで与える荷重値は基準荷重𝑁 に固定して解で示した指標変位𝑇析を行なう。

、𝑇 に弾性係数を合わせるための均質化解析を行なって調整済み弾

𝐿性係数を求める方法をとる。

𝐶𝐿

本研究では、後述するライナおよび中芯の𝐶弾性係数を設計変数として繰り返し解析を行𝐶なうことで、逆問題を解くようにして得られる修正済み弾性係数を用いる。

𝑦

𝐶

𝑥
𝑧
Fig.4 Boundary conditions of macro model.

4.2.2 マクロモデル
ここでは、式(1)、(2)の周期境界条件を適用した式(3)のマクロ構造の周期境界条件を解く

5. 調整済み弾性係数の算出

ことで変位を得る。

本研究では、両面段ボールに均質化解析を

Fig.4 に試験片全体、すなわちマクロ構造を適用するため、測定された変位と解析結果のソリッド要素で簡易的にモデル化したマクロ変位が同値となるよう調整済み弾性係数を求モデルを示すが、曲げ試験に用いた試験片とめることにした。そのため、抄紙方向および同様の寸法となるようモデル化し、抄紙方向

垂直方向の境界条件を用いた解析結果の変位と垂直方向の境界条件を定義する。また、段がそれぞれ指標変位𝑇

ボールの抄紙方向および垂直方向に関しては、

ライナと中芯の弾性係数を設計空間として探マクロモデルに Fig.3 に示したミクロモデルの索する必要がある。一方で、抄紙方向と垂直

座標系をそのまま適用した。Fig.4 より、抄紙方向は異方性を有するため、同一の段ボール方向に対して解析を行う場合の境界条件は、シートから作成された試験片であれば、2 方

𝐶 および𝐶 面上の𝑥𝑦𝑧方向を完全拘束し、𝐿 線

向の解析結果の変位は設計空間上の点で交わ上の鉛直方向に荷重を与える。一方、垂直方る。そこで、抄紙方向および垂直方向に対し向に対して解析を行う場合の境界条件は、𝐶て探索を行い、ここで得られた解析結果の変および𝐶 面上の𝑥𝑦𝑧方向を完全拘束し、𝐿 線上

位をサンプル点として、抄紙方向と垂直方向の鉛直方向に荷重を加える。また、マクロモの変位‐弾性係数応答で構成される応答曲面

デルでは、Fig.3 のミクロモデルから式(2)によを作成し、この 2 方向の応答曲面の交点座標って得られる均質化弾性係数を用いて、上述を 調 整 済 み 弾 性 係 数 と し て 算 出 し た 。 Fig.5

した 2 方向の境界条件から得られた解析結果には、応答曲面による交点の断面を示していの変位と 3 章で示した指標変位𝑇、𝑇 が同る。また、段ボールの剛性に糊の影響は無視

値なることを目的としている。そのため、こできない要因であるが、ここではライナと中－ 326 －

、 𝑇 と合うように日本包装学会誌 Vol.28 No.5 (2019)

芯しかモデル化していないものの測定された

ら、抄紙方向に 10 データ、垂直方向に 10

変位が合うように弾性係数を調整しているのデータを生成した。

で、得られた弾性係数には糊の影響も含まれていると考えられる。一方で、ここで得られ6.2 階層ベイズモデル
ベイズの定理は次式に示すことができる。

た調整済み弾性係数は決定値であり、測定値のばらつきを考慮していない。そこで、次節

𝜋 𝜽|𝒙 ∝ 𝑝 𝒙|𝜽 𝜋 𝜽から調整済み弾性係数のばらつきに対する評

4

価を行う。
こ こ で 、 𝜽は 調 整 済 み 弾 性 係 数 、 𝒙は 観 測 デCorrugated medium [MPa]ー タ 、 𝜋 𝜽|𝒙 は 事 後 確 率 密 度 関 数 ( 以 下 、 事
後 分 布 ) 、 𝑝 𝒙|𝜽 は 尤 度 、 𝜋 𝜽 は 事 前 分 布 でCross directionある。ベイズの定理では、尤度と事前分布を乗ずることで事後分布を求める。事前分布は、推定以前に解析者の把握している調整済み弾

Machine direction性係数に関する類似データや過去のデータといった情報にもとづき定義することが可能でAdjust elastic modulus

ある。

Liner [MPa]

本 研究 では、 情報 量基準 WAIC を 用い てFig.5 Adjusted elastic modulus.

事前分布を推定するためのハイパーパラメー
タを 最適化する 。情報量基 準 WAIC を 活 用

6. ベイズ推定によるばらつきの評価
調整済み弾性係数は曲げ試験による測定値のばらつきを含んでいる。ここでは、ベイズ

す る に は 、 式 (4) を 次 式 に 示 す 階 層 ベ イ ズ モ
デルへ拡張しなければならない。

推定を用いて、これらのばらつきを評価する

𝜋 𝜽|𝒙, 𝜆 ∝ 𝑝 𝒙|𝜽 𝜋 𝜽|𝜆 ℎ 𝜆方法を述べる。

5

ここで、ℎ 𝜆 はハイパーパラメータである。
6.1 観測データ

このハイパーパラメータの利用によりばらつ

本研究で定めた基準荷重𝑁20 は回帰直線上の値であることから各測定値の平均となる。
そのため、各測定値の基準荷重𝑁20 に対する変位をベイズ推定に用いるための観測データ

きに関する各々の因子をモデル化することが可能であるが、全ての因子を考慮するのは事実上不可能である。ここでは、試験機による測定精度の指標をパラメータとする。

𝒙とし、ベイズ推定では少数の観測データからでも比較的良い推定結果を得られることか－ 327 －

両面段ボールに対する均質化解析とベイズ推定を用いた確率的妥当性評価技術の開発

きを含むため、このハイパーパラメータに対

6.3 調整済み弾性係数評価モデル
調整済み弾性係数は、抄紙方向および垂直して試験機の測定精度を表す指標として、式方向の応答曲面から求めた値であり、確率的(8) の 形 状 パ ラ メ ー タ 𝛼と 尺 度 パ ラ メ ー タ 𝛽をにばらつきをもつ。従って、抄紙方向と垂直

パ ラ メ ー タ に も つ ガ ン マ 分 布 𝐺𝑎 𝛼, 𝛽 で 仮 定

方向の変位‐弾性係数応答に次式の多変量正する。測定精度の指標は負の値を含まないこ規分布𝑁 𝒙|𝜽, 𝜮 を仮定し、抄紙方向および垂とから、確率密度分布が正の値のみをもつガ

直方向の多変量正規分布を乗じたものを尤度

ンマ分布を採用した。また、このパラメータ関数𝑝 𝒙|𝜽 とする。

は、零に近いほど測定精度が高いことを表す。

𝑝 𝒙|𝜽
1
1
2𝜋 2 |𝜮|ℎ 𝜆 ~𝐺𝑎 𝛼, 𝛽𝑁 𝒙|𝜽, 𝜮
1
𝒙
2

𝑒𝑥𝑝

𝜽 𝑇𝜮 1 𝒙𝜽

86
以上より、式(6)、(7)、(8)を式(5)の右辺に代入して計算した事後分布𝜋 𝜽, 𝜆|𝒙 から調整済み

こ こ で 、 𝑇は 転 置 、 𝜽は 均 質 化 解 析 に よ る 変弾性係数のばらつきに対する評価を行う。尚、

位‐弾性係数応答の各サンプル点の集合から事後分布を推定するための計算には、後述すなるベクトル、 |𝜮|は共分散行列𝜮の行列式でるマルコフ連鎖モンテカルロ法 18) を用いる。

あり、変位‐弾性係数応答のサンプル点から求められる。また、調整済み弾性係数に関す

6.4 ハイパーパラメータの最適化
情報量基準 WAIC を用いてハイパーパラ

る事前分布𝜋 𝜽|𝜆 に対しても次式の多変量正メータを最適化する方法を提案する。

規分布を仮定する。
𝜋 𝜽|𝜆𝑁 𝒙|𝜽, 𝜮𝟎 , 𝜎

𝑒𝑥𝑝

1
2𝜎

1
2𝜋𝜎

1

6.4.1 情報量基準 WAIC 1 1 ) ， 1 2 ) による評価|𝜮 |2

𝑥 , 𝑥 , 𝑥 , ⋯ 𝑥 を 独 立 な 観 測 デ ー タ 、 𝑝 𝒙|𝜽
𝒙

𝜽𝟎 𝜮

𝒙

𝜽𝟎

7

を尤度関数、𝑝 𝜽|𝜆 を調整済み弾性係数に対す る 事前 分布 と した とき 、𝑍を 正 規化 定 数とすると、WAIC は次式により求められる。

ここで、事前分布の平均𝜽 は調整済み弾性係
数 の 値 、 𝜎2 は 事 前 分 布 の 共 分 散 行 列 𝜮0 1 を 規
定 す る ハ イ パ ー パ ラ メ ー タ で あ り 、 式 (5) の

𝑇

1
𝑛

log𝑝

1
𝑍

𝑝 𝑥 |𝜽

階層ベイズモデルにおける右辺のℎ 𝜆 に対応する。このハイパーパラメータは調整済み弾性係数のばらつきを推定するための精度向上となる。調整済み弾性係数は測定値のばらつ－ 328 －

𝑝 𝑥 |𝜽 𝑝 𝜽 𝜆 𝑑𝜽

9日本包装学会誌 Vol.28 No.5 (2019)

1
𝑍

𝑉
1
𝑍

log𝑝 𝑥 |𝜽

log𝑝 𝑥 |𝜽

𝑝 𝑥 |𝜽 𝑝 𝜽 𝜆 𝑑𝜽10𝑝 𝑥 |𝜽 𝑝 𝜽 𝜆 𝑑𝜽

𝛽

𝑥̅
𝛼

14

ここで、𝑥̅ は 観測データの平均である。式(13)に式(14)を代入することで、対数尤度関

𝑊𝐴𝐼𝐶

𝑇

𝑉
𝑛

11

数を次式に書き直すことができる。

log 𝑝 𝑥 𝛼, 𝛽

こ こ で 、 𝜆は ハ イ パ ー パ ラ メ ー タ の 最 尤 推 定

𝑛 𝛼

1 log 𝑥𝑛𝛼 log 𝑥̅

𝑛 log Γ 𝛼𝑛𝛼 log 𝛼𝑛𝛼15値である。ハイパーパラメータの最適化を行
うには周辺尤度を最大化することで、WAIC

こ こ で 、 log 𝑥 は 対 数 尤 度 の 平 均 で あ る 。 式を最小化する必要がある。従って、周辺尤度(15) よ り 、 𝛼に つ い て の 尤 度 方 程 式 を 解 く こを最大化するパラメータを点推定量として求と で最 尤推定 値𝛼 を 得るこ とが できる 。従 っめる。

て 、 式 (15) を 最 大 化 す る 𝛼を 数 値 的 に 求 め る 。
こ こで は、最 尤推 定値𝛼 を 求め るため 局所 近6.4.2 ハイパーパラメータの最適化

似 最 大 化 法 1 9 ) を 用 い る 。 そ の た め 、 式 (15)ガンマ分布の周辺尤度を最大化するパラメを下記に示す局所近似式で表す。

ータの推定を行なう。ハイパーパラメータℎ 𝜆 を規定するガンマ分布の確率分布は下記

𝑐log 𝑝 𝑥 𝛼, 𝛽𝑐 𝛼𝑐 log 𝑎16のように定義される。
ただし、𝑐 は定数、𝑐 および𝑐 は次式となる。
𝑝 𝑥|𝛼, 𝛽

𝑥
𝑒𝑥𝑝
𝛤 𝛼 𝛽

𝛼
𝛽

12
𝑐

𝑐
𝛼

17

𝛼 𝑝 𝛼

18

𝑝 𝛼

観測データ𝒙が与えられたとき、式(12)の対𝑐

数尤度関数は次式のように表される。
log 𝑝 𝑥|𝛼, 𝛽𝑛 log Γ 𝑥𝛼1

また、式(16)を𝑔 𝛼 と表し、微分することで下log 𝑥𝑛𝛼 log 𝛽1𝛽記の式を得る。
𝑥13𝑔 𝛼

式(13)より𝛽に関する最尤方程式を解くこと
𝑔 𝛼により最尤推定値𝛽を得ることができる。

－ 329 －

𝑐
𝛼

𝑐
𝑐
𝛼

19

20

両面段ボールに対する均質化解析とベイズ推定を用いた確率的妥当性評価技術の開発

ここで、式(15)の対数尤度関数は上に極限をここで、𝛹 𝛼 はディガンマ関数である。式持つため𝑝 𝛼(22)に式(25)、(26)を代入して整理すると最尤て、𝑐0の 関 係 が 成 り 立 つ 。 従 っ

0と な り 𝑔 𝛼 も 同 様 に 上 に 極 限 を 持

推定値𝛼を得るための更新式を得る。

つ。すなわち、𝑔 𝛼 を最適化するには
𝑔 𝛼0を満たす𝛼を次式により求める。
𝛼

𝑐
𝑐

1
𝛼

1
𝛼

log 𝑥

21

𝛼

𝛹 𝛼

log 𝑥̅

1
𝛼

𝛹 𝛼

log 𝛼

27

本研究では、一様となるガンマ分布から更新式(21)に式(17)、(18)を代入して整理すると次さ れ た 𝛼の 推 定 値 を 𝛼 と し て 、 観 測 デ ー タ 𝑥

式を得る。
1
1
𝛼 𝛼

を 用い て更新 する ことで 最尤 推定値𝛼 の 値 を
𝑝 𝛼
𝛼 𝑝 𝛼

22

ここで、𝑝 𝑎 および𝑝 𝛼 について、式(15)の対数尤度関数におけるα log 𝛼に着目すると下に極限を持つため、𝛼 における接線は次式の関係を得る。
𝑓 𝑎

求めた。

6.5 事後分布推定
ベイズ推定によって得られた事後分布から調整済み弾性係数のばらつきを考察する。

6.5.1 推定条件
1 𝛼

log 𝛼

𝛼

事後分布の計算には多変量であることや積23𝛼

分布可能な確率分布を考慮することで、尤度式(15)に式(23)を代入すると対数尤度関数の関数の計算負荷が極めて高くなる。そこで、下限は次式の関係で表される。

マルコフ連鎖モンテカルロ法(Markov Chain
Monte Carlo Methods 、 MCMC) 1 8 ) を 用 い て 、log 𝑝 𝑥 𝛼, 𝛽𝑛𝛼 log 𝑥̅

𝑛 𝛼

1 log 𝑥𝑛𝛼 𝑙𝑜𝑔𝛼𝑛 𝛼

事 後 分 布 を 推 定 す る 。 MCMC で は 、 初 期 値

𝑛 log 𝛤 𝛼

𝛼 𝑙𝑜𝑔𝛼𝑛𝛼24の影響があるサンプル冒頭部の破棄区間
burn-in を 20000 step とし、その後の 50000

式(23)の下限が𝑑𝛼⁄𝑑𝑥

step を事後分布とした。また、収束の判定指0となる導関数𝑝 𝑎

標 と し て 𝑅 20) を 用 い た 。 𝑅 は 同 時 に 複 数 のおよび𝑝 𝛼 は次式となる。

MCMC に よ っ て 推 定 さ れ た 定 常 状 態 後 の 各
𝑝 𝑎

log 𝑥

𝛹 𝛼

log 𝑥̅

log 𝛼

25

サンプルデータを用いて得られる値である。
こ の 複 数 の MCMC に よ る 推 定 の こ と を

𝑝 𝛼

𝑛𝛹 𝛼

𝑛
𝛼

26

Chains と呼び、4Chains から以下の式より𝑅を求めた。

－ 330 －

日本包装学会誌 Vol.28 No.5 (2019)

も的確に捉えられている。このような事後分
𝑅

𝑘

1
𝑘

1𝐵
𝑘𝑊

28

布が推定できたのは、ハイパーパラメータの最適化により確信度の高い事前分布のパラメ

ここで、W は Chains 内の分散、B は Chains 間ータを得られたためであると考えられる。

の分散、𝑘は barn-in 後のサンプル数である。𝑅ま た、 Fig.7 に調整 済み 弾性係数に おけるの値が 1 に近いほど収束していることを示し

ライナと中芯の各事後分布を示す。点線内のており、値が 1.1 以下であれば定常分布に収束

区間は 95 パーセント確信区間を表わしていしているといえる。本研究では𝑅の値が 1.1 以る 。 95 パ ー セ ン ト 確 信 区 間 は 、 こ の 区 間 内

下になるように解析条件を設定している 。

に含まれる調整済み弾性係数の信頼性を示しており、区間内の値を用いることで解析結果の妥当性を保証できる。 Fig.7 より、ライ ナ

6.5.2 ばらつきの推定結果

および中芯の調整済み弾性係数の事後分布はれた調整済み弾性係数の事後分布を Fig.6 に上に凸の形状をしており、ピークが存在する。

示す。図中の白点は、5 章で得られた調整済

また、各事後分布の 95 パーセント確信区間

み弾性係数、分布の色は確率密度関数に含まれる面積も大きく、的確に両袖を捉え(Probability density function 、 PDF) を 表 している。従って、各事後分布においても確信ている。

度の高い推定結果を得られたといえる。

PDF

ベイズ推定の計算に MCMC を用いて得ら

Ad j u s t e l a s t i c m o d u l u s

Fig.6 Posterior distribution

Fig.6 よ り、事 後分布の全 体形状は、 広範
囲に分布することなく、調整済み弾性係数が含むばらつきを表しているといえる。加えて、調整済み弾性係数に対する確信度の高い区間－ 331 －

Fig.7 95% credible interval for
liner and corrugated medium.

両面段ボールに対する均質化解析とベイズ推定を用いた確率的妥当性評価技術の開発

このようにベイズ推定を用いることで、信析結果は、回帰直線と概ね一致していること

頼 性 の あ る パ ラ メ ー タ を 、 95 パ ー セ ン ト 確

から、調整済み弾性係数が正しく機能してお信区間をもって評価できることは、従来になり、均質化解析に適用可能な調整済み弾性係

い妥当性評価法である。これらのことから、数の有効性を確認することができた。一方、

本研 究で提案し た情報量基 準 WAIC を 用 い

95 パ ー セ ン ト 確 信 区 間 で は 、 上 限 値 お よ びてハイパーパラメータを最適化することで、下限値を用いて数値解析を行なった結果が、調整済み弾性係数が含むばらつきを精度よく調整済み弾性係数を用いた解析結果を内包し推定 できたと言 え、情報量 基 準 WAIC を 用て い る 。 従 っ て 、 95 パ ー セ ン ト 確 信 区 間 のいたハイパーパラメータの最適化法の有用性

調整済み弾性係数も正しく機能しており、上を示すことができた。

限値と下限値を用いて解析結果が保証できる範囲を明示的に示すことができた。

6.5.3 妥当性検証
調整済み弾性係数の有効性を検証するため、推定された事後分布の平均値と 95 パーセン
ト確信区間を用いた均質化解析の結果と測定
値および測定値の回帰直線と比較を行った。
調 整 済 み 弾 性 係 数 の 上 限 値 𝑈お よ び 下 限 値 𝐿はそれぞれ次式により求めた。
𝑈

𝐿

𝜇̅

𝜇̅

1.96

1.96

𝑠
√𝑛
𝑠
√𝑛

29
(a) Machine direction30

ここで、𝜇̅ は 事後分布の平均、sは事後分布の標 準 偏 差 、𝑛は 観 測 デ ータ 数 で あ る。 ま た 、
この上限値と下限値を用いた数値解析は、信頼性が保証できる解析結果の範囲を意味する。
Fig.8 に抄紙方向および垂直方向の比較結 果を示す。図中の凡例以外の各線は Fig.2 の拡
大 図 で 示 し た 各 測 定 値 (Measuring values) 、点線は測定値から求めた回帰直線を表している。 Fig.8 より調整済み弾性係数を用いた 解－ 332 －(b) Cross direction
Fig.8 Validation of estimated elastic
modulus. Each line shows measurement values in each experiment .

日本包装学会誌 Vol.28 No.5 (2019)

7. 結言

下限値を均質化解析に用いた結果に調

本論文では、均質化法の数学的な理論の拡

整済み弾性係数を用いた解析結果が内張を行わずに、両面段ボール対して均質化解包 さ れ て い る こ と か ら 、 95 パ ー セ ン ト

析の適用が可能となる調整済み弾性係数の提確信区間が正しく機能していることが案および調整済み弾性係数が含むばらつきの確認でき、解析結果の信頼性が保証で評価法を提案した。調整済み弾性係数においきる範囲を明示的に示した。

ては、試験と均質化解析から導けることを示

本研究で示した一連の手法により、両面段した。また、ベイズ推定を用いた調整済み弾

ボールに対する均質化解析が可能となり、且性係数におけるばらつきの評価については、

つその妥当性を保障することで、解析に基づ

ハイパーパラメータに対して情報量基準

く包装設計の信頼性向上に役立てることがで

WAIC を用いて最適化を行ったパラメータかきる。本論では線形範囲に限定して議論を進ら推定された事後分布より、調整済み弾性係

めているが、非線形解析に拡張することで、数のばらつきが精度よく評価できることを示

今まで、試行錯誤に依存していた緩衝材等のした。得られた知見を以下にまとめる。

設計においても均質化解析を活用した信頼性(1) 基準荷重における指標変位に弾性係数が評価を伴う設計に寄与できるものと考えられ合うように均質解析に行うことで、両る。

面段ボールに均質化解析が適用可能な調整済み弾性係数を算出することがで＜参考文献＞きた。

1)

松島理、松島成夫、面圧縮を受ける段ボ(2) 調整済み弾性係数は指標変位に合わせる

ール中芯の形状と弾性強度(楕円対象曲

ことで、ライナと中芯を接着する糊な線 部 よ り な る 中 芯 波 形 の 場 合 ) 、 Vol.5,

どの影響を考慮した値として算出される。

No.2, pp.107-118 (1996)(3) 情報量基準 WAIC を用いてハイパーパ

2)

中 川 幸 臣 、 丹 羽 一 邦 、 包 装 設 計 の FEM

ラメータを最適化することで、確信度解析における段ボールのモデル化の検の高い事前分布を得ることができた。

討 、 日 本 包 装 学 会 誌 、 Vol.14, No.5,

その結果、推定された事後分布が調整

pp.317-327 (2005)

済み弾性係数のばらつきを捉えている3)

中川幸臣、丹羽一邦、斎藤勝彦、段ボーことから、最適化パラメータによってル構造体の落下衝撃シミュレーション、

調整済み弾性係数のばらつきを精度よ日本包装学会誌、Vol.15, No.6, pp.333-

く推定できることを明らかにした。

343 (2006)(4) 調整済み弾性係数を用いた均質化解析と4)

三上英一、茨木孝昌、小林孝男、小高功、

測定値および測定値の回帰直線におい有限要素法による段ボールシートの垂て、概ね一致する結果を示した。また、

直圧縮強さの解析、日本包装学会誌、

95 パ ー セ ン ト 確 信 区 間 の 上 限 値 お よ び

Vol.13, No.4, pp.243-252 (2004)－ 333 －

両面段ボールに対する均質化解析とベイズ推定を用いた確率的妥当性評価技術の開発

5)
6)

寺田賢二郎、菊池昇、均質化法入門、

corrugated

丸善出版(2003)

Engineering Journal, Vol.3, No.4(2016)

高野直樹、座古勝、菊池昇、均質化法にMechanical14) 高野直樹、均質化法による新しい数値シ

よる単層波状コアサンドイッチ板の解ミュレーション、日本複合材料学会誌、

析 、 日 本 材 料 学 会 、 Vol.44, No.502,

Vol.7, No.1, pp.4-11(2001)15) H. Anders and R. Matti, Large strain

pp.905-910(1995)

elasto-plastic

7) R. V. Kohn and M. Vogelius, A new

model

of

paper

and

rapidly

corrugate board, International Journal ofvarying thickness, International Journal

Solids and Structures, Vol.45, pp.3334-

of Solids and Structure, Vol.20, No.4,

3352(2008)

model

for

thin

plateswith16) Z. W. Wang, Y. P. E: Energy absorption

pp.333-350(1984)

properties of multi-layered corrugated8) 寺田賢二郎、平山紀夫、山本晃司、松原

9)

cardboard,

成志朗、複合板のマルチスケール解析のpaperboardinための数値平板試験、日本計算工学会、

humidities,

Materials

No.20150001(2015)

Vol.32, pp.3476-3485(2011)

ベイズ統計学概説、松原望、培風館、

variousambient

andDesign,17) 奥村大、大野信忠、小野田孝行、曽我将(2010)

孝、弾性積層段ボールの圧縮強度に微

10) M. B. Christopher, Pattern recognition

視構造の及ぼす影響の均質化解析、計and machine learning, Springer, pp.359-

算 数 理 工 学 論 文 集 、 Vol.12 、 No. 18-378(2006)

121212(2012)11) S. Watanabe, Asymptotic equivalence18) 久保裕也、データ解析のための統計モデ

of Bayes cross validation and widely

リ ン グ 入 門 、 岩 波 書 店 、 pp.170-191

applicable(2012)

information

criterion

inJournalof

singularlearning

theory,MachineLearning

Research,

19)

P.

Minka,

Distribution,

Vol.14,

Estimatea

Microsoft

Gamma
Research

Technical Report (2002)

pp.3571-3591(2010)12)T.

likelihood20) A. Gelman, J. B. Carlin, H. S. Stern, D.

Identification, of Gaussian autoregress-

B. Dunson, A. Vehtari and D. B. Rubin,

ive moving average model, Biometrika,

Bayesian data analysis, Journal of theVol.60, No.2, pp.255-265(1973)

AmericanH.

Akaike,Maximum13) M. Shimamura, Y. Ezawa, Y. Tamura, S.
Takashimizu

and

D.

Satou,

Statistical

Association,

Vol.109, No.507, p.297(2014)

High

accurate analysis by experimental and（原稿受付2019 年 4 月 10 日）

simulation using Bayesian inference for（審査受理2019 年 8 月 21 日）－ 334 －