Python - Recreate Minitab normal probability plot

Question

Essentially same question as was asked here, but I want to do it in Python. I have used scipy stats to get a probplot, but I want to recreate the confidence interval curves and I'm not sure how to proceed. Can anyone point me in a direction??

This is where I'm at:

This is where I want to be:

Answer 1

I have an answer for the first part of the task but I am not sure how minitab calculates the confidence interval. None of the definitions I found yields something similar. Here is the code for the base plot and the fit:

import numpy as np
import scipy.stats as stats
from matplotlib import scale as mscale
from matplotlib import transforms as mtransforms
from matplotlib.ticker import Formatter, Locator


class PPFScale(mscale.ScaleBase):
    name = 'ppf'

    def __init__(self, axis, **kwargs):
        mscale.ScaleBase.__init__(self)

    def get_transform(self):
        return self.PPFTransform()

    def set_default_locators_and_formatters(self, axis):
        class PercFormatter(Formatter):
            def __call__(self, x, pos=None):
                # \u00b0 : degree symbol
                return "%d %%" % (x*100)

        class PPFLocator(Locator):
            def __call__(self):
                return np.array([1,5,10,20,30,40,50,60,70,80,90,95,99])/100.0

        axis.set_major_locator(PPFLocator())
        axis.set_major_formatter(PercFormatter())
        axis.set_minor_formatter(PercFormatter())

    def limit_range_for_scale(self, vmin, vmax, minpos):
        return max(vmin, 1e-6), min(vmax, 1-1e-6)

    class PPFTransform(mtransforms.Transform):
        input_dims = 1
        output_dims = 1
        is_separable = True

        def ___init__(self, thresh):
            mtransforms.Transform.__init__(self)
            self.thresh = thresh

        def transform_non_affine(self, a):
            out = stats.norm.ppf(a)
            return out


        def inverted(self):
            return PPFScale.IPPFTransform()

    class IPPFTransform(mtransforms.Transform):
        input_dims = 1
        output_dims = 1
        is_separable = True

        def transform_non_affine(self, a):

            return stats.norm.cdf(a)

        def inverted(self):
            return PPFScale.PPFTransform()

mscale.register_scale(PPFScale)


if __name__ == '__main__':
    import matplotlib.pyplot as plt
    from statsmodels.tools.tools import ECDF

    size = 20

    #generate some data
    pf = stats.norm(loc=9, scale=2.0)
    values = pf.rvs(size=size)     
    values.sort()

    #calculate empirical CDF
    cumprob = ECDF(values)(values)

    #fit data
    loc, scale = stats.norm.fit(values)
    pffit = stats.norm(loc=loc,scale=scale)

    x = np.linspace(values.min(),values.max(),3)
    ax = plt.subplot(111)
    ax.plot(values,cumprob, 'go', alpha=0.7, markersize=10)
    ax.plot(x,pffit.cdf(x),'-',label='mean: {:.2f}'.format(loc))
    ax.set_yscale('ppf')
    ax.set_ylim(0.01,0.99)
    ax.grid(True)
    ax.legend(loc=0)
    plt.show()

The nonlinear axis class is taken from one of the matplotlib examples. This gives the following plot: