tuse computed percentile values - cosmo - front and backend for Markov-Chain Monte Carlo inversion of cosmogenic nuclide concentrations
 (HTM) git clone git://src.adamsgaard.dk/cosmo
 (DIR) Log
 (DIR) Files
 (DIR) Refs
 (DIR) README
 (DIR) LICENSE
       ---
 (DIR) commit a26ef7c98f2bbcebeea1c38d677a7cd508b27f12
 (DIR) parent d4be54b2ee3e227b7fd3ec6d95ffaa1ace77dbd4
 (HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
       Date:   Wed, 18 Nov 2015 19:55:10 +0100
       
       use computed percentile values
       
       Diffstat:
         M matlab/generate_plots.m             |      53 +++++++++++++++++++++++--------
       
       1 file changed, 39 insertions(+), 14 deletions(-)
       ---
 (DIR) diff --git a/matlab/generate_plots.m b/matlab/generate_plots.m
       t@@ -14,7 +14,21 @@ titles = 0;
        
        % Burn-in and convergence overview MCMC plot, fh(1)
        fh(1) = figure('visible', show_figures);
       -    
       +
       +% generate matrices for percentiles
       +epsilon_int_25 = zeros(Nwalkers,1);
       +epsilon_int_50 = zeros(Nwalkers,1);
       +epsilon_int_75 = zeros(Nwalkers,1);
       +epsilon_gla_25 = zeros(Nwalkers,1);
       +epsilon_gla_50 = zeros(Nwalkers,1);
       +epsilon_gla_75 = zeros(Nwalkers,1);
       +record_threshold_25 = zeros(Nwalkers,1);
       +record_threshold_50 = zeros(Nwalkers,1);
       +record_threshold_75 = zeros(Nwalkers,1);
       +E_25 = zeros(Nwalkers,1);
       +E_50 = zeros(Nwalkers,1);
       +E_75 = zeros(Nwalkers,1);
       +
        for iwalk=1:min(4,Nwalkers) %only up to the first four walks
          QsBurnIns(:,iwalk)=Ss{iwalk}.QsBurnIn;
          Qss(:,iwalk)       =Ss{iwalk}.Qs;
       t@@ -222,28 +236,38 @@ for i1 = 1:M % for each model parameter
            
            if i1 == 1
                title(['MCMC walker ' num2str(iwalk)])
       -    end
       -    
       -    if i1 == 1
                %xlabel('Interglacial erosion rate [mm/yr]')
                xlabel('Interglacial erosion rate [m/Myr]')
                text(0.02,0.98,'a', 'Units', ...
                    'Normalized', 'VerticalAlignment', 'Top')
       +        epsilon_int_25 = prctile(Ss{iwalk}.ms(i1,:)*1000., 25);
       +        epsilon_int_50 = prctile(Ss{iwalk}.ms(i1,:)*1000., 50);
       +        epsilon_int_75 = prctile(Ss{iwalk}.ms(i1,:)*1000., 75);
       +        
            elseif i1 == 2
                %xlabel('Glacial erosion rate [mm/yr]')
                xlabel('Glacial erosion rate [m/Myr]')
                text(0.02,0.98,'b', 'Units', ...
                    'Normalized', 'VerticalAlignment', 'Top')
       +        epsilon_gla_25 = prctile(Ss{iwalk}.ms(i1,:)*1000., 25);
       +        epsilon_gla_50 = prctile(Ss{iwalk}.ms(i1,:)*1000., 50);
       +        epsilon_gla_75 = prctile(Ss{iwalk}.ms(i1,:)*1000., 75);
       +        
            elseif i1 == 3
                xlabel('Timing of last deglaciation [yr]')
                text(0.02,0.98,'c', 'Units', ...
                    'Normalized', 'VerticalAlignment', 'Top')
       +        
            elseif i1 == 4
                %xlabel('$\delta^{18}$O$_\mathrm{threshold}$ [$^o/_{oo}$]', ...
                    %'Interpreter', 'LaTeX')
                xlabel(['\delta^{18}O_{threshold} [' char(8240) ']'])
                text(0.02,0.98,'d','Units', ...
                    'Normalized', 'VerticalAlignment', 'Top')
       +        record_threshold_25 = prctile(Ss{iwalk}.ms(i1,:), 25);
       +        record_threshold_50 = prctile(Ss{iwalk}.ms(i1,:), 50);
       +        record_threshold_75 = prctile(Ss{iwalk}.ms(i1,:), 75);
       +        
            else
                disp(['Using mname for i1 = ' i1])
                xlabel(fixed_stuff.mname{i1})
       t@@ -694,35 +718,36 @@ html = [html, ...
            '    <tbody>\n'...
            '        <tr>\n'...
            '            <td>&nbsp;</td>\n'...
       -    '            <td>25%%</td>\n'];
       +    '            <td align="center">25%%</td>\n'];
        
        for i=1:Nwalkers
       -    html = [html, '            <td>w', num2str(i), ',25%%</td>\n'];
       +    html = [html, '            <td>',num2str(epsilon_int_25(i)),'</td>\n'];
        end
            
       -html = [html, '            <td>wavg,25%%</td>\n'...
       +html = [html, '            <td>', num2str(sum(epsilon_int_25)/Nwalkers),'</td>\n'...
            '        </tr>\n'...
            '        <tr>\n'...
       -    '            <td>&epsilon;<sub>int</sub></td>\n'...
       -    '            <td>50%%</td>\n'];
       +    '            <td align="center">&epsilon;<sub>int</sub></td>\n'...
       +    '            <td align="center">50%%</td>\n'];
        
        for i=1:Nwalkers
       -    html = [html, '            <td>w', num2str(i), ',50%%</td>\n'];
       +    html = [html, '            <td>',num2str(epsilon_int_50(i)),'</td>\n'];
        end
        
        
       -html = [html, '            <td>wavg,50%%</td>\n'...
       +html = [html, '            <td>', num2str(sum(epsilon_int_50)/Nwalkers),'</td>\n'...
            '        </tr>\n'...
            '        <tr>\n'...
            '            <td>&nbsp;</td>\n'...
       -    '            <td>75%%</td>\n'];
       +    '            <td align="center">75%%</td>\n'];
        
        for i=1:Nwalkers
       -    html = [html, '            <td>w', num2str(i), ',75%%</td>\n'];
       +    html = [html, '            <td>',num2str(epsilon_int_75(i)),'</td>\n'];
        end
        
       -html = [html, '            <td>wavg,75%%</td>\n'...
       +html = [html, '            <td>', num2str(sum(epsilon_int_75)/Nwalkers),'</td>\n'...
            '        </tr>\n'...
       +    '        <tr style="border-bottom:1px solid black"><td colspan="100%"></td></tr>\n'...
            '    </tbody>\n'...
            '</table>\n'...
            ];