計数規準型一回抜取検査(JIS Z 9002)の拡張

#計数規準型一回抜取検査(JIS Z 9002)の拡張

###ざっくり説明
拡張した部分は、

• JISで使われている二項分布以外に、ポアソン分布、超幾何分布を使えるようにしたこと
• JISよりも、$\alpha$と$\beta$をもうちょっと厳密に守るようにしたこと

です。

###用語など
一般的に品質保証/品質管理分野で通用する用語がわからないので、こんな意味で使ってます、という定義を書いておきます。

• $\alpha$　生産者危険、本来合格になるはずの不適合品率のロットが不合格とされてしまう確率。　
• $\beta$　消費者危険、本来不合格となるはずの不適合品率のロットが合格とされてしまう確率。
• $AQL$　なるべく合格としたい不適合率の上限。$p_0$のこと
• $RQL$　なるべく不合格としたい不適合率の下限。　$p_1$のこと、なお$AQL<RQL$です。
• $n$ サンプルサイズ
• $c$　合格判定個数
• $x$　不適合品の数　

###アルゴリズム
与えられた$\alpha$、$\beta$、$AQL$、$RQL$のもとで、次の式を同時に満たす$n$,$c$を探索的に求めます。

$L(AQL)$を不適合品率$AQL$のロットが合格する確率、$L(RQL)$を不適合品率$RQL$のロットが合格する確率、$P(x)$をサイズ$n$のサンプル中に不適合品が$x$個存在する確率として

\alpha = 1-L(AQL)=\sum_{x=c+1}^{n}P(x)=1-\sum_{x=0}^{c}P(x)\\
\beta = L(RQL)=\sum_{x=0}^{c}P(x)

$n$について小さい方から大きい方へ探し、最初に見つかった$n$,$c$を抜き取り方式とします。

これをjavaで書いてみました。
apache commons math3を使っています。

SamplingInspectionByAttribute.java

package QA;

import java.util.ArrayList;
import java.util.List;
import java.util.TreeMap;
import org.apache.commons.math3.distribution.AbstractIntegerDistribution;
import org.apache.commons.math3.distribution.BinomialDistribution;
import org.apache.commons.math3.distribution.HypergeometricDistribution;
import org.apache.commons.math3.distribution.PoissonDistribution;

/**
 * JIS Z 9002(計数規準型抜取検査)の拡張
 *
 * @author gsgs
 */
public class SamplingInspectionByAttribute {

    double AQL; //合格としたい品質水準
    double RQL; //不合格としたい品質水準

    double alpha; //生産者リスク
    double beta; //消費者リスク

    int max_sample; //最大サンプル数
    int interval;   //n探索の細かさ

    DistributionType distribution;

    int n; //最小サンプル数
    int c; //許容基準

    boolean fin = false;

    /**
     * 分布<br>
     * BINOMIAL 二項分布　JISに同じ。通常はこれを使う<br>
     * POISSON ポアソン分布　拡張。非常に大きいロットサイズかつ非常に低い不良率の時に使う<br>
     * HYPERGEOMETRIC　超幾何分布　拡張。有限母集団からの非復元抽出の時に使う<br>
     */
    public static enum DistributionType {
        /**
         * 二項分布
         */
        BINOMIAL,
        /**
         * ポアソン分布
         */
        POISSON,
        /**
         * 超幾何分布
         */
        HYPERGEOMETRIC
    }

    /**
     * JIS Z 9002タイプの計数規準型抜取検査<br>
     * α=0.05,β=0.1,二項分布を用いて必要サンプル数(最小値)と許容範囲を求めます。
     *
     * @param AQL なるべく合格としたい最大不良率(0.0-1.0で与える)　JISではp0
     * @param RQL なるべく不合格としたい最小不良率(0.0-1.0で与える) JISではp1
     */
    public SamplingInspectionByAttribute(double AQL, double RQL) {
        this(AQL, RQL, 10000, DistributionType.BINOMIAL);
    }

    /**
     * JIS Z 9002タイプの計数規準型抜取検査<br>
     * α=0.05,β=0.1,指定した分布とロットサイズを用いて必要サンプル数(最小値)と許容範囲を求めます。
     *
     * @param AQL なるべく合格としたい最大不良率(0.0-1.0で与える)　JISではp0
     * @param RQL なるべく不合格としたい最小不良率(0.0-1.0で与える) JISではp1
     * @param max_sample　サンプルサイズ。超幾何分布の時は厳密なロットサイズを与える(他の分布の時には必要サンプル数の探索上限を与える)
     * @param distribution　二項分布、ポアソン分布、超幾何分布から選ぶ
     */
    public SamplingInspectionByAttribute(double AQL, double RQL, int max_sample, DistributionType distribution) {
        this(AQL, RQL, 0.05, 0.1, max_sample, 1, distribution);
    }

