function [dp, exp_dp] = dprime_lik(data, figon)
% Usage: 
% dp: ふつうに計算したd'
% exp_dp: 尤度から計算したp(hit), p(fa)に基づいたd'
% data: [hit miss fa cr]のtrial数
% figon: 1のときはfigureを描画
% 
% example
% [dp, exp_dp] = dprime_lik([10 0 5 5], 1)

hit  = data(1);
miss = data(2);
fa   = data(3);
cr   = data(4);

ph = (0:hit + miss)/(hit + miss);
pf = (0:fa + cr)  / (fa  + cr);

likh = ph.^hit .* (1-ph).^miss;
likf = pf.^fa .* (1-pf).^cr;
% normalized
nlikh = likh / sum(likh(:));
nlikf = likf / sum(likf(:));

% dataから直で計算したp(hit), p(fa)
dat_ph = hit / (hit + miss);
dat_pf = fa  / (fa  + cr);

% 尤度使って期待値を推定した p(hit), p(fa)
exp_ph = ph * nlikh';
exp_pf = pf * nlikf';

% Inf, -Infの処理
if dat_ph == 0 || dat_pf == 1
    dp = -Inf;
elseif dat_ph == 1 || dat_pf == 0
    dp = Inf;
else
    dp = norminv(dat_ph) -norminv(dat_pf);
end

% こっちはInf, -Infにはならないのでそのまま
exp_dp = norminv(exp_ph) -norminv(exp_pf);

if figon
    % figonでだけmeshgridを作る
    lik = likh' * likf;
    
    binh1 = (hit + miss) / 10;
    binh2 = (fa + cr) / 10;
    
    % exp_dp = norminv(phd) -norminv(pfd);
    % この式を使って、exp_dp が既知のときの
    % pfd=0:0.01:1に振ったときのphdの値を計算してやる。
    % norminv(phd) = exp_dp + norminv(pfd);
    % phd = normcdf(exp_dp + norminv(pfd), 0, 1);
    invpfd = -4:0.1:4;
    phd = normcdf(exp_dp + invpfd, 0, 1);
    pfd = normcdf(invpfd, 0, 1);
    phd = [0 phd 1]; 
    pfd = [0 pfd 1];

    figure;
    set(gcf, 'PaperOrientation', 'landscape')
    set(gcf, 'Paperposition', [2 2 12 5])
    
    subplot(2,2,1)
    hold on
    plot(ph, nlikh, '.-')
    axis tight
    yl1 = get(gca,'ylim');
    plot([exp_ph exp_ph], [0 yl1(2)], 'r:')
    ss1 = sprintf('hit %3d / miss %3d / dat %0.2f / exp %0.2f',...
        hit, miss, dat_ph, exp_ph);
    title(ss1)
    ww = get(gca, 'XLabel'); set(ww, 'string', 'p(hit)')
    set(gca, 'ytick', yl1)
    set(gca, 'yticklabel',...
        [{0}, {sprintf('%0.2f', yl1(2))}])

    subplot(2,2,3)
    hold on
    plot(pf, nlikf, '.-')
    axis tight
    yl2 = get(gca,'ylim');
    plot([exp_pf exp_pf], [0 yl2(2)], 'r:')
    ss2 = sprintf('FA %3d / CR %3d / dat %0.2f / exp %0.2f',...
        fa, cr, dat_pf, exp_pf);
    title(ss2)
    ww = get(gca, 'XLabel'); set(ww, 'string', 'p(FA)')
    set(gca, 'ytick', yl2)
    set(gca, 'yticklabel',...
        [{0}, {sprintf('%0.2f', yl2(2))}])

    subplot(1,2,2)
    imagesc(lik)
    axis xy
    axis tight
    axis square
    set(gca,'xtick', 1:binh2:length(pf))
    set(gca,'xticklabel', 0:0.1:1)
    set(gca,'ytick', 1:binh1:length(ph))
    set(gca,'yticklabel', 0:0.1:1)
    hold on
    xx = pfd * (length(pf) - 1) + 1;
    yy = phd * (length(ph) - 1) + 1;
    plot(xx, yy, 'w-');

    xx0 = exp_pf * (length(pf) - 1) + 1;
    yy0 = exp_ph * (length(ph) - 1) + 1;
    plot(xx0, yy0, 'wo');
    
    % plot(1, 1, 'wo') -- [0,0]
    % plot(length(pf), length(ph), 'wo') -- [1,1]

    ww = get(gca, 'XLabel'); set(ww, 'string', 'p(FA)')
    ww = get(gca, 'YLabel'); set(ww, 'string', 'p(hit)')

    ss3 = sprintf('d-prime(dat) %1.2f / d-prime(exp) %1.2f',...
    dp, exp_dp);
    title(ss3)
    
    % print -dpng s_dprime_lik1

end

return

