quantile distribution plot

#sample_dims - spatial dimensions 
#samples_hierarchical_errosion - samples from poisson mean

mu = samples_hierarchical_errosion[f'lm_er_{s}'].mean(axis=0) # [tensor dims locations x genes]
k = iss_data[s] # [observed values locations x genes]
u = stats.uniform.rvs(size=mu.shape[0]*mu.shape[1]).reshape(mu.shape) 
ppois = u * stats.poisson.cdf(k, mu) + (1-u) * stats.poisson.cdf(k-1, mu) # qusntiles of observations given estimated poisson mean  
qnorm = stats.norm.ppf(ppois) * mask[s].reshape(-1,1) #assuming pseudo normal distribution of quantiles 

plt.imshow(qnorm[:,g].reshape(*sample_dims[s]).T[::-1, :], vmin=-3, vmax=3) # plot per gene