    /**
     * JIS Z 9002タイプの計数規準型抜取検査<br>
     * 指定したパラメータ、指定した探索粒度で必要サンプル数(最小値ではないかもしれない)と許容範囲を求めます。
     *
     * @param AQL なるべく合格としたい最大不良率(0.0-1.0で与える)　JISではp0
     * @param RQL なるべく不合格としたい最小不良率(0.0-1.0で与える) JISではp1
     * @param alpha　生産者危険　0.0-1.0で与える
     * @param beta　消費者危険　0.0-1.0で与える
     * @param max_sample　サンプルサイズ。超幾何分布の時は厳密なロットサイズを与える(他の分布の時には必要サンプル数の探索上限を与える)
     * @param interval
     * 探索粒度。1以上かつmax_sample以下で与える。大きい数字ほど、ざっくり調べる(2以上だと必要サンプル数が最小でなくなる可能性がある、大きすぎると全数検査が必要とされやすくなる)
     * @param distribution　二項分布、ポアソン分布、超幾何分布から選ぶ
     */
    public SamplingInspectionByAttribute(double AQL, double RQL, double alpha, double beta, int max_sample, int interval, DistributionType distribution) {
        this.AQL = AQL;
        this.RQL = RQL;
        this.alpha = alpha;
        this.beta = beta;
        this.max_sample = max_sample;
        this.interval = interval;
        this.interval = this.interval < 1 ? 1 : this.interval;
        this.distribution = distribution;

        fin = false;

        for (n = interval; n < max_sample && !fin; n += interval) {
            IntegerDistribution a_dist = new IntegerDistribution(max_sample, n, AQL, distribution);
            IntegerDistribution b_dist = new IntegerDistribution(max_sample, n, RQL, distribution);

            int min_c = a_dist.inverseCumulativeProbability(1 - alpha);
            int max_c = b_dist.inverseCumulativeProbability(beta) - 1;
            if (min_c <= max_c && a_dist.cumulativeProbability(min_c) >= 1 - alpha && b_dist.cumulativeProbability(min_c) <= beta) {
                c = min_c;
                fin = true;
                break;
            }
        }
    }

    /**
     * 全数検査が必要かどうかを返します
     *
     * @return 全数検査が必要なときtrueを返す
     */
    public boolean needTotalInspection() {
        return !fin;
    }

    /**
     * 計算された最小抜取サンプル数を返す
     *
     * @return
     */
    public int getN() {
        return fin ? n : max_sample;
    }

    /**
     * 計算された許容規準を返す。Ac
     *
     * @return
     */
    public int getC() {
        return fin ? c : (int) (AQL * max_sample);
    }

    /**
     * サンプリングする対象番号の一例を作成する。
     *
     * @return サンプリングする製品の番号のリスト
     */
    public Set<Integer> createSampleList() {
        Set<Integer> sampleList = new TreeSet<>();
        while (sampleList.size() < n) {
            int x = (int) (Math.random() * max_sample);
            sampleList.add(x + 1);
        }
        return sampleList;
    }

    /**
     * 抜取検査計画をテキストで返す
     *
     * @return
     */
    public String getResultReport() {
        IntegerDistribution a_dist = new IntegerDistribution(max_sample, n, AQL, distribution);
        IntegerDistribution b_dist = new IntegerDistribution(max_sample, n, RQL, distribution);
        StringBuilder buf = new StringBuilder();
        buf.append("計数抜取検査の計画\n")
            .append("\n")
            .append("\n")
            .append("測定タイプ: 合格/不合格\n")
            .append("不良率(%)によるロット品質\n")
            .append("ロットサイズ: ").append(max_sample).append("\n")
            .append("指定された分布: ").append(distribution.toString()).append(" Distribution\n")
            .append("\n")
            .append("方法\n")
            .append("\n")
            .append("合格品質水準 (AQL): ").append(AQL * 100).append("%\n")
            .append("生産者リスク (α): ").append(alpha * 100).append("%\n")
            .append("\n")
            .append("不合格となる品質水準 (RQLまたはLTPD): ").append(RQL * 100).append("%\n")
            .append("消費者リスク (β): ").append(beta * 100).append("%\n")
            .append("\n")
            .append("生成された計画\n")
            .append("\n");

        if (fin) {
            buf.append("サンプルサイズ: ").append(n).append("\n")
                .append("許容数: ").append(c).append("\n")
                .append("\n")
                .append("不良率(%)\t合格確率\n")
                .append(AQL * 100).append("\t").append(String.format("%.4f", a_dist.cumulativeProbability(c))).append("\n")
                .append(RQL * 100).append("\t").append(String.format("%.4f", b_dist.cumulativeProbability(c))).append("\n")
                .append("\n");
        } else {
            buf.append("全数検査\n")
                .append("\n");
        }
        return buf.toString();
    }

