試行回数と信頼性

ロボットの確率・統計　第５章

以下の２つのグラフを再現したい

clear; close all;

ts = 0:0.01:1;
as = [1 0 0 1 1];
bs = [1 1 1 1 1];

p_t_a1_i_1 = zeros(length(as), length(ts));
p_t_a1_i_1(1,:) = 1/length(ts);

ps_a = compute_posteriors(ts, as, p_t_a1_i_1);
ps_b = compute_posteriors(ts, bs, p_t_a1_i_1);
for idx = 1:length(as)
    subplot(1,length(as),idx);
    plot(ts,ps_a(idx,:));
    ylim([0 0.06]);
    hold on;
    plot(ts,ps_b(idx,:));
    ylim([0 0.06]);
    p = prob_b_is_greather_than_a(idx, ts, ps_a, ps_b);
    title(sprintf('Pr(b<a)=%.2f',p));
end

big;

function p = prob_b_is_greather_than_a(n,ts,ps_a,ps_b)
p = 0;
a = ps_a(n,:);
b = ps_b(n,:);
for idx=2:length(ts)
    if a(idx) < b(idx)
        p = p + a(idx);
    end
end
end


function posteriors =  compute_posteriors(ts, as, p_t_a1_i_1)
qts = zeros(length(as), length(ts));
posteriors = zeros(length(as), length(ts));

for idx = 1:length(as)
    a = as(idx);
    if a==0 %失敗
        p_ai_t = 1-ts;
    elseif a == 1 %成功
        p_ai_t = ts;
    end
    qts(idx,:) = p_ai_t .* p_t_a1_i_1(idx,:);
    eta = 1/sum(qts(idx,:));
    posteriors(idx,:) = eta * qts(idx,:);
    p_t_a1_i_1(idx+1,:) = posteriors(idx,:);
end

end

同じく以下のグラフを再現したい

clear; close all;

ts = 0:0.01:1;
Ns = 5:5:100;

success_rates = [0.7 0.8 0.9 1.0];
for jdx = 1:length(success_rates)
    success_rate = success_rates(jdx);
    pN = [];
    for N=Ns
        as = generate_sample(N,0.6);
        bs = generate_sample(N,success_rate);
        % as = [1 0 0 1 1];
        % bs = [1 1 1 1 1];

        p_t_a1_i_1 = zeros(length(as), length(ts));
        p_t_a1_i_1(1,:) = 1/length(ts);

        ps_a = compute_posteriors(ts, as, p_t_a1_i_1);
        ps_b = compute_posteriors(ts, bs, p_t_a1_i_1);
        ps = [];
        for idx = 1:length(as)
            p = prob_b_is_greather_than_a(idx, ts, ps_a, ps_b);
            ps(end+1) = p;
        end
        pN(end+1) = ps(end);
    end
    plot(Ns, pN, 'o-');
    hold on;
end
yline(0.05);


big;

function xs = generate_sample(N, success_rate)
successes = round(success_rate * N);
failures = N - successes;
xs = horzcat(ones(1,successes), zeros(1,failures));
end


function p = prob_b_is_greather_than_a(n,ts,ps_a,ps_b)
p = 0;
a = ps_a(n,:);
b = ps_b(n,:);
for idx=2:length(ts)
    if a(idx) < b(idx)
        p = p + a(idx);
    end
end
end


function posteriors =  compute_posteriors(ts, as, p_t_a1_i_1)
qts = zeros(length(as), length(ts));
posteriors = zeros(length(as), length(ts));

for idx = 1:length(as)
    a = as(idx);
    if a==0 %失敗
        p_ai_t = 1-ts;
    elseif a == 1 %成功
        p_ai_t = ts;
    end
    qts(idx,:) = p_ai_t .* p_t_a1_i_1(idx,:);
    eta = 1/sum(qts(idx,:));
    posteriors(idx,:) = eta * qts(idx,:);
    p_t_a1_i_1(idx+1,:) = posteriors(idx,:);
end

end

少し改良した

clear; close all;

ts = 0:0.01:1;
Ns = 5:5:100;

success_rates = [0.7 0.8 0.9 1.0];
pNs = zeros(length(success_rates),length(Ns));
for jdx = 1:length(success_rates)
    success_rate = success_rates(jdx);
    pN = [];
    for N=Ns
        as = generate_sample(N,0.6);
        bs = generate_sample(N,success_rate);

        p_t_a1_i_1 = zeros(length(as), length(ts));
        p_t_a1_i_1(1,:) = 1/length(ts);

        ps_a = compute_posteriors(ts, as, p_t_a1_i_1);
        ps_b = compute_posteriors(ts, bs, p_t_a1_i_1);
        ps = [];

        subplot(length(success_rates), 1, jdx);
        cla;
        plot(ts,ps_a(N,:));
        hold on;
        plot(ts,ps_b(N,:));
        ylim([0 0.2]);
        title(sprintf('試行回数: %d', N)); 
        legend(['before=60(％)'],[sprintf('after=%d(％)',success_rate*100)]);
        big;
        drawnow;       
        pause(0.001);

        for idx = 1:length(as)
            p = prob_b_is_greather_than_a(idx, ts, ps_a, ps_b);
            ps(end+1) = p;
        end
        pN(end+1) = ps(end);
    end
    pNs(jdx,:) = pN;
end
figure();
plot(Ns, pNs, 'o-');

yline(0.05);


big;

function xs = generate_sample(N, success_rate)
successes = round(success_rate * N);
failures = N - successes;
xs = horzcat(ones(1,successes), zeros(1,failures));
xs = xs(randperm(length(xs))); %shuffle
end


function p = prob_b_is_greather_than_a(n,ts,ps_a,ps_b)
p = 0;
a = ps_a(n,:);
b = ps_b(n,:);
for idx=2:length(ts)
    if a(idx) < b(idx)
        p = p + a(idx);
    end
end
end


function posteriors =  compute_posteriors(ts, as, p_t_a1_i_1)
qts = zeros(length(as), length(ts));
posteriors = zeros(length(as), length(ts));

for idx = 1:length(as)
    a = as(idx);
    if a==0 %失敗
        p_ai_t = 1-ts;
    elseif a == 1 %成功
        p_ai_t = ts;
    end
    qts(idx,:) = p_ai_t .* p_t_a1_i_1(idx,:);
    eta = 1/sum(qts(idx,:));
    posteriors(idx,:) = eta * qts(idx,:);
    p_t_a1_i_1(idx+1,:) = posteriors(idx,:);
end

end