    /**
     * OC曲線を返す
     *
     * @param min　OC曲線のx軸(ロット内の不良率)の最小値を与える 0.0-1.0
     * @param max　OC曲線のx軸(ロット内の不良率)の最大値を与える 0.0-1.0
     * @param interval　OC曲線のロット内の不良率の計算間隔を与える 0.0-1.0
     * @return double[OC曲線の点インデクス][x値=0,y値=1]
     */
    public double[][] getOCCurve(double min, double max, double interval) {
        int n = (int) ((max - min) / interval) + 1;
        double[][] ret = new double[n][2];

        if (needTotalInspection()) {
            for (int i = 0; i < n; i++) {
                ret[i][0] = min + interval * i;
                ret[i][1] = ret[i][0] <= AQL ? 1 : 0;
            }
        } else {
            for (int i = 0; i < n; i++) {
                ret[i][0] = min + interval * i;
                if (min + interval * i == 0d) {
                    ret[i][1] = 1.0d;
                } else {
                    IntegerDistribution dist = new IntegerDistribution(max_sample, getN(), ret[i][0], distribution);
                    ret[i][1] = dist.cumulativeProbability(getC());
                }
            }
        }
        return ret;
    }

    /**
     * このクラス内で使用する統計分布の抽象的なオブジェクト
     */
    public static class IntegerDistribution {

        AbstractIntegerDistribution dist;

        public IntegerDistribution(int samples, int trial, double prob, DistributionType distribution) {
            switch (distribution) {
                case BINOMIAL:
                    dist = new BinomialDistribution(trial, prob);
                    break;
                case POISSON:
                    dist = new PoissonDistribution(prob * trial);
                    break;
                case HYPERGEOMETRIC:
                    dist = new HypergeometricDistribution(samples, (int) (samples * prob), trial);
                    break;
                default:
                    throw new AssertionError();
            }
        }

        public double cumulativeProbability(int x) {
            return dist.cumulativeProbability(x);
        }

        public int inverseCumulativeProbability(double p) {
            return dist.inverseCumulativeProbability(p);
        }
    }
}

###使い方

内部監査などで、例えば書類としてある品質記録10000枚に問題ありそうかなさそうかをまずは判断しなければいけないとき、あるベータのもとでの抜き取り検査方式をその場でリアルタイムに設計・計画できれば素早く詳細点検の必要性を判断できるのでは？と考えて、パラメタを比較的自由に設定できるように書きました。
(そのため、抜き取りNo.があれば便利だろうと、余計なおせっかいをやっています。)

動機とは異なりますが、普通の使い方での一例を次の「使い方」に書いておきます。

注意点として・・
createSampleList()を使って抜き取る試料番号を決定する場合は、max_sampleは確率分布の種類を問わず正確に与えてください。
(二項分布とポアソン分布で抜き取り検査方式が存在するかどうかだけのチェックならば、max_sampleは探索するnの上限で構いません)

使い方.java

    public static void main(String[] args) {
        //条件設定
        double AQL = 0.005;
        double RQL = 0.05;
        double alpha = 0.05;
        double beta = 0.1;
        int max_sample = 1000;
        int interval = 1;
        DistributionType distribution = SamplingInspectionByAttribute.DistributionType.BINOMIAL;

        //インスタンス生成
        SamplingInspectionByAttribute m = new SamplingInspectionByAttribute(
            AQL, 
            RQL, 
            alpha, 
            beta, 
            max_sample, 
            interval, 
            distribution
        );
        
        //結果の表示
        System.out.println(m.getResultReport());

        //OC曲線(本当はチャートライブラリへ突っ込んで
        System.out.println("\nOC曲線");
        double[][] occ = m.getOCCurve(0, RQL * 1.1, RQL / 50);
        for (double[] occ1 : occ) {
            System.out.printf("%.4f\t%.4f%n", occ1[0] * 100, occ1[1] * 100);
        }

        //便利機能：の抜き取るサンプル番号(乱数で生成、ソート済み)
        System.out.println("\n抜き取りNo\n");
        Set<Integer> sampleList = m.createSampleList();
        int cnt = 0;
        for (Integer samp : sampleList) {
            System.out.print("□ " + samp + "\t");
            cnt++;
            if (cnt % 5 == 0) {
                System.out.println("");
            }
        }
        
        System.out.println("------");
    }

###実行結果

run:
計数抜取検査の計画


測定タイプ: 合格/不合格
不良率(%)によるロット品質
ロットサイズ: 1000
指定された分布: BINOMIAL Distribution

方法

合格品質水準 (AQL): 0.5%
生産者リスク (α): 5.0%

不合格となる品質水準 (RQLまたはLTPD): 5.0%
消費者リスク (β): 10.0%

生成された計画

サンプルサイズ: 105
許容数: 2

不良率(%)	合格確率
0.5	0.9839
5.0	0.0992



OC曲線
0.0000	100.0000
0.1000	99.9826
0.2000	99.8712
0.3000	99.5970
0.4000	99.1141
0.5000	98.3947
0.6000	97.4253
0.7000	96.2039
0.8000	94.7366
0.9000	93.0363
1.0000	91.1201
1.1000	89.0084
1.2000	86.7235
1.3000	84.2887
1.4000	81.7276
1.5000	79.0632
1.6000	76.3182
1.7000	73.5137
1.8000	70.6698
1.9000	67.8051
2.0000	64.9366
2.1000	62.0796
2.2000	59.2480
2.3000	56.4540
2.4000	53.7081
2.5000	51.0198
2.6000	48.3967
2.7000	45.8454
2.8000	43.3714
2.9000	40.9789
3.0000	38.6710
3.1000	36.4503
3.2000	34.3182
3.3000	32.2756
3.4000	30.3225
3.5000	28.4587
3.6000	26.6832
3.7000	24.9946
3.8000	23.3914
3.9000	21.8715
4.0000	20.4328
4.1000	19.0727
4.2000	17.7887
4.3000	16.5781
4.4000	15.4381
4.5000	14.3657
4.6000	13.3582
4.7000	12.4126
4.8000	11.5261
4.9000	10.6957
5.0000	9.9187
5.1000	9.1923
5.2000	8.5139
5.3000	7.8807
5.4000	7.2904
5.5000	6.7404

抜き取りNo

□ 2	□ 8	□ 12	□ 40	□ 45	
□ 57	□ 58	□ 66	□ 74	□ 85	
□ 87	□ 96	□ 107	□ 113	□ 120	
□ 122	□ 134	□ 142	□ 146	□ 150	
□ 163	□ 172	□ 175	□ 186	□ 205	
□ 207	□ 213	□ 237	□ 247	□ 251	
□ 253	□ 261	□ 276	□ 278	□ 286	
□ 293	□ 306	□ 307	□ 310	□ 316	
□ 319	□ 320	□ 337	□ 366	□ 374	
□ 385	□ 392	□ 407	□ 419	□ 421	
□ 429	□ 457	□ 471	□ 479	□ 518	
□ 528	□ 533	□ 546	□ 564	□ 572	
□ 586	□ 588	□ 598	□ 599	□ 611	
□ 617	□ 628	□ 629	□ 649	□ 650	
□ 674	□ 677	□ 683	□ 689	□ 691	
□ 692	□ 710	□ 713	□ 722	□ 739	
□ 744	□ 752	□ 756	□ 761	□ 794	
□ 807	□ 809	□ 830	□ 843	□ 863	
□ 881	□ 888	□ 889	□ 898	□ 902	
□ 906	□ 913	□ 917	□ 931	□ 949	
□ 964	□ 984	□ 993	□ 995	□ 1000	
------
ビルド成功(合計時間: 0秒)

###お願い
実用とするにはGUIが必要ですので、お好きにかぶせてください。
間違いに気づかれた方は、どうか、教えてください。
n,cを探索せずとも、一発で求められるよ、という方もぜひ教えてください。